# Tag Info

32

Here's an attempt in which I start with a set of "void points", which will be the centres of the gaps between filaments. The stars are then created as an initially random distribution, and are repeatedly nudged away from their nearest void point. Or, to look at it another way, they are attracted towards the edges of the Voronoi cells defined by the void ...

29

This sounds very similar to a traveling salesman problem for which we have the built-in (as of v8) FindShortestPath function. However, this function minimizes the overall path length (which is the sum of all sub-paths), whereas we need a functions that maximizes the path probability (which is the product of the transition probabilities along the path). ...

26

Mathematica v8 does not provide support for automated random number generation from multivariate distributions, specified in terms of its probability density functions, as you have already discovered it. At the Wolfram Technology conference 2011, I gave a presentation "Create Your Own Distribution", where the issue of sampling from custom distribution is ...

25

If you use the code take the hessian of it and plot the map of the largest eigenvalues you get a nice filamentary map like the bottom right panel. see this reference (specifically pp 28 of the phd). In mathematica it can be coded as follows nn = 256;u = GaussianRandomField[nn, 2, Function[k, k^-4]]//GaussianFilter[#, 4] & // Chop; Clear[f]; ...

25

To answer your question: I don't think it's a bad or good idea to use If. It depends on how you do it. To demonstrate I'll use If combined very powerfully with Mathematica 10's ability to tell if a point is inside a specified region or not. step[position_, region_] := Module[{randomStep}, randomStep = RandomChoice[{{-1, 0}, {1, 0}, {0, -1}, {0, 1}}]; ...

23

So your problem reduces to: SeedRandom[8396] RandomReal[NormalDistribution[0, 1], {4, 4, 4}, WorkingPrecision -> 6] It does appear to be caused by reducing WorkingPrecision to 6, because it goes away when you stop forcing Mathematica to behave like a bad pocket calculator. I can't see any reason for you to do this... Better to leave WorkingPrecision ...

21

This is not an efficient answer but it is fun to play with so I thought I'd post it. For efficiency the use of Nearest might provide a good starting point. g[n_, {low_, high_}, minDist_, step_: 1] := Block[{data = RandomReal[{low, high}, {n, 2}], temp, happy, sdata, hdata}, While[True, temp = ((Nearest[data][#, 2][[-1]] & /@ data)); happy ...

21

You are trying to implement Euler-Maruyama simulation method for a 2-stage short-term interest rate model which is given by the following system of SDEs: \begin{eqnarray} \mathrm{d} \theta_t &=& -\lambda_\theta \left( \theta_t - \bar\theta\right) \mathrm{d}t + \sigma_\theta \mathrm{d}W_{\theta,t} \\ \mathrm{d} \pi_t &=& ...

20

Metropolis algorithm Update: ~15x speedup with Compile! I propose an original solution, which consists in using the Metropolis algorithm. It is a very general approach, which is applicable for any probability density function in any dimensions. Metropolis /: Random`DistributionVector[ Metropolis[pdf_, u0_, s_: 1, n0_: 100, chains_: 200], n_Integer, ...

20

The programming style you are using is not very fitting for Mathematica. Here's a better way (shorter, much faster): n = 1000000; (* number of points to use *) octantVolume = N[ Total@UnitStep[1 - Norm /@ RandomReal[1, {n, 3}]]/n ] The reason why you get the error you mention is that for some x, y, the expression 1 - x^2 - y^2 is negative, thus its ...

19

Edit: this answer is now structured in two sections. The first deals about creating a candidate RNG from audio data. The second demonstrates some testing I performed on this RNG. Creating the RNG Okay, I'll got at it another way then. I recorded 10 seconds of ambient noise on my MacBook Pro internal speakers. I was possibly in the worst conditions for ...

19

@Simon Woods points out in a comment that: In fact the delay on the initial run is caused by compiling code to provide the Poisson distribution :-) You can look at ImageColorOperationsDumpiImageEffectPoissonNoise to see how it works internally. Now, although PoissonDistribution can't be compiled, there's nothing stopping the use of my own C++ ...

18

If you want a random vector just because you need some arbitrary vector and you don't really care what it is, then Mr.Wizard's method of picking three random coordinates in [-1,1] will work. But if you care about the statistical properties of your vector, and in particular if you want it drawn from a uniform distribution over the surface of the sphere, then ...

17

I finally found some time to investigate this. I think it warrants a detailed response. In places I will repeat what others have pointed out, but I wanted something that ties together the various threads as best I can discern them. I'm not certain what is meant by the Rule 30 RNG having "an extremely small effective size". Possibly it refers to taking bits ...

17

This answer is going to be a bit of a sprawl. Please read on. I am going to present several methods of simulation, hopefully in increasing order of performance. Method 1 We can carry out the filling of seats, at least as I understand the puzzle, quite literally like this: fillseats[seats_List] := ReplacePart[seats, {{1}, {2}} + RandomChoice @ ...

17

A nice question. Sampling from tCopula is done in stages. First a sample is generated from the copula with uniform marginal distributions, and then quantiles of appropriate marginal distributions are applied to the respective slots. Most of the time goes into evaluation of these quantiles, and they are expensive to compute. Being interested in \$\geqslant ...

17

Lets call your plot res. res = RegionPlot[And @@ Table[ Dot[{Phi1, Phi2}, Eta[b]] <= Norm[NI[Pi] - Eta[b]]^2 + 2, {b, 0, 2 Pi, 2 Pi/10}], {Phi1, -7, 7}, {Phi2, -7, 7}]; Lets extract the mesh Mathematica is generating by default. Use more PlotPoints to get more triangular mesh of your 2D region. pts = res[[1, 1]]; (* Vertices *) {triangles, qd} = ...

17

Since this is rather long, some might prefer a teaser of what is coming: Introduction First of all, I don't really know why you make your figure inconsistent. I mean, from the first big triangle you separate three smaller triangles. Why don't you just repeat this process and inscribe a circle in each of the new triangles and again separate three new ...

16

The mathematica help is very thorough and is very indicative of what you should do next. By way of the histogram diagram obtained, you can compare your data against the proposed distribution. Show[Histogram[w[[2, 1]], Automatic, "ProbabilityDensity"], Plot[PDF[h["FittedDistribution"], x], {x, 0, 1500}, PlotStyle -> Thick]] The reference points you ...

15

I don't have much time right now to write a longer post (and anyway, this is not really an area in which I would trust myself as an "expert") but I think it is better not to use the words "high quality" and "low quality" too loosely. Essentially, when dealing with quasi-random number generators there are at least two senses in which these terms are used, and ...

15

David's answer has given the methods for producing random points that are uniformly distributed over the surface of the sphere. Of course, there are other probability distributions on the sphere that are of interest, as well as a number of methods for generating them. For instance, here is how to generate a random unit vector which follows the von ...

15

May be the most compact approach: Hold[i1++, i1--, i2++, i2--][[RandomInteger[{1, 4}]]]

15

When n is large it's much faster to operate on a 3 x n array than to process each of the n 3-vectors separately. This is one of the standard "tricks" to speed things up. n = 10^6; (* Isn't that easier to read than 1000000 ? *) AbsoluteTiming @ N[ Total@UnitStep[ 1. - Norm/@RandomReal[1,{n,3}] ]/n ] (* {4.555842, 0.524302} *) AbsoluteTiming @ N[ ...

15

Method of random number generation is also significant: Default: n = 10^6; AbsoluteTiming[N@Mean@UnitStep[1. - Total[RandomReal[1, {3, n}]^2]] - π/6] {0.197896, 0.000649224} Niederreiter low-discrepancy sequence (see "methods" here): SeedRandom[Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 3}}]; ...

15

A simple numerical maximization using NMaximize as suggested by b.gatessucks: pts = Array[{x[#], y[#]} &, 10]; mindist2 = Min[#.# & /@ Subtract @@@ Subsets[pts, {2}]]; vars = Flatten[pts]; constraints = Thread[0 <= vars <= 1]; {md2, rules} = NMaximize[{mindist2, constraints}, vars]; minimaldistance = Sqrt[md2] (* 0.381759 *) ...

14

If you don't need high precision, you can do something along these lines: canvas = Image@ConstantArray[0, {100, 100}]; distance = 6; {img, {pts}} = Reap[Nest[ ImageCompose[#, SetAlphaChannel[#, #] &@Image@DiskMatrix[distance], Sow@RandomChoice@ Position[Transpose@ImageData[#, DataReversed -> True], 0.]] &, ...

14

In general you should use RandomVariate for distributions and RandomReal for uniforms. Often RandomVariate calls RandomReal or RandomInteger under the hood but it varies on a distribution by distribution basis. After loading any necessary symbols, on evaluation, any timing differences should be negligible. RandomVariate is intended to give the flexibility ...

14

To my knowledge UniformDistribution and NormalDistribution are the only distributions that are directly compilable for RandomVariate. Consider that sampling from a UniformDistribution is what RandomReal was originally designed to do. This code is likely written deep down in C and so compiles without any special effort. In order to hook up RandomVariate ...

14

Here's a reorganization of GaussianRandomField[] that works for any valid dimension, without the use of casework: GaussianRandomField[size : (_Integer?Positive) : 256, dim : (_Integer?Positive) : 2, Pk_: Function[k, k^-3]] := Module[{Pkn, fftIndgen, noise, amplitude, s2}, Pkn = Compile[{{vec, _Real, 1}}, With[{nrm = Norm[vec]}, ...

14

As ssch alluded to in a comment this can be easily done with NestWhileList: NestWhileList[# + RandomChoice[{-1, 1}] &, 0, # < 10 &] This is hardly the only way to proceed however. You could for example Sow and Reap values from within a recursion: f[x_] /; x < 10 := f @ Sow[x + RandomChoice[{-1, 1}]]; Reap[f[0]][[2, 1]] This actually ...

Only top voted, non community-wiki answers of a minimum length are eligible