Problem using NIntegrate within FixedPoint
- To: mathgroup at smc.vnet.net
- Subject: [mg24519] Problem using NIntegrate within FixedPoint
- From: Albert Maydeu-Olivares <amaydeu at tinet.fut.es>
- Date: Mon, 24 Jul 2000 03:04:02 -0400 (EDT)
- Organization: Universitat de Barcelona
- Sender: owner-wri-mathgroup at wolfram.com
Hello,
I am not able to use NIntegrate within FixedPoint in one piece of code.
I would very much appreciate if anyone has a clue of what the problem is
and how to solve it.
These are univariate and bivariate normal CDF functions
pdf1[h_] := N[1/(E^(h^2/2)*Sqrt[2*Pi]), 16];
pdf2[h_, k_, x_] := N[E^((-(k*(k/(1 - x^2) - (h*x)/(1 - x^2))) -
h*(h/(1 - x^2) - (k*x)/(1 - x^2)))/2)/(2*Pi*Sqrt[1
-x^2]), 16];
cdf1[h_] := Chop[N[(1 + Erf[h/Sqrt[2]])/2, 16]];
cdf2[h_, k_, r_] := Chop[Module[{x}, NIntegrate[
E^((-(k*(k/(1 - x^2) - (h*x)/(1 - x^2))) -
h*(h/(1 - x^2) - (k*x)/(1 - x^2)))/2)/(2*Pi*
Sqrt[1 - x^2]), {x, 0, r}, PrecisionGoal -> 10] +
cdf1[h] cdf1[k]]];
der[h_, k_, x_] :=
N[-(E^((h^2 + k^2 - 2*h*k*x)/(-2 + 2*x^2))*(h^2*x - h*k*(1 + x^2) +
x*(-1 + k^2 + x^2)))/(2*Pi*(1 - x^2)^(5/2)), 16];
In case you are curious, the objective is to estimate the underlying
correlation in a standard BVN distribution that has been dichotomized at
t1 & t2 using MLE
via Newton's method.
These are the prob in a 2x2 table according to the model
t1=1.;t2=-1.;r=.3;
obs={cdf2[t1,t2,r],
cdf1[t1]-cdf2[t1,t2,r],cdf1[t2]-cdf2[t1,t2,r],1-cdf1[t1]-cdf1[t2]+cdf2[t1,t2,r]}
Clear[r];
Out[61]=
{0.148338, 0.693007, 0.010317, 0.148338}
Given this "data" and assuming I know t1 and t2, this is the MLE
solution
exp:={cdf2[t1,t2,r], cdf1[t1]-cdf2[t1,t2,r],
cdf1[t2]-cdf2[t1,t2,r],1-cdf1[t1]-cdf1[t2]+cdf2[t1,t2,r]};
first:={pdf2[t1,t2,r], -pdf2[t1,t2,r], -pdf2[t1,t2,r],pdf2[t1,t2,r]};
second:={der[t1,t2,r], -der[t1,t2,r], -der[t1,t2,r],der[t1,t2,r]};
f:=Apply[Plus,obs * first / exp];
grad:=Apply[Plus,Table[(obs[[i]]*first[[i]]^2+obs[[i]]*exp[[i]]*second[[i]])/exp[[i]]^2,{i,4}]];
step[r_] :=r- f/grad;
FixedPoint[step,.2,SameTest\[Rule](Abs[#1-#2]<10.^(-6)&)]
But I get this error which I find completely uninformative.
>From In[60]:=
NIntegrate::"nlim": "\!\(x$60\) = \!\(r\) is not a valid limit of \
integration."
However, everything seems to be in order, see
r=.3;
grad
f
step[r]
Clear[r]
Out[52]=
0.178779
Out[53]=
\!\(7.331185010880925`*^-6\)
Out[54]=
0.3
Any comments/suggestions?
As a matter of fact, this problem has a much simpler solution using any
of the 4 cells, which works well giving me the true r in just 4
iterations. For example using the first cell, the solution is
In[58]:=
step[r_] :=r-(cdf2[t1,t2,r]-obs[[1]])/pdf2[t1,t2,r];
FixedPointList[step,.2,SameTest->(Abs[#1-#2]<10.^(-6)&)]//Timing
Out[59]=
{0. Second, {0.2, 0.293054, 0.299959, 0.3, 0.3}}
I'm really puzzled about why the first piece of code does not work and
the second one does.
Albert