2

I am trying to generate points uniformly distributed on a region of a single n-sphere surface (all angles in hyperspheric coordinates are between 0 and Pi/2).

I decided to use von Neumann's method to generate such a distribution from an uniform distribution. Thus, I need to generate values with probability density of

norm[dimension_?NumberQ]=:=1  /  (  dimension*Pi^(dimension/2)/Gamma[1+dimension/2]/2^dimension  )  ;
newDens[arr__]:=Module [ {var=arr}, norm[var]*Product[  (  Sin [ var[[k]]   ] )^(k-1)  , {k,2,  Length[var] }  ]  ];

(see https://en.wikipedia.org/wiki/N-sphere)

So, I generate "angles" uniformely in (n-1)cuboid [0, Pi/2] and values between 0 and norm[n] and then compare them with newDens[angles].

Here is my code:

generateNeumann[count_?IntegerQ, dimension_?NumberQ]:=(
  norm=1  /  (  dimension*Pi^(dimension/2)/Gamma[1+dimension/2]/2^dimension  )  ;
  normN=N [norm];
  newDens[arr__]:=norm*Module [ {var=arr}, Product[  (  Sin [ var[[k]]   ] )^(k-1)  , {k,2,  Length[var] }  ]  ];

  pointsDistr={};

  For[i=1, i<=count, i++, {
    angles=RandomReal[{0,Pi/2},(dimension-1)  ],
    value=RandomReal[normN],
    surface=newDens[angles],
    If[value<surface, AppendTo[pointsDistr, angles]]

   }
 ]

);

Obviously, Length[pointsDistr]/count should be about (Pi/2)^(1-dimension)*norm^(-1).

Let's try it for dimension=3.

Timing [  generateNeumann[10^3,3]  ]
Length[  pointsDistr  ]
normN

(*  {0.046875,Null}  *)
(*655*)
(*0.63662*)

It seems to be good:

 decartePointsDistr=Table[{ Sin[pointsDistr[[i,2]]]*Sin[pointsDistr[[i,1]]],   Sin[pointsDistr[[i,2]]]*Cos[pointsDistr[[i,1]]],  Cos[ pointsDistr[[i, 2]]   ]  },{i,1,Length[pointsDistr]}];
 Show[  SphericalPlot3D[1, theta, phi],ListPointPlot3D[decartePointsDistr]  ]

Trying to increase the number of points:

Timing [  generateNeumann[10^4,3]  ]
Length[  pointsDistr  ]

(*{0.546875,Null}*)
(*6443*)

OK too.

But:

Timing [  generateNeumann[10^5,3]  ]
Length[  pointsDistr  ]

(*{20.484375,Null}*)
(*21808*)

It's very slow and not correct.

What's wrong with large number of points? In R, such a procedure works very fast even when count=10^6.

Maybe I'm doing something wrong?

Thanks for your help!

Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
Sasha
  • 195
  • 6
  • 4
    First likely problem: AppendTo within a For loop. Replace both with Table or Array and report your findings. – Mr.Wizard Feb 12 '15 at 12:33
  • What Mr.wizard says. Anyway, it looks like the commas in the body of your For loop should be semicolumns. – Sjoerd C. de Vries Feb 12 '15 at 13:01
  • @Sjoerd the body is contained with { } so I think that part is OK. – Mr.Wizard Feb 12 '15 at 13:15
  • @Mr.Wizard what should I do with conditions (If[]) when using Table? – Sasha Feb 12 '15 at 13:19
  • 2
    @Sasha My method is here: (3705) If most of the elements will be dropped it is more efficient to use Sow and Reap along with Do. – Mr.Wizard Feb 12 '15 at 13:20
  • 1
    @Mr.Wizard Yes, it works. I've just used pointsDistr=If[newDens[#]>RandomReal[normN], #, ##&[]]&/@RandomReal[{0,Pi/2}, {count, dimension-1}] instead of the loop. What's the trouble with AppendTo within the loop? – Sasha Feb 12 '15 at 14:56
  • 1
    In Mathematica lists are internally implemented as C arrays (as I understand it) and changing the length of the list results in reallocation of the array. The longer the array the slower the append operation. – Mr.Wizard Feb 12 '15 at 15:07
  • 1
    @Mr.Wizard Yes, but I don't like the style. It kind of builts up a return vector (which of course isn't retuned as For doesn't do that). The canonical way of making a For body is (IMHO) using no brackets and separating statements using semicolons. – Sjoerd C. de Vries Feb 12 '15 at 15:42

0 Answers0