|
|
Re: joint normal distribution integral over the unit sphere
Posted:
Sep 8, 2004 11:04 PM
|
|
Ray Koopman wrote:
>
> Here's some Mathematica code that seems to give 7-place accuracy.
Earlier today I posted a message very similar to the following,
but it seems to have gotten lost.
Here's some code that is simpler but more accurate than "vol",
and increasing PrecisionGoal will make it more accurate yet.
Experiments (not reported here) indicate that integrating w.r.t. r
seems to give more accurate results than integrating w.r.t. r^2,
even though the latter is simpler.
In[1]:= newvol[sx_,sy_,sz_] := With[{a = -.25(sx^-2 + sy^-2),
b = .25(sx^-2 - sy^-2)}, NIntegrate[r Exp[a*r^2] *
Erf[Sqrt[.5(1.-r^2)]/sz] BesselI[0,b*r^2],{r,0,1}]/(sx*sy)]
In[2]:= InputForm[ ColumnForm[ Map[ newvol[ Sequence@@#]&,
Permutations@{.3,.4,.5}]]]
Out[2]//InputForm=
0.8835309249169815
0.8835309249174408
0.8835309249169815
0.8835309249181131
0.8835309249174408
0.8835309249181131
In[3]:= With[{s = 1.3}, InputForm @ ColumnForm @
{1-GammaRegularized[ 3/2,.5/s^2],newvol[s,s,s]}]
Out[3]//InputForm=
0.10167383475354619
0.10167383475389877
|
|
