version 1.0.0.0 (542 Bytes) by
Roger Stafford

Randomly and uniformly distributes points throughout a hypersphere.

8 Downloads

Updated 28 Dec 2005

No License

This creates a set of random points defined by Cartesian coordinates and uniformly distributed over the interior of an n-dimensional hypersphere of radius r with center at the origin.

The 'randn' function is first used to create independent multivariate normally distributed sets of n random variables, each representing points in n-dimensional space. The incomplete gamma function, 'gammainc', is then used to map these points radially to the interior of an n-dimensional hypersphere of finite radius r in such a way as to be spatially uniformly distributed.

Roger Stafford (2021). Random Points in an n-Dimensional Hypersphere (https://www.mathworks.com/matlabcentral/fileexchange/9443-random-points-in-an-n-dimensional-hypersphere), MATLAB Central File Exchange. Retrieved .

Created with
R10

Compatible with any release

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Sarafa Iyaniwuracheulhee kwonOndrej PanekThanks a lot!

shdotcom shdotcomWhat does S2 mean?

How to make the points distribute around the a center point ?

Elizabeth HouNever mind. My bad I had a plotting typo.

Elizabeth HouThe data seems to lie on a disk in n-dimensional space. Is there a way to make it an actual sphere?

Fjedor GaedeSreejita Ghoshsheren makhoulThank you:)

Tian LIUnice function

ErfanThank you very much for sharing information!

P.-O.Job% Super elegant!

% For the specific 3D case, howeverm, there is a 4x faster solution:

radius = r*rand(n,1).^(1/3);

phi = 2*pi*rand(n,1);

costheta = 1 - 2*rand(n,1);

radsintheta = radius.*sin(acos(costheta));

X = [...

radsintheta.*cos(phi), ...

radsintheta.*sin(phi), ...

radius.*costheta];

JobKuan-SungHi,

In AMB's code, we should generate X again before we perform Laurent's approach. Then the results from 2 approaches are similar.

Thanks,

Kuan-Sung

AMBHi - This post is in reference to Laurent's comments.

Using this sample code, we see that Laurent's code substitution is faster, but the result is different from the original. Laurent's modification appears to produce data that spans a different radius than the original and is comprised of data which when viewed on the screen in a 2D projection appears more centrally clustered. The original code produces more defined boundaries.

Both ideas may have their own virtues depending on application.

m = 30000;n = 2; r = 2;

X = randn(m,n);

s2 = sum(X.^2,2);

tic

X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);

toc

subplot(1,2,1);plot(X(:,1),X(:,2),'r.','markersize',1);

axis equal;zoom off; zoom on;drawnow;shg;

ax = axis;

tic

X = X.*repmat(r*(rand(m,1).^(1/n))./sqrt(s2),1,n);

toc

subplot(1,2,2);plot(X(:,1),X(:,2),'r.','markersize',1);

axis equal;zoom off; zoom on;drawnow;shg;

axis(ax);

pause

m = 30000;n = 3; r = 2;

X = randn(m,n);

s2 = sum(X.^2,2);

tic

X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);

toc

subplot(1,2,1);plot3(X(:,1),X(:,2),X(:,3),'r.','markersize',1);

axis equal;rotate3d off; rotate3d on;drawnow;shg;

ax = axis;

tic

X = X.*repmat(r*(rand(m,1).^(1/n))./sqrt(s2),1,n);

toc

subplot(1,2,2);plot3(X(:,1),X(:,2),X(:,3),'r.','markersize',1);

axis equal;rotate3d off; rotate3d on;drawnow;shg;

axis(ax);

pause

m = 30000;n = 20; r = 2;

X = randn(m,n);

s2 = sum(X.^2,2);

tic

X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);

toc

subplot(1,2,1);plot3(X(:,1),X(:,2),X(:,3),'r.','markersize',1);

axis equal;rotate3d off; rotate3d on;drawnow;shg;

ax = axis;

tic

X = X.*repmat(r*(rand(m,1).^(1/n))./sqrt(s2),1,n);

toc

subplot(1,2,2);plot3(X(:,1),X(:,2),X(:,3),'r.','markersize',1);

axis equal;rotate3d off; rotate3d on;drawnow;shg;

axis(ax);

Vincent DubourgHi again,

I googled a little more and finally found a short explanation in the book by Rubinstein (2008, p. 69). I leave the link here in case it deserves someone else...

http://books.google.com/books?id=1-ffZVmazvwC&pg=PA69&lpg=PA69&dq=random+hypersphere+rubinstein&source=bl&ots=-q859Ytx3q&sig=q-PtG2P62_kTCgzi-J7QuqvPOZA&hl=fr&ei=zZ25TcfUJ4SChQeIzY3_Dg&sa=X&oi=book_result&ct=result&resnum=1&ved=0CBwQ6AEwAA

Cheers,

Vincent

Vincent DubourgHi!

Thank you for sharing this neat function! I was wondering if you could give me any reference explaining the theoretical concepts underlying this script?

Thanks,

Vincent

LaurentThe convenient use of the incomplete gamma function luckily provided by Matlab can be avoided, so as to make the function probably faster, more elegant and reliable. No non-elementary function is necessary here.

In the last line, simply replace

gammainc(s2/2,n/2)

by

rand(m,1)

This amount to sampling the radius randomly, independently of the point "X/sqrt(s2)" that is uniform on the surface of the sphere. Indeed, one can check that if u is uniform on [0,1] then r*u^(1/n) has same distribution as the norm of a vector chosen uniformly inside the ball B(0,r) in n-dimensional space.

I don't have matlab right here to check for speed compared to using gammainc, but in any case this method does not rely on some blackbox approximate integration, which makes it, to my taste, more elegant. It is also more portable into other programming languages.

chen ?Thanks a lot!

I have another question.If I wanna to get a uniform distribution on the surface of the hypersphere (radius equals 1), whether is it correct to normalize the generated vectors to be length 1. Thank you very much!

JiaqiFlorianThe elegancy is breathtaking! I wish you would post more functions!

molc morlocMathMooseAree WitoelarThanks a lot for the elegant code! I had trouble generating this distribution. Before this, my solution was to sample and reject; but for high dimensions it was too crude and consuming.

Steve LordVery nice piece of code. I have just two small suggestions for enhancements to the M-file help:

1) Include the calling syntax early in the help. I will admit I sometimes forget the calling syntax for functions I don't use often; my first resort is to type "help <functionname>", either to refresh my memory or to copy and paste the syntax into my code. In a case like this, I can see users perhaps forgetting the order of the parameters.

2) Include an example. Having a example that users can copy and paste included in the help gives people who may not know how to use your function a "known good" state -- from there, they can tweak. For a straightforward function like this, it may not seem necessary, but I think it's a good habit. It could be as simple as the two lines of code to generate the sample image.

Lee Joneselegant procedure, surprisingly

Peter NaveThe compactness of the code is stunning.

John D'ErricoI'll admit it took me a minute or so to verify this code. Its

mathematically elegant. Its efficient. It even works in the

special case of n = 1.

Peter (PB) BodinSimple, yet very clever solution.