home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / STATISTI.PAK / INVERSES.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  11.8 KB  |  407 lines
  1.  
  2. (*:Version: Mathematica 2.0 *)
  3.  
  4. (*:Name: Statistics`InverseStatisticalFunctions` *)
  5.  
  6. (*:Title: Evaluating the Inverse of Erf[ ], GammaRegularized[ ], and
  7.         BetaRegularized[ ] *)
  8.  
  9. (*:Author:
  10.   Jerry B. Keiper (Wolfram Research, Inc.), March, 1989.
  11.   Current version April, 1991.
  12. *)
  13.  
  14. (*:Summary:
  15.   Defines the inverses of the following statistical functions:
  16.   Erf[ ], Erfc[ ], GammaRegularized[ ], and BetaRegularized[ ].
  17.   The algorithms are similar to Newton's method except that, rather
  18.   than using a linear approximation, higher-order approximations are
  19.   used leading to cubic, and in the case of InverseErf[ ], quartic
  20.   convergence.
  21. *)
  22.  
  23. (*:Keywords: inverse function, Newton's method, InverseErf[ ],
  24.   InverseErfc[ ], InverseGammaRegularized[ ], InverseBetaRegularized[ ]
  25. *)
  26.  
  27. (*:Requirements: no special system requirements. *)
  28.  
  29. (*:Warnings: none. *)
  30.  
  31. (*:Sources:
  32.   M. Abramowitz & I. Stegun, Handbook of Mathematical Functions,
  33.     Dover, New York, 1972, chapter 26
  34. *)
  35.  
  36. BeginPackage["Statistics`InverseStatisticalFunctions`"]
  37.  
  38. InverseErf::usage = 
  39. "InverseErf[p] gives the value of x such that p == Erf[x].  InverseErf[x0,p]
  40. gives the value of x such that p == Erf[x0,x].  InverseErf[x0,-p] gives the
  41. value of x such that p == Erf[x,x0].  x0 and p must be real, and p must be
  42. between Erf[x0, -Infinity] and Erf[x0, Infinity]"
  43.  
  44. InverseErfc::usage =
  45. "InverseErfc[p] gives the value of x such that p == Erfc[x].  p must satisfy
  46. 0 <= p <= 2."
  47.  
  48. InverseGammaRegularized::usage =
  49. "InverseGammaRegularized[a,x0,p] gives the value of x such that p ==
  50. GammaRegularized[a,x0,x].  InverseGammaRegularized[a,x0,-p] gives the value
  51. of x such that p == GammaRegularized[a,x,x0].  a, x0, and p must be real,
  52. a > 0, 0 <= x <= Infinity, and GammaRegularized[a,x,0] <= p <=
  53. GammaRegularized[a,x,Infinity]."
  54.  
  55. InverseBetaRegularized::usage =
  56. "InverseBetaRegularized[x0,p,a,b] gives the value of x such that p ==
  57. BetaRegularized[x0,x,a,b].  InverseBetaRegularized[x0,-p,a,b] gives the
  58. value of x such that p == BetaRegularized[x,x0,a,b].  a, b, x0, and p must
  59. be real, a > 0, b > 0, 0 <= x <= 1, and BetaRegularized[x,0,a,b] <= p <=
  60. BetaRegularized[x,1,a,b]."
  61.  
  62. Begin["Statistics`InverseStatisticalFunctions`Private`"]
  63.  
  64. InverseErf[p_] := InverseErf[0,p];
  65. InverseErfc[p_] := InverseErf[Infinity,-p];
  66. InverseErf[x_, 0] := x;
  67.  
  68. InverseErf[x_,p_?NumberQ] :=
  69.     Module[{e},
  70.     e /; (e = inverf[x,p]) =!= $Failed
  71.     ] /; (Im[p]==0 && ((NumberQ[x] && Im[x]==0) ||
  72.         SameQ[x,Infinity] || SameQ[x,-Infinity]))
  73.  
  74. InverseGammaRegularized[a_?NumberQ,x_,p_?NumberQ] :=
  75.     Module[{e},
  76.     e /; (e = invgam[a,x,p]) =!= $Failed
  77.     ] /; (Im[p]==0 && Im[a]==0 && a>0 && (SameQ[x,Infinity] ||
  78.             (NumberQ[x] && Im[x]==0 && x >= 0)))
  79.  
  80. InverseBetaRegularized[x_?NumberQ,p_?NumberQ,a_?NumberQ,b_?NumberQ] :=
  81.     Module[{e},
  82.     e /; (e = invbet[x,p,a,b]) =!= $Failed
  83.     ] /; (Im[x]==0 && x >= 0 && x <= 1 && Im[p]==0 &&
  84.             Im[a]==0 && a>0 && b>0 && Im[b]==0) 
  85.  
  86. mp[j_] := Max[j, $MachinePrecision]
  87. precis[x_] := If[MachineNumberQ[x], $MachinePrecision, Precision[x]]
  88.  
  89. invgam[a_,x_,pp_] :=
  90.     Module[{prec=Min[precis[a],precis[x],precis[pp]],p,x0,
  91.         gamx,gama,cflag,j,pa,temp,app,dy,x1,x2},
  92.         If[prec == $MachinePrecision, prec--];
  93.         If[SameQ[prec,Infinity], Return[$Failed]];
  94.         gama = N[Gamma[a],prec];
  95.         If[SameQ[Head[x],DirectedInfinity],
  96.             cflag = True;
  97.             gamx = 0;
  98.             p = pp+1,
  99.             (* else *)
  100.             If[cflag = (x>a),
  101.                 gamx = Gamma[a,N[x,prec],Infinity]/gama;
  102.                 p = pp+1 - gamx,
  103.                 (* else *)
  104.                 gamx = Gamma[a,0,N[x,prec]]/gama;
  105.                 p = pp + gamx
  106.                 ]
  107.             ];
  108.         (* p == GammaRegularized[a,0,x0] where x0 is still unknown
  109.         gamx is the value of the incomplete gamma function from
  110.         0 to x or x to Infinity depending on cflag. *)
  111.         If[p==0, Return[p]];
  112.         If[p<0 || p>1, Return[$Failed]];
  113.         If[a == 1,
  114.             p = If[cflag, gamx - pp, 1 - p];
  115.             Return[-N[Log[p],prec]]];
  116.         app = N[a];
  117.         pa = (p a gama)^(1/a);
  118.         temp = pa/(app+1);
  119.         x0 = pa(1+temp(1+temp((5+3app)/(2(app+2))+
  120.             temp((31+app(33+8app))/(3(2+app)(3+app))+temp(
  121.             (2888+app(5661+app(3971+app(1179+125app))))/
  122.                 (24(2+app)^2(3+app)(4+app))+
  123.             temp((9074+app(20997+app(18375+app(7575+app(
  124.                 1471+108app)))))/
  125.             (10(2+app)^2(3+app)(4+app)(5+app))))))));
  126.         If[x0 > 5, (* if x0 is large it may be very wrong *)
  127.             If[cflag,
  128.                 temp = inverf[Infinity,2(pp - gamx)],
  129.                 (* else *)
  130.                 temp = inverf[-Infinity,2 N[p]]];
  131.             x1 = app (1 - 1/(9app) + temp/Sqrt[4.5 app])^3;
  132.             (* if x1 is small it may be very wrong *)
  133.             If[x1 > 5, x0 = x1]];
  134.         If[x0 > app,
  135.             p = If[cflag, gamx - pp, cflag = True; 1 - p],
  136.             (* else *)
  137.             cflag = False];
  138.         If[p<0, Return[$Failed]];
  139.     (* now use Taylor series to quartically converge to x0. *)
  140.         p *= gama;
  141.         j = 0;
  142.         While[j < prec,
  143.             If[j > 0 && NumberQ[x0],
  144.                 j = Max[j,Min[4j,
  145.                     Max[1,If[x1 == 0,
  146.                         Max[4j,4],
  147.                         (* else *)
  148.                         Round[-1.7 Log[Abs[N[x1]]]]
  149.                     ]]]];
  150.                 If[j > prec, j = If[j > 2prec, prec, prec-1]];
  151.                 x0 = SetPrecision[x0,mp[j]],
  152.                 (* else *)
  153.                 j = 1
  154.                 ];
  155.             dy = p-If[cflag,Gamma[a,x0,Infinity],Gamma[a,0,x0]];
  156.             While[dy == 0 && j < prec,
  157.                 (* adding a fake zero would corrupt everything *)
  158.                 If[j < $MachinePrecision,
  159.                     j = 2 $MachinePrecision, j *= 2];
  160.                 If[j > prec, j = prec];
  161.                 x0 = SetPrecision[x0,mp[j]];
  162.                 dy = p-If[cflag,
  163.                     Gamma[a,x0,Infinity],
  164.                     Gamma[a,0,x0]]
  165.                 ];
  166.             If[dy == 0, Return[x0]];
  167.             dy = SetPrecision[If[cflag, -dy, dy],mp[j]];
  168.             x1 = Exp[x0] x0^-a dy;
  169.             x2 = x1((1-a+x0)/2+x1((1+a(-3+2a)+(4-4a+2x0)x0)/6+
  170.                 x1(1+a(-6+a(11-6a))+x0(11+a(-29+18a)+
  171.                 x0(18-18a+6x0)))/24));
  172.             x1 *= x0 (1+SetPrecision[x2,mp[j]]);
  173.             x0 += SetPrecision[x1,mp[prec]]
  174.             ];
  175.         x0
  176.         ];
  177.     
  178. invbet[x_,pp_,a_,b_] :=
  179.     Module[{prec=Min[precis[a],precis[b],precis[x],precis[pp]],
  180.         p,pc,x0,x1,x2,betx,betab,c=a/(a+b),cflag,ar,br,it,xa,xb},
  181.         If[prec == $MachinePrecision, prec--];
  182.         If[SameQ[prec,Infinity], Return[$Failed]];
  183.         betab = N[Beta[a,b],prec];
  184.         If[cflag = (x>c),
  185.             betx = Beta[N[1-x,prec],b,a]/betab;
  186.             p = pp+1 - betx;
  187.             pc = betx - pp,
  188.             (* else *)
  189.             betx = Beta[N[x,prec],a,b]/betab;
  190.             p = pp + betx;
  191.             pc = 1-betx-pp,
  192.             ];
  193.     (* p == BetaRegularized[0,x0,a,b] where x0 is still unknown
  194.       pc == BetaRegularized[x0,1,a,b] *)
  195.         If[p==0, Return[p]];
  196.         If[pc==0, Return[1-pc]];
  197.         If[p<0 || pc<0, Return[$Failed]];
  198.         If[a==1, Return[1-pc^(1/b)]];
  199.         If[b==1, Return[p^(1/a)]];
  200.  
  201.         xa = N[a p betab]^(1/a);
  202.         x1 = (b-1.)/(a+1.);
  203.         x2 = x1(-4+a(-1+a)+b(5+3a))/(2(a+1)(a+2));
  204.         xa *= 1+xa(x1+xa(x2+xa x1(18+a(4+a(-10+a(-1+a)))+
  205.             b(-47+a(-30+a(11+6a))+b(31+a(33+8a))))/
  206.             (3(a+1)^2(a+2)(a+3))));
  207.         If[xa < 0, xa = 10    (* just make it invalid. *)];
  208.  
  209.         xb = N[b pc betab]^(1/b);
  210.         x1 = (a-1.)/(b+1.);
  211.         x2 = x1(-4+b(-1+b)+a(5+3b))/(2(b+1)(b+2));
  212.         xb *= 1+xb(x1+xb(x2+xb x1(18+b(4+b(-10+b(-1+b)))+
  213.             a(-47+b(-30+b(11+6b))+a(31+b(33+8b))))/
  214.             (3(b+1)^2(b+2)(b+3))));
  215.         If[xb < 0, xb = 10    (* just make it invalid *)];
  216.  
  217.         If[a<1 || b<1 || Min[xa,xb]<0.1,
  218.             If[xa < xb,    (* choose the better initial value. *)
  219.                 x0 = xa; ar = a; br = b; cflag = False,
  220.                 (* else *)
  221.                 x0 = xb; ar = b; br = a; p = pc; cflag = True],
  222.             (* else *)
  223.             (* try AS 26.5.22 since both xa and xb may be wrong. *)
  224.             x1 = Sqrt[2.] If[p < pc, inverf[Infinity,-2p],
  225.                         inverf[-Infinity,2pc]];
  226.             x2 = 2/(1/(2a-1)+1/(2b-1));
  227.             ar = (x1^2-3)/6;
  228.             br = (1/(2a-1)-1/(2b-1))(ar+5/6-2/(3x2));
  229.             br += x1 Sqrt[x2+ar]/x2;
  230.             If[br > 0,
  231.                 x0 = a/(a+b Exp[2br]);
  232.                 ar = a; br = b; cflag = False,
  233.                 (* else *)
  234.                 x0 = b/(b+a Exp[-2br]);
  235.                 ar = b; br = a; p = pc; cflag = True]
  236.                 ];
  237.         p *= betab;
  238.  
  239.         (* use binary search to improve it. *)
  240.         betx = p-Beta[x0,ar,br];
  241.         If[betx > 0,
  242.             xb = If[x0 < .25, 2x0, (1+2x0)/3];
  243.             x2 = p-Beta[xb,ar,br];
  244.             While[x2 > 0,
  245.                 x0 = xb;
  246.                 betx = x2;
  247.                 xb = If[x0 < .25, 2x0, (1+2x0)/3];
  248.                 x2 = p-Beta[xb,ar,br]
  249.                 ];
  250.             xa = x0,
  251.             (* else *)
  252.             xa = If[x0 < .75, 2x0/3, 2x0-1];
  253.             x2 = p-Beta[xa,ar,br];
  254.             While[x2 < 0,
  255.                 x0 = xa;
  256.                 betx = x2;
  257.                 xa = If[x0 < .75, 2x0/3, 2x0-1];
  258.                 x2 = p-Beta[xa,ar,br]
  259.                 ];
  260.             xb = x0
  261.             ];
  262.  
  263.         j = 0;
  264.         While[j < prec,
  265.             If[j > 0 && NumberQ[x0],
  266.                 j = Max[j,Min[3j,
  267.                     Max[1,If[x1 == 0,
  268.                         Max[j*3,4],
  269.                         (* else *)
  270.                         Round[-1.3 Log[Abs[N[x1]]]]
  271.                         ]
  272.                     ]]];
  273.                 If[j > prec, j = If[j > 2prec, prec, prec-1]];
  274.                 x0 = SetPrecision[x0,mp[j]],
  275.                 (* else *)
  276.                 j = 1
  277.                 ];
  278.             dy = p-Beta[x0,ar,br];
  279.             While[dy == 0 && j < prec,
  280.                 (* adding a fake zero would corrupt everything. *)
  281.                 If[j < $MachinePrecision,
  282.                     j = 2 $MachinePrecision, j *= 2];
  283.                 If[j > prec, j = prec];
  284.                 x0 = SetPrecision[x0,mp[j]];
  285.                 dy = p-Beta[x0,ar,br]
  286.                 ];
  287.             If[dy == 0, Return[If[cflag, 1-x0, x0]]];
  288.             dy = SetPrecision[dy,mp[j]];
  289.             x1 = x0^-ar (1-x0)^-br dy;
  290.             x2 = x1((1-ar+x0(ar+br-2))/2+x1(1+ar(-3+2ar)+
  291.                 x0(-6+ar(10-4ar)+4br(1-ar)+
  292.                 x0(6+ar(-7+2ar)+br(-7+4ar+2br))))/6);
  293.             x1 *= x0 (1-x0) (1+x2);
  294.             x0 += x1
  295.             ];
  296.  
  297.         If[cflag, 1-x0, x0]
  298.         ];
  299.  
  300. inverf[xx_,pp_] :=
  301.     Module[{prec=Min[precis[xx],precis[pp]],
  302.         p,x=Abs[xx],x0,sgnxx,erfxx,cflag,rflag,i=16/3},
  303.         If[prec == $MachinePrecision, prec--];
  304.         If[SameQ[prec,Infinity], Return[$Failed]];
  305.         If[SameQ[Head[xx],DirectedInfinity],
  306.             cflag = True;
  307.             erfxx = 0;
  308.             sgnxx = Sign[xx[[1]]];
  309.             p = pp+sgnxx,
  310.             (* else *)
  311.             sgnxx = Sign[xx];
  312.             If[cflag = (x>1/2),
  313.                 erfxx = Erf[N[Abs[xx],prec],Infinity];
  314.                 p = pp+sgnxx(1-erfxx),
  315.                 (* else *)
  316.                 erfxx = Erf[N[Abs[xx],prec]];
  317.                 p = pp+sgnxx erfxx
  318.                 ]
  319.             ];
  320.         If[p==0, Return[p]];
  321.         If[rflag = (p<0), p = -p];
  322.         If[p<.9,
  323.             x0 = (0.8862244286 + (-1.220051315 +
  324.                 (0.2194688168 + 0.1329123619p)p)p)p/
  325.                 (1 - (1.37690286 + (0.01110212605 +
  326.                 (-0.4939796583 + 0.0964477892p)p)p)p);
  327.             If[cflag = (x0 > 1/2), p = 1-p],
  328.             (* else *)
  329.             (* p = 1-p , with roundoff error minimized. *)
  330.             If[cflag,
  331.                 If[rflag,
  332.                     (* p = 1 + (pp+sgnxx(1-erfxx)) *)
  333.                     p = pp-sgnxx erfxx;
  334.                     If[sgnxx==1, p += 2],
  335.                     (* else *)
  336.                     (* p = 1 - (pp+sgnxx(1-erfxx)) *)
  337.                     p = sgnxx erfxx-pp;
  338.                     If[sgnxx==-1, p += 2]
  339.                     ],
  340.                 (* else *)
  341.                 p = 1-p
  342.                 ];
  343.             cflag = True;
  344.             If[p<0, Return[$Failed]];
  345.             If[p==0, Return[If[rflag,-Infinity,Infinity]]];
  346.             x0 = N[Sqrt[(-0.555075007823429663  + p(
  347.                     -22943097299079.457 + p(
  348.                     -3.90911954269625125*10^21 + p(
  349.                     -4.20106086583717546*10^26 + p(
  350.                     -2.42091291012776922*10^29 + p(
  351.                      7.8542986097983114*10^29))))))/
  352.                 (1 + p(41894651750846.3689 + p(
  353.                      7.23863012895805264*10^21 + p(
  354.                      7.9066722226489883*10^26 + p(
  355.                      4.63843915091736701*10^29 + p(
  356.                     -1.56919477544184394*10^30))))))
  357.                 -Log[p]-Log[-Log[p]]/2]]
  358.             ];
  359.  
  360.         erfxx = p;
  361.         p = N[Sqrt[Pi],prec]/2;
  362.         i=0;
  363.         While[i<prec,
  364.             If[i > 0 && NumberQ[x0],
  365.                 i = Min[3i,Max[1,If[x == 0,
  366.                         Max[i*3,10],
  367.                         (* else *)
  368.                         Round[-8 Log[Abs[N[x]]]]
  369.                         ]
  370.                     ]];
  371.                 If[i > prec, i = If[i > 2prec, prec, prec-1]];
  372.                 x0 = SetPrecision[x0,mp[i]],
  373.                 (* else *)
  374.                 i = 1
  375.                 ];
  376.             x = If[cflag,
  377.                 (erfxx-Erf[x0,Infinity])p Exp[x0^2],
  378.                 (* else *)
  379.                 (Erf[x0]-erfxx)p Exp[x0^2]
  380.                 ];
  381.             x = SetPrecision[x,mp[i]];
  382.             x *= 1 - x x0/(1 + 2x x0);
  383.             x0 -= x
  384.             ];
  385.  
  386.         If[rflag, -x0, x0]
  387.         ];
  388.  
  389. End[]  (* Statistics`InverseStatisticalFunctions`Private *)
  390.  
  391. SetAttributes[ InverseErf ,ReadProtected];
  392. SetAttributes[ InverseErfc ,ReadProtected];
  393. SetAttributes[ InverseGammaRegularized ,ReadProtected];
  394. SetAttributes[ InverseBetaRegularized ,ReadProtected];
  395.  
  396. Protect[InverseErf, InverseErfc, InverseGammaRegularized, 
  397.     InverseBetaRegularized];
  398.  
  399. EndPackage[ ] (* Statistics`InverseStatisticalFunctions` *)
  400.  
  401. (*:Limitations:
  402.     The inverse functions defined only work for real arguments within
  403.     ranges which give real results.
  404. *)
  405.  
  406.  
  407.