# Von Neumann's method to generate points uniformly distributed on a region of a n-sphere surface

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] }  ]  ];


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?

• 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
• @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
• @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