MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • 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.