26 ms·

Generating random points inside a sphere

252 points by kkaranth 8y ago

acidburnNSA 8y ago
Correct uniform sampling in things like this comes up a lot in Monte Carlo radiation transport problems and is covered very well in your basic first year graduate Monte Carlo class in a good nuclear engineering program. If you want to go deeper, check out the literature around that.
- chrisseaton 8y ago
  > in your basic first year graduate Monte Carlo class How many years of graduate classes are there on Monte Carlo methods for there to be a 'first year' one?!
  - vortico 8y ago
    I think he/she means the Monte Carlo class you take in your first year of nuclear grad courses.
    - acidburnNSA 8y ago
      Yeah, generally there's a lot of coursework like this in the first year. Then if you specialize in it you're going to take like 2-3 advanced classes plus relevant math and computer science stuff.
  - cjslep 8y ago
    Depends if your graduate focus is on core neutronics or core thermodynamics.
- banachtarski 8y ago
  Monte carlo based sampling shows up also in finance, graphics programming, materials science, etc. Graphics programming and physics were areas where uniform sampling of a geometrical surfaces were important for me.
yiyus 8y ago
This problem is equivalent to generating random orientations, which can be easily done generating random unit quaternions (see, for example, [1]). Each quaternion found will correspond to a point in the surface of a 4D sphere, which can be mapped to a point inside a 3D sphere. [1] http://mathproofs.blogspot.com/2005/05/uniformly-distributed-random-unit.html http://mathproofs.blogspot.com/2005/05/uniformly-distributed...
- jacobolus 8y ago
  If you take the stereographic projection of a uniform random unit quaternion onto 3-space, then to get a uniform distribution in the ball you need to apply a non-linear scaling function to the radius. I’m not sure this really makes much sense as a way to generate random points in the ball (especially if you care about the details of rounding etc.), but maybe fun for someone to implement. I would expect the rejection method to be noticeably faster. * * * Another cute thing you could do is generate 3 normally distributed random numbers (giving a normally distributed random point in a 3-space), and then apply the appropriate non-linear scaling to the radius.
tambourine_man 8y ago
Pleasantly surprised by the effect when scrolling over the sphere
- pilaf 8y ago
  I couldn't see any spheres at all so I was confused by your comment, which led me to realize WebGL wasn't working, which in turn reminded me I had disabled hardware acceleration in Chrome long ago because of some bug which made browsing unbearably slow. I've now re-enabled it and it seems the bug has since been fixed. So, thank you!
  - dzhiurgis 8y ago
    It’s the forced graphics card switching that causes system crash in every 20th case.
- undersuit 8y ago
  Same, you can also drag them around to a new rotation.
alstange 8y ago
So which algorithm ended up being fastest?
- jacobolus 8y ago
  Depends on how fast your PRNG is, but I would expect the rejection method to be best under pretty much any circumstances. The code is also simplest. If PRNG speed becomes a bottleneck just use a faster PRNG.
- Cyph0n 8y ago
  To my untrained eye, it looks like the first naive algorithm would be fastest. Two of them require trig functions while the other relies on sqrt and cube root, all of which are expensive functions when compared to the first approach. So I'd wager that the 50% reject rate would still be cheaper at scale than evaluating trigonometric functions.
  - Buge 8y ago
    But you might want to consider average case time and worst case time. The first algorithm might have a good average case time, but its worst case time is infinity, which might be unacceptable.
    - the8472 8y ago
      The probability of experiencing bad cases decays exponentially with the number of iterations. 0.48^n. Your CPU would spontaneously undergo fusion by tunneling before hitting on a case where it's stuck on a bad case for more than 10 microseconds.
      - glhaynes 8y ago
        Could you say a little more about what relevance 0.48 has here?
        mrastro 8y ago
andrew3726 8y ago
This is also very important in many global illumination techniques which requires sampling the hemisphere around a point. As case in point a simple Monte Carlo path tracer uses this to sample random directions, see e.g. [0]. In addition to this there are many 'tricks' one can use to reduce the variance/noise using a differently weighted sampling (to help the converge of MC). A nice technique to get a cosine-weighted hemisphere sampling is to generate a random point on the unit disc and then project them onto the (hemi)sphere. Seems somewhat relevant here. [0] https://en.wikipedia.org/wiki/Path_tracing#Algorithm https://en.wikipedia.org/wiki/Path_tracing#Algorithm
machiavelli1024 8y ago
Notice how the rejection method becomes useless for higher dimensions. https://en.m.wikipedia.org/wiki/Curse_of_dimensionality https://en.m.wikipedia.org/wiki/Curse_of_dimensionality
- peteretep 8y ago
  While that's interesting, I'm having a complete lack of imagination fail on why one would be trying to solve this problem for higher dimensions, other than academic pursuits. Could you give a practical example?
  - machiavelli1024 8y ago
    “Sampling from the uniform distribution on the N-dimensional Euclidean ball (the N-ball) and its surface (the N-sphere) is a tool useful in many diverse research fields, for instance Monte Carlo integration in physics, generating random directions for MCMC sampling on bounded regions, generating random correlation matrices, Monte Carlo analysis of random packing on the sphere, generating random rotations in cryptography as well as various simulation studies in statistics” https://www.sciencedirect.com/science/article/pii/S0047259X10001211 https://www.sciencedirect.com/science/article/pii/S0047259X1...
DanAndersen 8y ago
Related: selecting points on the surface of a sphere https://www.jasondavies.com/maps/random-points/ https://www.jasondavies.com/maps/random-points/ Poisson disc sampling can be quite useful.
- 2dvisio 8y ago
  I guess what it would look like in 3D :)
mwkaufma 8y ago
"about 48% of the points chosen are discarded. This may not be a huge problem if the number of points that need to be generated is small" ... or if the number of points isn't small -- it's only 2x the cost of the simple linear case without a loss of uniformity or adding any costly trig functions. That's still pretty good!
- chris_mc 8y ago
  2x the compute resources is still horrible for our environment, which isn't pretty good to me. Computer Scientists/Software Engineers need to get it into their head that wasting compute resources is directly increasing our footprint on this planet. We need to simplify the web not because of bandwidth use, but because all that compute resource moving around useless ads and scripts is largely a waste. We need to stop using block-chain for everything where a centralized ledger would suffice because block-chain is literally wasting resources in those problems for no reason. We need to start thinking about our software's impact on the environment due to wasted compute cycles for things that are unnecessary or which could be improved by reducing the cycles needed and stop doing that lazy or fancy stuff.
  - jpatokal 8y ago
    If you read the comments, you'll see that the "wasteful" method is actually the cheapest computationally.
    - chris_mc 8y ago
      I'm responding to the general lack of concern by the commenter above and most engineers I know. You and this other guy commenting are missing the point, I'm not talking specifically about this algorithm, but about engineers in general who don't care about wasted cycles.
      - tomc1985 8y ago
        Yes but that's the entire point of these comments, wasting cycles is a great optimization but unless the goal is purely intellectual exercise, you lose more time with expensive trig functions. The real world favors the naive approach, in this context + case.
woopwoop 8y ago
A note about the rejection method: in high dimensions, this becomes very expensive. This is because the volume of a ball inscribed in a unit cube goes to zero (very quickly) as the dimension goes to infinity, so you start rejecting points with high probability.
- Pxtl 8y ago
  I think that might be premature optimization to worry about finding random points in an n-dimensional hypersphere where n is a huge number.
  - woopwoop 8y ago
    Sure, but if you can't prematurely optimize in HN comments, when can you?
  - tantalor 8y ago
    No, n does not need to be huge, n=8 is high enough.
    - Retric 8y ago
      In this context 8 is a huge number.
      - cultus 8y ago
        Not at all. There are plenty of circumstances in all kinds of numerics where one would want points on the sphere in high dimensions. It's not just a geometric object that has limited practical value for d>3. It is the set of points where the norm is unity, which is a pretty fundamental concept in all sorts of places.
        Retric 8y ago
        Not as part of a graphic's library. You might want to do 4d or 5d. But past that is very unlikely. In terms of actually solving to proble you just want a vector random distance and a function that to map random distance to actual distance. But if you never deal in 6+ dimensions that's a waste of time.
vzaliva 8y ago
When the author mentions "normally distributed numbers" it sounds like he is sampling a normal (Gaussian) distribution. In fact, he is sampling from a uniform distribution.
- kkylin 8y ago
  I think he does mean sampling from a gaussian. You can tell because the example calls randn(). The gaussian distribution is invariant under rotations in 3d (and more generally invariant under any orthogonal transformation in any dimension), so if one takes a random vector X whose components are independent N(0,1) random variables and compute X/|X| (where |X| is the usual euclidean norm of X) the result is guaranteed to be uniform on the unit sphere.
  - kkaranth 8y ago
    Indeed, I've used randn() wherever a normally distributed random number is needed. This is computed using the Box-Muller method.
charmides 8y ago
The graphics are nice, but this article really could have used some complexity analysis and some probability theory. The author neither discusses the asymptomatic complexity of the algorithms, nor the run-time of his implementation (which is more pertinent here), nor gives any proofs of why some of the algorithms sample uniformly from the unit ball and why some of them don't. Also, it would have been really nice to generalize this problem to n dimensions. I assume that there is some small value of n where naive rejection sampling is worse in practice than one of the more sophisticated methods.
- jpatokal 8y ago
  Does he need to? All algorithms given except the first appear to be O(N), and their performance is going to depend almost entirely on how fast the local implementations of trig operations involved are.
  - gamegoblin 8y ago
    Isn’t the first method still O(N)? You will reject a constant percentage of points on average (the author says 48%), but constants fall out of big-O.
    - jpatokal 8y ago
      That's true if you assume a while pending on rand() is constant, but that wasn't covered in my CS 101 on complexity analysis...
      - Sharlin 8y ago
        Average-case constant. Worst case? Never halts. Generating a single point doesn’t depend at all on our choice of n (the number of points generated). But yeah, complexity analysis of stochastic algorithms is fairly advanced stuff.
        heavenlyblue 8y ago
nayuki 8y ago
I think the section "Using normally distributed random numbers" deserves an additional explanation. The code may look simple, but like most of the other methods mentioned on the page, generating normally (Gaussian) distributed random numbers still requires the use of rejection sampling and sine/cosine functions. See https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
- strainer 8y ago
  I implemented Marsaglia's Box Muller method in a library (^) and it seems quite quick. It doesn't involve trig functions, just a sqrt, a log, and an average of about two rands per normal outputted. It seems to work out around the same speed as a trig function - in javascript. For some purposes just averaging a few rands could be used to get a good enough approximation of a normal. I wonder if a short Taylor series could be mined that could quite accurately map equal to gaussian distribution. [^]http://strainer.github.io/Fdrandom.js/ http://strainer.github.io/Fdrandom.js/
nalourie 8y ago
For people who like the simplicity of the rejection algorithm, it might be worth checking out Ziggurat sampling which is based on a similar idea, but shapes the sampling region to reduce rejections: https://en.m.wikipedia.org/wiki/Ziggurat_algorithm https://en.m.wikipedia.org/wiki/Ziggurat_algorithm
- charmides 8y ago
  To add to that comment, importance sampling is also an attempt to reduce the number of rejections and could be interesting here. https://en.m.wikipedia.org/wiki/Importance_sampling https://en.m.wikipedia.org/wiki/Importance_sampling
  - throwaway313213 8y ago
    Importance sampling works when you're interested in some linear function[al] of the distribution; not for actual sampling.
cjhanks 8y ago
Maybe a mathematician can correct me. But could you not also model 3D vertices of uniform density into a 1D space and then uniformly sample the 1D space. Of course this may have issues with highly dense objects or very large ones.
- lixtra 8y ago
  I don’t get what you mean by 1D spaces but you could triangulate into tetrahedra and sample those[1] (taking their volume into account). You will write a lot of code and still have the triangulation approximation error. For high dimensional spheres you can possibly subdivide the space in hyper cubes and do the in/out approach inside each cube that contains a part of the sphere to get a better hit miss ratio. [1] http://vcg.isti.cnr.it/jgt/tetra.htm http://vcg.isti.cnr.it/jgt/tetra.htm
  - cjhanks 8y ago
    Yes, I was thinking you could triangulate the sphere at discreet shells. I have done 2D symmetric polygons where you needed to uniquely index into a 1D array of vertex values. Ie, based on the magnitude you can estimate which shell. By knowing the shell, you know the number of vertices for the shell, by knowing both you can map it to a specific vertex, as you can define a priori the absolute traversal order. I have no idea if the 2D case extends to higher dimensions.
herbertkane 8y ago
Section 3 of this paper explains this topic pretty well in N-dimensions: http://jmlr.org/papers/volume17/blaser16a/blaser16a.pdf http://jmlr.org/papers/volume17/blaser16a/blaser16a.pdf To get a random point within the sphere (rather than on the surface of the sphere, as described in the paper) you also need to pick a random distance from the origin to randomly scale the point on the surface.
throwaway313213 8y ago
This is similar to the method used for generating random numbers using Box-Muller transform. https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
arnarbi 8y ago
The main issue with the second method, which isn't discussed at all, is that it's more likely to pick a direction towards the corners of the unit cube. You can even see that bias in the point cloud. The problem with choosing d uniformly is also better explained by pointing out that shells of larger radius have larger surface, thus they'll be underrepresented.
sgentle 8y ago
I thought this was kinda fun, and I was trying to figure out why you couldn't just take a random x and use it to find a y that preserves x²+y²=r² (ie by setting y to a random number between 0 and ±√(r²-x²)). In the end I figured out that you need to generate x the same way, but if you have, say: n=rand(r), x=rand(±√(r²-n²)), y=rand(±√(r²-x²)) it seems to work. Here's a version I extended to three dimensions, designed to be run in the JS console of the original article (it'll show up underneath the first sphere): new SphereSimulation(document.getElementById("spheres1"), () => { const rand = Math.random const root = (x) => rand() < 0.5 ? Math.sqrt(x) : -Math.sqrt(x) const invdist = (a=0, b=0) => root(RADIUS*RADIUS - a*a - b*b) const x1 = rand() * invdist() const x2 = rand() * invdist(x1) const y2 = rand() * invdist(x2) const x3 = rand() * invdist(x2, y2) const y3 = rand() * invdist(x3, y2) const z3 = rand() * invdist(x3, y3) return new THREE.Vector3(x3, y3, z3) }) It looks okay to me - am I missing something? Is this a reasonable approach?
- goldenkey 8y ago
  I posted a comment where I described this method as my preferred solution. You did not need to patch it up though..its fine to just generate the 3 numbers and then choose a random assignment from them to x y and z. This fixes the bias that would be present if you always generated x first.
- jpeanuts 8y ago
  Your missing the requirement that the random sample be uniform in the sphere. That is, if the sphere is cut into any number of pieces of equal volume, the sample is equally likely to appear in any piece. Your approach samples uniformly from the radius, so you're going to end up with far too many samples close to the center.
  - sgentle 8y ago
    Interesting... thanks for the pointer. I did some testing and the points definitely do seem to cluster along the axis... just increasing the number of iterations of √(r²-n²) seemed to improve it though. At 5, I can't tell the difference anymore: https://codepen.io/anon/pen/xJJOEq https://codepen.io/anon/pen/xJJOEq
goldenkey 8y ago
The method I immediately think of is just based on the relation. Sqrt(x^2 + y^2 +z^2) <= 1 reduces to x^2 + y^2 +z^2 <= 1 Generate the random value for x^2, any random value <= 1 Now we have y^2 +z^2 <= 1-x^2 Generate a random value for y^2, any random value <= 1-x^2 Finally generate a random value for z^2, any random value <= 1-x^2-y^2 Take the square roots of each number we generated to get the final coordinates. But wait! Theres another step..the most important one! There are 6 orderings of (x,y,z) We must choose a random one because otherwise x will be biased towards higher values compared to y and z. So I kind of lied. All we needed was 3 random values via the subtraction methods above. Then we need to randomly choose their ordering as (#,#,#) This seems like the easiest most trivial solution that would be uniformly distributed. Let me know if I am wrong.
- kgwgk 8y ago
  You are wrong.
- machiavelli1024 8y ago
  >Let me know if I am wrong. You're wrong on multiple levels. 1. Generating random values for x^2, y^2, z^2 and taking their square root will only give you values in the x,y,z > 0 octant. (But let's say you "fix" this by randomly multiplying them with -1.) 2. Taking the square root of a uniformly distributed random variable is no longer uniformly distributed. 3. Randomly reordering the coordinates won't fix your bias. Here's a demonstration in 2D: https://jsfiddle.net/tz85wnxy/59/ https://jsfiddle.net/tz85wnxy/59/ Uncomment line 17, 18, 20 to see how it's still not uniform even if you randomly multiply the coordinates by -1 and reorder them.
  - goldenkey 8y ago
    Works fine as long as you generate x,y,z randomly and not their squares. Same process but minor tweak. https://jsfiddle.net/w8zLvsy2/ https://jsfiddle.net/w8zLvsy2/
    - machiavelli1024 8y ago
enugu 8y ago
For the second example, sampling the radius uniformly leads to an error with a higher point density in the centre of sphere. The radius, r, has a weight proportional to r^2 or the surface area. So, one can take a uniform random sample, s between 0 and R^2, and take the square root of s. That will give higher probability density around larger radius. Then, the second picture should look like the first one.
- nwellnhof 8y ago
  Seems like you stopped reading after the second example.
  - enugu 8y ago
    Yes, didnt look at the details of the next methods. He adjusts for this in the next examples by taking the cube root of the uniform sample(not square root as in my post), which is the correct adjusting factor, as a small shell will have weight in proportion to difference between cubic powers, which is proportional to r^2.
- deleted 8y ago
  [deleted]
anewhnaccount2 8y ago
Similar problem which ends up not being quite as simple in the general case as you might think: sampling from a unit simplex: http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf
jjgreen 8y ago
http://blog.geomblog.org/2013/01/a-sampling-gem-sampling-from-ellp-balls.html http://blog.geomblog.org/2013/01/a-sampling-gem-sampling-fro...
expuexto 8y ago
I don't agree with the way he normalizes the radius in the last 2 methods. I think it should be a square root, not a cubic one. This is because the volume diferential r^2•sin(phi)•d_phi•d_theta•dr is proportional to the square of r. This mistake is also made in a source he cites: https://math.stackexchange.com/questions/87230/picking-random-points-in-the-volume-of-sphere-with-uniform-probability/2872878 https://math.stackexchange.com/questions/87230/picking-rando... and I tried to point it out there as well.
- improbable22 8y ago
  That's the correct volume element, but not what you want. The volume inside radius r grows as r^3. It's the inverse of this _cumulative_ density function which maps uniform-in-r numbers to uniform-in-volume ones.
- antidesitter 8y ago
  Did you empirically check your hypothesis? What's the integral of r^2 dr? Think about whether your claim makes sense when n = 1.
- expuexto 8y ago
  You both are right. Thank you for the comments.
plg 8y ago
That (last method) is a lot of sines and cosines, cube root and an acos... I bet for 3D the rejection method is not that much slower, in wall clock time. OK so I coded up the first one and the last one in C and compiled on my mac using whatever gcc links to. I found 1e7 points in the sphere: first routine: 0.521s last routine: 0.850s [edit]: C code here: https://pastebin.com/zH0BehSM https://pastebin.com/zH0BehSM
- froogie 8y ago
  The first routine is about 3 rands per point minus rate of discarded, about half? So 6 rands per point. Now, since the second one has three rands in it, and plenty of trig, it really becomes a question if you can make the trigonometry faster than three rands. But as you saw, straightforward attempt will pretty much certainly end up having worse performance. But what if one used optimizations like lookup tables and trigonometric approximations instead? Could they end up costing less cycles than three rands? Maybe not LUTs, maybe other trickery someone knows of? And what about the cost of those three rands? Could one use a fast xorshift with good enough results here? Does one even need thee random vectors? Maybe two + modulo trickery will do?
  - acqq 8y ago
    Typically the generation of the uniform random numbers can be made quite cheaply, especially if the resolution is not critical (e.g. if uniform floats are enough and not uniform doubles etc) or if the quality of the randomness is not critical, and in these cases discarding half of these will win versus other costly operations.
- improbable22 8y ago
  Note that his last way does cos(acos(cosϕ))) after generating cosϕ uniformly. Which is a waste, you can just set sinϕ = sqrt(1-cosϕ^2). After fixing that, this way is faster for me. Some trig, but no branches. Edit: https://pastebin.com/UExZ5Siw https://pastebin.com/UExZ5Siw
  - shpx 8y ago
hatsunearu 8y ago
So what ends up being the general purpose method for uniformly distributed random numbers in some arbitrary shape? (other than the rejection method of course) If arbitrary shapes always requires the rejection method, what are some "classes" that have easier algorithms (the sphere is obviously included here)
- smaddox 8y ago
  Not my area of expertise, but based on what I've read, you general try to get as close as possible through pure mappings, and then use rejection sampling to throw out the samples that don't fall under the desired distribution. This can work even if the distribution is dynamically changing.
wruza 8y ago
Interestingly, it seems to be all about bending some distribution and then unbending it again. Orthogonal dimensions are independent, so choosing three random numbers generates uniform in-cube distribution. But once you bend it via polar coordinates or vectors, you “wrap” many points of polar/vector abstract space to near the cartesian origin. And then you have to bend back, to make it uniform again. That raises interesting question on whether it relates to topology problems. Since generating N randoms is a fastest primitive operation, is there a quick and uniform translation from a cube to the sphere, such that all points are still uniform? Can we apply space-filling curve magic here (like Hilbert’s)? I feel like there is much more geometry power under that, but sadly I’m just a guy with no math bg.
- 3pt14159 8y ago
  It's an interesting idea, but I think it would prove to be equivalent. As another commenter pointed out, the efficiency of this operation really only matters in higher dimensions. I bet a generative model would work, though it's tricky enough for me not to be able to think it through on the spot. Something like: 1. Pick a subset of points thus far generated at random. 2. From this subset, pick the point with the furthest distance from the origin. 3. Form a new unit sphere where r is 1-d of the picked point (ie, make a new sphere that is a subset of the parent sphere with a single touching point of the parent sphere). 4. Generate a new point in this subsphere using the normalized vector method. Not exactly the above, but something like it. The downside of the above approach is that I think it will create a distribution of points that may be unnaturally uniform. It would be interesting to define a set of equations that could benchmark different approaches for how well they compare to the in-cube strategy.
chias 8y ago
Interestingly, the method in which you determine your "random point" can greatly affect the end result, so swapping out the methods (e.g. "generate three coordinates and discard" method for the "generate a direction and a length" method) may produce wildly different results. An interesting illustration of this is Bertrand's Paradox, which asks the question: "if you draw a random chord inside a circle, what is the probability that its length is greater than sqrt(3)?". The correct answer is 1/2. The correct answer is also 1/3. Oh, and the correct answer is also 1/4. It depends on three equally valid (and uniformly distributed) ways of defining your random chord. See: http://web.mit.edu/tee/www/bertrand/ http://web.mit.edu/tee/www/bertrand/
- improbable22 8y ago
  Paradox seems like a strong word for this. "Draw a random chord" is almost a pun, draw from what distribution? Until you tell me, there really is no one answer.
  - chias 8y ago
    > "Draw a random chord" is almost a pun No more so than 'pick a random point in a circle'. Note in all cases we're talking the uniform distribution. The question of using Cartesian coordinates, parametric coordinates, etc. is what causes the differences.
    - improbable22 8y ago
      I'm complaining that there isn't one uniform distribution. For example uniform in x on [1, 2] is different to uniform in 1/x on the same interval, and which one you draw from may matter. You have to be told. In fact every (one-variable) distribution is uniform in something! "pick a random point inside a circle" comes with an implicit "uniformly according to the ordinary euclidean metric on the paper". I guess it's a judgement call that this is obvious enough, but "a random chord" is clearly not... otherwise we wouldn't be here!
- deleted 8y ago
  [deleted]
Myrmornis 8y ago
I find it unclear whether the author is aiming for Uniform, or just “looks pretty uniform”.

252 points by kkaranth 8y ago

Correct uniform sampling in things like this comes up a lot in Monte Carlo radiation transport problems and is covered very well in your basic first year graduate Monte Carlo class in a good nuclear engineering program. If you want to go deeper, check out the literature around that.

chrisseaton 8y ago
> in your basic first year graduate Monte Carlo class How many years of graduate classes are there on Monte Carlo methods for there to be a 'first year' one?!
- vortico 8y ago
  I think he/she means the Monte Carlo class you take in your first year of nuclear grad courses.
  - acidburnNSA 8y ago
    Yeah, generally there's a lot of coursework like this in the first year. Then if you specialize in it you're going to take like 2-3 advanced classes plus relevant math and computer science stuff.
- cjslep 8y ago
  Depends if your graduate focus is on core neutronics or core thermodynamics.
banachtarski 8y ago
Monte carlo based sampling shows up also in finance, graphics programming, materials science, etc. Graphics programming and physics were areas where uniform sampling of a geometrical surfaces were important for me.

yiyus 8y ago

This problem is equivalent to generating random orientations, which can be easily done generating random unit quaternions (see, for example, [1]). Each quaternion found will correspond to a point in the surface of a 4D sphere, which can be mapped to a point inside a 3D sphere. [1] http://mathproofs.blogspot.com/2005/05/uniformly-distributed-random-unit.html http://mathproofs.blogspot.com/2005/05/uniformly-distributed...

jacobolus 8y ago
If you take the stereographic projection of a uniform random unit quaternion onto 3-space, then to get a uniform distribution in the ball you need to apply a non-linear scaling function to the radius. I’m not sure this really makes much sense as a way to generate random points in the ball (especially if you care about the details of rounding etc.), but maybe fun for someone to implement. I would expect the rejection method to be noticeably faster. * * * Another cute thing you could do is generate 3 normally distributed random numbers (giving a normally distributed random point in a 3-space), and then apply the appropriate non-linear scaling to the radius.

tambourine_man 8y ago

Pleasantly surprised by the effect when scrolling over the sphere

pilaf 8y ago
I couldn't see any spheres at all so I was confused by your comment, which led me to realize WebGL wasn't working, which in turn reminded me I had disabled hardware acceleration in Chrome long ago because of some bug which made browsing unbearably slow. I've now re-enabled it and it seems the bug has since been fixed. So, thank you!
- dzhiurgis 8y ago
  It’s the forced graphics card switching that causes system crash in every 20th case.
undersuit 8y ago
Same, you can also drag them around to a new rotation.

alstange 8y ago

So which algorithm ended up being fastest?

jacobolus 8y ago
Depends on how fast your PRNG is, but I would expect the rejection method to be best under pretty much any circumstances. The code is also simplest. If PRNG speed becomes a bottleneck just use a faster PRNG.
Cyph0n 8y ago
To my untrained eye, it looks like the first naive algorithm would be fastest. Two of them require trig functions while the other relies on sqrt and cube root, all of which are expensive functions when compared to the first approach. So I'd wager that the 50% reject rate would still be cheaper at scale than evaluating trigonometric functions.
- Buge 8y ago
  But you might want to consider average case time and worst case time. The first algorithm might have a good average case time, but its worst case time is infinity, which might be unacceptable.
  - the8472 8y ago
    The probability of experiencing bad cases decays exponentially with the number of iterations. 0.48^n. Your CPU would spontaneously undergo fusion by tunneling before hitting on a case where it's stuck on a bad case for more than 10 microseconds.
    - glhaynes 8y ago
      Could you say a little more about what relevance 0.48 has here?
      - mrastro 8y ago

andrew3726 8y ago

This is also very important in many global illumination techniques which requires sampling the hemisphere around a point. As case in point a simple Monte Carlo path tracer uses this to sample random directions, see e.g. [0]. In addition to this there are many 'tricks' one can use to reduce the variance/noise using a differently weighted sampling (to help the converge of MC). A nice technique to get a cosine-weighted hemisphere sampling is to generate a random point on the unit disc and then project them onto the (hemi)sphere. Seems somewhat relevant here. [0] https://en.wikipedia.org/wiki/Path_tracing#Algorithm https://en.wikipedia.org/wiki/Path_tracing#Algorithm

machiavelli1024 8y ago

Notice how the rejection method becomes useless for higher dimensions. https://en.m.wikipedia.org/wiki/Curse_of_dimensionality https://en.m.wikipedia.org/wiki/Curse_of_dimensionality

peteretep 8y ago
While that's interesting, I'm having a complete lack of imagination fail on why one would be trying to solve this problem for higher dimensions, other than academic pursuits. Could you give a practical example?
- machiavelli1024 8y ago
  “Sampling from the uniform distribution on the N-dimensional Euclidean ball (the N-ball) and its surface (the N-sphere) is a tool useful in many diverse research fields, for instance Monte Carlo integration in physics, generating random directions for MCMC sampling on bounded regions, generating random correlation matrices, Monte Carlo analysis of random packing on the sphere, generating random rotations in cryptography as well as various simulation studies in statistics” https://www.sciencedirect.com/science/article/pii/S0047259X10001211 https://www.sciencedirect.com/science/article/pii/S0047259X1...

DanAndersen 8y ago

Related: selecting points on the surface of a sphere https://www.jasondavies.com/maps/random-points/ https://www.jasondavies.com/maps/random-points/ Poisson disc sampling can be quite useful.

2dvisio 8y ago
I guess what it would look like in 3D :)

mwkaufma 8y ago

"about 48% of the points chosen are discarded. This may not be a huge problem if the number of points that need to be generated is small" ... or if the number of points isn't small -- it's only 2x the cost of the simple linear case without a loss of uniformity or adding any costly trig functions. That's still pretty good!

chris_mc 8y ago
2x the compute resources is still horrible for our environment, which isn't pretty good to me. Computer Scientists/Software Engineers need to get it into their head that wasting compute resources is directly increasing our footprint on this planet. We need to simplify the web not because of bandwidth use, but because all that compute resource moving around useless ads and scripts is largely a waste. We need to stop using block-chain for everything where a centralized ledger would suffice because block-chain is literally wasting resources in those problems for no reason. We need to start thinking about our software's impact on the environment due to wasted compute cycles for things that are unnecessary or which could be improved by reducing the cycles needed and stop doing that lazy or fancy stuff.
- jpatokal 8y ago
  If you read the comments, you'll see that the "wasteful" method is actually the cheapest computationally.
  - chris_mc 8y ago
    I'm responding to the general lack of concern by the commenter above and most engineers I know. You and this other guy commenting are missing the point, I'm not talking specifically about this algorithm, but about engineers in general who don't care about wasted cycles.
    - tomc1985 8y ago
      Yes but that's the entire point of these comments, wasting cycles is a great optimization but unless the goal is purely intellectual exercise, you lose more time with expensive trig functions. The real world favors the naive approach, in this context + case.

woopwoop 8y ago

A note about the rejection method: in high dimensions, this becomes very expensive. This is because the volume of a ball inscribed in a unit cube goes to zero (very quickly) as the dimension goes to infinity, so you start rejecting points with high probability.

Pxtl 8y ago
I think that might be premature optimization to worry about finding random points in an n-dimensional hypersphere where n is a huge number.
- woopwoop 8y ago
  Sure, but if you can't prematurely optimize in HN comments, when can you?
- tantalor 8y ago
  No, n does not need to be huge, n=8 is high enough.
  - Retric 8y ago
    In this context 8 is a huge number.
    - cultus 8y ago
      Not at all. There are plenty of circumstances in all kinds of numerics where one would want points on the sphere in high dimensions. It's not just a geometric object that has limited practical value for d>3. It is the set of points where the norm is unity, which is a pretty fundamental concept in all sorts of places.
      - Retric 8y ago
        Not as part of a graphic's library. You might want to do 4d or 5d. But past that is very unlikely. In terms of actually solving to proble you just want a vector random distance and a function that to map random distance to actual distance. But if you never deal in 6+ dimensions that's a waste of time.

vzaliva 8y ago

When the author mentions "normally distributed numbers" it sounds like he is sampling a normal (Gaussian) distribution. In fact, he is sampling from a uniform distribution.

kkylin 8y ago
I think he does mean sampling from a gaussian. You can tell because the example calls randn(). The gaussian distribution is invariant under rotations in 3d (and more generally invariant under any orthogonal transformation in any dimension), so if one takes a random vector X whose components are independent N(0,1) random variables and compute X/|X| (where |X| is the usual euclidean norm of X) the result is guaranteed to be uniform on the unit sphere.
- kkaranth 8y ago
  Indeed, I've used randn() wherever a normally distributed random number is needed. This is computed using the Box-Muller method.

charmides 8y ago

The graphics are nice, but this article really could have used some complexity analysis and some probability theory. The author neither discusses the asymptomatic complexity of the algorithms, nor the run-time of his implementation (which is more pertinent here), nor gives any proofs of why some of the algorithms sample uniformly from the unit ball and why some of them don't. Also, it would have been really nice to generalize this problem to n dimensions. I assume that there is some small value of n where naive rejection sampling is worse in practice than one of the more sophisticated methods.

jpatokal 8y ago
Does he need to? All algorithms given except the first appear to be O(N), and their performance is going to depend almost entirely on how fast the local implementations of trig operations involved are.
- gamegoblin 8y ago
  Isn’t the first method still O(N)? You will reject a constant percentage of points on average (the author says 48%), but constants fall out of big-O.
  - jpatokal 8y ago
    That's true if you assume a while pending on rand() is constant, but that wasn't covered in my CS 101 on complexity analysis...
    - Sharlin 8y ago
      Average-case constant. Worst case? Never halts. Generating a single point doesn’t depend at all on our choice of n (the number of points generated). But yeah, complexity analysis of stochastic algorithms is fairly advanced stuff.
      - heavenlyblue 8y ago

nayuki 8y ago

I think the section "Using normally distributed random numbers" deserves an additional explanation. The code may look simple, but like most of the other methods mentioned on the page, generating normally (Gaussian) distributed random numbers still requires the use of rejection sampling and sine/cosine functions. See https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

strainer 8y ago
I implemented Marsaglia's Box Muller method in a library (^) and it seems quite quick. It doesn't involve trig functions, just a sqrt, a log, and an average of about two rands per normal outputted. It seems to work out around the same speed as a trig function - in javascript. For some purposes just averaging a few rands could be used to get a good enough approximation of a normal. I wonder if a short Taylor series could be mined that could quite accurately map equal to gaussian distribution. [^]http://strainer.github.io/Fdrandom.js/ http://strainer.github.io/Fdrandom.js/

nalourie 8y ago

For people who like the simplicity of the rejection algorithm, it might be worth checking out Ziggurat sampling which is based on a similar idea, but shapes the sampling region to reduce rejections: https://en.m.wikipedia.org/wiki/Ziggurat_algorithm https://en.m.wikipedia.org/wiki/Ziggurat_algorithm

charmides 8y ago
To add to that comment, importance sampling is also an attempt to reduce the number of rejections and could be interesting here. https://en.m.wikipedia.org/wiki/Importance_sampling https://en.m.wikipedia.org/wiki/Importance_sampling
- throwaway313213 8y ago
  Importance sampling works when you're interested in some linear function[al] of the distribution; not for actual sampling.

cjhanks 8y ago

Maybe a mathematician can correct me. But could you not also model 3D vertices of uniform density into a 1D space and then uniformly sample the 1D space. Of course this may have issues with highly dense objects or very large ones.

lixtra 8y ago
I don’t get what you mean by 1D spaces but you could triangulate into tetrahedra and sample those[1] (taking their volume into account). You will write a lot of code and still have the triangulation approximation error. For high dimensional spheres you can possibly subdivide the space in hyper cubes and do the in/out approach inside each cube that contains a part of the sphere to get a better hit miss ratio. [1] http://vcg.isti.cnr.it/jgt/tetra.htm http://vcg.isti.cnr.it/jgt/tetra.htm
- cjhanks 8y ago
  Yes, I was thinking you could triangulate the sphere at discreet shells. I have done 2D symmetric polygons where you needed to uniquely index into a 1D array of vertex values. Ie, based on the magnitude you can estimate which shell. By knowing the shell, you know the number of vertices for the shell, by knowing both you can map it to a specific vertex, as you can define a priori the absolute traversal order. I have no idea if the 2D case extends to higher dimensions.

herbertkane 8y ago

Section 3 of this paper explains this topic pretty well in N-dimensions: http://jmlr.org/papers/volume17/blaser16a/blaser16a.pdf http://jmlr.org/papers/volume17/blaser16a/blaser16a.pdf To get a random point within the sphere (rather than on the surface of the sphere, as described in the paper) you also need to pick a random distance from the origin to randomly scale the point on the surface.

throwaway313213 8y ago

This is similar to the method used for generating random numbers using Box-Muller transform. https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

arnarbi 8y ago

The main issue with the second method, which isn't discussed at all, is that it's more likely to pick a direction towards the corners of the unit cube. You can even see that bias in the point cloud. The problem with choosing d uniformly is also better explained by pointing out that shells of larger radius have larger surface, thus they'll be underrepresented.

sgentle 8y ago

I thought this was kinda fun, and I was trying to figure out why you couldn't just take a random x and use it to find a y that preserves x²+y²=r² (ie by setting y to a random number between 0 and ±√(r²-x²)). In the end I figured out that you need to generate x the same way, but if you have, say: n=rand(r), x=rand(±√(r²-n²)), y=rand(±√(r²-x²)) it seems to work. Here's a version I extended to three dimensions, designed to be run in the JS console of the original article (it'll show up underneath the first sphere): new SphereSimulation(document.getElementById("spheres1"), () => { const rand = Math.random const root = (x) => rand() < 0.5 ? Math.sqrt(x) : -Math.sqrt(x) const invdist = (a=0, b=0) => root(RADIUS*RADIUS - a*a - b*b) const x1 = rand() * invdist() const x2 = rand() * invdist(x1) const y2 = rand() * invdist(x2) const x3 = rand() * invdist(x2, y2) const y3 = rand() * invdist(x3, y2) const z3 = rand() * invdist(x3, y3) return new THREE.Vector3(x3, y3, z3) }) It looks okay to me - am I missing something? Is this a reasonable approach?

goldenkey 8y ago
I posted a comment where I described this method as my preferred solution. You did not need to patch it up though..its fine to just generate the 3 numbers and then choose a random assignment from them to x y and z. This fixes the bias that would be present if you always generated x first.
jpeanuts 8y ago
Your missing the requirement that the random sample be uniform in the sphere. That is, if the sphere is cut into any number of pieces of equal volume, the sample is equally likely to appear in any piece. Your approach samples uniformly from the radius, so you're going to end up with far too many samples close to the center.
- sgentle 8y ago
  Interesting... thanks for the pointer. I did some testing and the points definitely do seem to cluster along the axis... just increasing the number of iterations of √(r²-n²) seemed to improve it though. At 5, I can't tell the difference anymore: https://codepen.io/anon/pen/xJJOEq https://codepen.io/anon/pen/xJJOEq

goldenkey 8y ago

The method I immediately think of is just based on the relation. Sqrt(x^2 + y^2 +z^2) <= 1 reduces to x^2 + y^2 +z^2 <= 1 Generate the random value for x^2, any random value <= 1 Now we have y^2 +z^2 <= 1-x^2 Generate a random value for y^2, any random value <= 1-x^2 Finally generate a random value for z^2, any random value <= 1-x^2-y^2 Take the square roots of each number we generated to get the final coordinates. But wait! Theres another step..the most important one! There are 6 orderings of (x,y,z) We must choose a random one because otherwise x will be biased towards higher values compared to y and z. So I kind of lied. All we needed was 3 random values via the subtraction methods above. Then we need to randomly choose their ordering as (#,#,#) This seems like the easiest most trivial solution that would be uniformly distributed. Let me know if I am wrong.

kgwgk 8y ago
You are wrong.
machiavelli1024 8y ago
>Let me know if I am wrong. You're wrong on multiple levels. 1. Generating random values for x^2, y^2, z^2 and taking their square root will only give you values in the x,y,z > 0 octant. (But let's say you "fix" this by randomly multiplying them with -1.) 2. Taking the square root of a uniformly distributed random variable is no longer uniformly distributed. 3. Randomly reordering the coordinates won't fix your bias. Here's a demonstration in 2D: https://jsfiddle.net/tz85wnxy/59/ https://jsfiddle.net/tz85wnxy/59/ Uncomment line 17, 18, 20 to see how it's still not uniform even if you randomly multiply the coordinates by -1 and reorder them.
- goldenkey 8y ago
  Works fine as long as you generate x,y,z randomly and not their squares. Same process but minor tweak. https://jsfiddle.net/w8zLvsy2/ https://jsfiddle.net/w8zLvsy2/
  - machiavelli1024 8y ago

enugu 8y ago

For the second example, sampling the radius uniformly leads to an error with a higher point density in the centre of sphere. The radius, r, has a weight proportional to r^2 or the surface area. So, one can take a uniform random sample, s between 0 and R^2, and take the square root of s. That will give higher probability density around larger radius. Then, the second picture should look like the first one.

nwellnhof 8y ago
Seems like you stopped reading after the second example.
- enugu 8y ago
  Yes, didnt look at the details of the next methods. He adjusts for this in the next examples by taking the cube root of the uniform sample(not square root as in my post), which is the correct adjusting factor, as a small shell will have weight in proportion to difference between cubic powers, which is proportional to r^2.
deleted 8y ago
[deleted]

anewhnaccount2 8y ago

Similar problem which ends up not being quite as simple in the general case as you might think: sampling from a unit simplex: http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf

jjgreen 8y ago

http://blog.geomblog.org/2013/01/a-sampling-gem-sampling-from-ellp-balls.html http://blog.geomblog.org/2013/01/a-sampling-gem-sampling-fro...

expuexto 8y ago

I don't agree with the way he normalizes the radius in the last 2 methods. I think it should be a square root, not a cubic one. This is because the volume diferential r^2•sin(phi)•d_phi•d_theta•dr is proportional to the square of r. This mistake is also made in a source he cites: https://math.stackexchange.com/questions/87230/picking-random-points-in-the-volume-of-sphere-with-uniform-probability/2872878 https://math.stackexchange.com/questions/87230/picking-rando... and I tried to point it out there as well.

improbable22 8y ago
That's the correct volume element, but not what you want. The volume inside radius r grows as r^3. It's the inverse of this _cumulative_ density function which maps uniform-in-r numbers to uniform-in-volume ones.
antidesitter 8y ago
Did you empirically check your hypothesis? What's the integral of r^2 dr? Think about whether your claim makes sense when n = 1.
expuexto 8y ago
You both are right. Thank you for the comments.

plg 8y ago

That (last method) is a lot of sines and cosines, cube root and an acos... I bet for 3D the rejection method is not that much slower, in wall clock time. OK so I coded up the first one and the last one in C and compiled on my mac using whatever gcc links to. I found 1e7 points in the sphere: first routine: 0.521s last routine: 0.850s [edit]: C code here: https://pastebin.com/zH0BehSM https://pastebin.com/zH0BehSM

froogie 8y ago
The first routine is about 3 rands per point minus rate of discarded, about half? So 6 rands per point. Now, since the second one has three rands in it, and plenty of trig, it really becomes a question if you can make the trigonometry faster than three rands. But as you saw, straightforward attempt will pretty much certainly end up having worse performance. But what if one used optimizations like lookup tables and trigonometric approximations instead? Could they end up costing less cycles than three rands? Maybe not LUTs, maybe other trickery someone knows of? And what about the cost of those three rands? Could one use a fast xorshift with good enough results here? Does one even need thee random vectors? Maybe two + modulo trickery will do?
- acqq 8y ago
  Typically the generation of the uniform random numbers can be made quite cheaply, especially if the resolution is not critical (e.g. if uniform floats are enough and not uniform doubles etc) or if the quality of the randomness is not critical, and in these cases discarding half of these will win versus other costly operations.
improbable22 8y ago
Note that his last way does cos(acos(cosϕ))) after generating cosϕ uniformly. Which is a waste, you can just set sinϕ = sqrt(1-cosϕ^2). After fixing that, this way is faster for me. Some trig, but no branches. Edit: https://pastebin.com/UExZ5Siw https://pastebin.com/UExZ5Siw
- shpx 8y ago

hatsunearu 8y ago

So what ends up being the general purpose method for uniformly distributed random numbers in some arbitrary shape? (other than the rejection method of course) If arbitrary shapes always requires the rejection method, what are some "classes" that have easier algorithms (the sphere is obviously included here)

smaddox 8y ago
Not my area of expertise, but based on what I've read, you general try to get as close as possible through pure mappings, and then use rejection sampling to throw out the samples that don't fall under the desired distribution. This can work even if the distribution is dynamically changing.

wruza 8y ago

Interestingly, it seems to be all about bending some distribution and then unbending it again. Orthogonal dimensions are independent, so choosing three random numbers generates uniform in-cube distribution. But once you bend it via polar coordinates or vectors, you “wrap” many points of polar/vector abstract space to near the cartesian origin. And then you have to bend back, to make it uniform again. That raises interesting question on whether it relates to topology problems. Since generating N randoms is a fastest primitive operation, is there a quick and uniform translation from a cube to the sphere, such that all points are still uniform? Can we apply space-filling curve magic here (like Hilbert’s)? I feel like there is much more geometry power under that, but sadly I’m just a guy with no math bg.

3pt14159 8y ago
It's an interesting idea, but I think it would prove to be equivalent. As another commenter pointed out, the efficiency of this operation really only matters in higher dimensions. I bet a generative model would work, though it's tricky enough for me not to be able to think it through on the spot. Something like: 1. Pick a subset of points thus far generated at random. 2. From this subset, pick the point with the furthest distance from the origin. 3. Form a new unit sphere where r is 1-d of the picked point (ie, make a new sphere that is a subset of the parent sphere with a single touching point of the parent sphere). 4. Generate a new point in this subsphere using the normalized vector method. Not exactly the above, but something like it. The downside of the above approach is that I think it will create a distribution of points that may be unnaturally uniform. It would be interesting to define a set of equations that could benchmark different approaches for how well they compare to the in-cube strategy.

chias 8y ago

Interestingly, the method in which you determine your "random point" can greatly affect the end result, so swapping out the methods (e.g. "generate three coordinates and discard" method for the "generate a direction and a length" method) may produce wildly different results. An interesting illustration of this is Bertrand's Paradox, which asks the question: "if you draw a random chord inside a circle, what is the probability that its length is greater than sqrt(3)?". The correct answer is 1/2. The correct answer is also 1/3. Oh, and the correct answer is also 1/4. It depends on three equally valid (and uniformly distributed) ways of defining your random chord. See: http://web.mit.edu/tee/www/bertrand/ http://web.mit.edu/tee/www/bertrand/

improbable22 8y ago
Paradox seems like a strong word for this. "Draw a random chord" is almost a pun, draw from what distribution? Until you tell me, there really is no one answer.
- chias 8y ago
  > "Draw a random chord" is almost a pun No more so than 'pick a random point in a circle'. Note in all cases we're talking the uniform distribution. The question of using Cartesian coordinates, parametric coordinates, etc. is what causes the differences.
  - improbable22 8y ago
    I'm complaining that there isn't one uniform distribution. For example uniform in x on [1, 2] is different to uniform in 1/x on the same interval, and which one you draw from may matter. You have to be told. In fact every (one-variable) distribution is uniform in something! "pick a random point inside a circle" comes with an implicit "uniformly according to the ordinary euclidean metric on the paper". I guess it's a judgement call that this is obvious enough, but "a random chord" is clearly not... otherwise we wouldn't be here!
deleted 8y ago
[deleted]

Myrmornis 8y ago

I find it unclear whether the author is aiming for Uniform, or just “looks pretty uniform”.