RE: Bivariate InverseErf funct
- To: mathgroup@smc.vnet.net
- Subject: [mg11465] RE: [mg11367] Bivariate InverseErf funct
- From: Ersek_Ted%PAX1A@mr.nawcad.navy.mil
- Date: Thu, 12 Mar 1998 01:34:32 -0500
Albert Maydeu-Olivares wrote:
|
|I need to perform an inverse interpolation of the bivariate normal
|distribution function.
|In other words, I want to estimate a tetrachoric correlation, so I need
|a bivariate version of the InverseErf function. A glance at the
|standard references for this problem make me believe that programming
|this from scratch may be rather time consuming. I'm interested in
|knowing whether anyone has done it before, or whether I can use some
of |the built-in Mathematica functions to do the job. |
|In a nutshell, this is the problem
|
|
|Consider a bivariate normal distribution with mean {0,0} and
correlation |matrix {{1,r},{r,1}} that has been dichotomized at m1 and
m2. Suppose |m1,m2 & r are known,
|
|m1 = 1; m2 =-1; r=.5;
|
|then we can compute
|
|p1 = N[1 - CDF[NormalDistribution[0, 1], m1]]; p2 = N[1 -
|CDF[NormalDistribution[0, 1], m2]];
|p21=NIntegrate[PDF[MultinormalDistribution[{0,0},{{1,r},{r,1}}],{x1,x2}],{x1
|,m1,\[Infinity]},{x2,m2,\[Infinity]}]; |
|int = NIntegrate[1/(E^((m1^2 + m2^2 - 2*m1*m2*x)/(2*(1 - x^2)))*Sqrt[1
|-x^2]), {x, 0, r}]/ (2*Pi);
|
|What I am interested in is the inverse problem, given numerical
|estimates of p1,p2 and p21 to estimate m1, m2 and r. |
|m1 and m2 can be estimated using the built-in InverseErf function. To
|estimate r from p21 for given m1 and m2 I would need something like a
|bivariate InverseErf function.
|
|
|Without it, this is something I've tried: An anternative expression for
|p21 is
|
|int = NIntegrate[1/(E^((m1^2 + m2^2 - 2*m1*m2*x)/(2*(1 - x^2)))*Sqrt[1
|-x^2]), {x, 0, r}]/ (2*Pi);
|p21= int + p1*p2;
|
Albert,
I think it will work very well if you use FindRoot. I start with what
you did.
In[1]:=
<<Statistics`MultinormalDistribution`
In[2]:=
m1 = 1;
m2 =-1;
In[3]:=
p1 =N[1 - CDF[NormalDistribution[0, 1], m1]]
Out[3]=
0.158655
In[4]:=
p2 = N[1 -CDF[NormalDistribution[0, 1], m2]]
Out[4]=
0.841345
In[5]:=
distr=1/(E^((m1^2 + m2^2 - 2*m1*m2*x)/(2*(1 - x^2)))*Sqrt[1
-x^2])/(2 Pi);
In[6]:=
int[r_]:=NIntegrate[Evaluate[distr], {x, 0,r}]
__________________________________
Then I use FindRoot to define the inverse function.
In[7]:=
InverseMultiNormal[p21_]:=
FindRoot[int[r]+p1*p2==p21,
{r,0.000001,0.9999999},Compiled->False]
Lets see how it works.
Now I find r for (p12=0.15) and the values given above for (p1), (p2),
(m1), (m2).
In[8]:=
soln=InverseMultiNormal[0.15]
Out[8]=
{r\[Rule]0.343179}
Now I make sure it works.
In[9]:=
int[r/.soln]+p1*p2
Out[9]=
0.15
Hope it helps.
Ted Ersek