Bivariate InverseErf function
- To: mathgroup@smc.vnet.net
- Subject: [mg11367] Bivariate InverseErf function
- From: "Albert Maydeu-Olivares" <amaydeu@tinet.fut.es>
- Date: Sat, 7 Mar 1998 02:06:39 -0500
- Organization: University of Barcelona
Hi everyone,
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;
then, in principle this should work...
Clear[m1,m2,r];
m1 = Sqrt[2]*InverseErf[0, 1 - 2*p1] m2 = Sqrt[2]*InverseErf[0, 1 -
2*p2]
j =
(2\[Pi])^-1*(1-r^2)^(-1/2)*Exp[(-1/2)*((m1^2-2*m1*m2*r+m2^2)(1-r^2)^-1)]
r = FindMinimum[p21 - p1*p2 - NIntegrate[1/(E^((m1^2 + m2^2
-2*m1*m2*x)/(2*(1 - x^2)))*Sqrt[1 - x^2]), {x, 0, r}]/
(2*Pi),{r,.5},Gradient->j,Method->Newton]
which fails miserably.
(Someone suggested to use the secant method in FindMinimum and did not
work either)
Albert
--
Albert Maydeu-Olivares Tel. +34 3 4021079 ext. 3099 Faculty of
Psychology Fax. +34 3 4021362 University of Barcelona E-Mail:
amaydeu@tinet.fut.es Passeig de la Vall d'Hebron, 171.
08035 - Barcelona (Spain)