       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, 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=
{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];

FixedPoint[step,.2,SameTest\[Rule](Abs[#1-#2]<10.^(-6)&)]

But I get this error which I find completely uninformative.

>From In:=
NIntegrate::"nlim": "\!\(x\$60\) = \!\(r\) is not a valid limit of \
integration."

However, everything seems to be in order, see

r=.3;
f
step[r]
Clear[r]
Out=
0.178779
Out=
\!\(7.331185010880925`*^-6\)
Out=
0.3

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:=
step[r_] :=r-(cdf2[t1,t2,r]-obs[])/pdf2[t1,t2,r];

FixedPointList[step,.2,SameTest->(Abs[#1-#2]<10.^(-6)&)]//Timing

Out=
{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

```

• Prev by Date: RE: Why Does AbsoluteOptions Not Tell about all Automatic?
• Next by Date: Re: Mathematica 3.0: reliability close to LogZero?
• Previous by thread: Re: Speeding up Replacement Rules
• Next by thread: HairyPooter IV.75 cont.