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