
(*:Version: Mathematica 2.0 *)

(*:Name: Statistics`InverseStatisticalFunctions` *)

(*:Title: Evaluating the Inverse of Erf[ ], GammaRegularized[ ], and
        BetaRegularized[ ] *)

(*:Author:
  Jerry B. Keiper (Wolfram Research, Inc.), March, 1989.
  Current version April, 1991.
*)

(*:Summary:
  Defines the inverses of the following statistical functions:
  Erf[ ], Erfc[ ], GammaRegularized[ ], and BetaRegularized[ ].
  The algorithms are similar to Newton's method except that, rather
  than using a linear approximation, higher-order approximations are
  used leading to cubic, and in the case of InverseErf[ ], quartic
  convergence.
*)

(*:Keywords: inverse function, Newton's method, InverseErf[ ],
  InverseErfc[ ], InverseGammaRegularized[ ], InverseBetaRegularized[ ]
*)

(*:Requirements: no special system requirements. *)

(*:Warnings: none. *)

(*:Sources:
  M. Abramowitz & I. Stegun, Handbook of Mathematical Functions,
    Dover, New York, 1972, chapter 26
*)

BeginPackage["Statistics`InverseStatisticalFunctions`"]

InverseErf::usage = 
"InverseErf[p] gives the value of x such that p == Erf[x].  InverseErf[x0,p]
gives the value of x such that p == Erf[x0,x].  InverseErf[x0,-p] gives the
value of x such that p == Erf[x,x0].  x0 and p must be real, and p must be
between Erf[x0, -Infinity] and Erf[x0, Infinity]"

InverseErfc::usage =
"InverseErfc[p] gives the value of x such that p == Erfc[x].  p must satisfy
0 <= p <= 2."

InverseGammaRegularized::usage =
"InverseGammaRegularized[a,x0,p] gives the value of x such that p ==
GammaRegularized[a,x0,x].  InverseGammaRegularized[a,x0,-p] gives the value
of x such that p == GammaRegularized[a,x,x0].  a, x0, and p must be real,
a > 0, 0 <= x <= Infinity, and GammaRegularized[a,x,0] <= p <=
GammaRegularized[a,x,Infinity]."

InverseBetaRegularized::usage =
"InverseBetaRegularized[x0,p,a,b] gives the value of x such that p ==
BetaRegularized[x0,x,a,b].  InverseBetaRegularized[x0,-p,a,b] gives the
value of x such that p == BetaRegularized[x,x0,a,b].  a, b, x0, and p must
be real, a > 0, b > 0, 0 <= x <= 1, and BetaRegularized[x,0,a,b] <= p <=
BetaRegularized[x,1,a,b]."

Begin["Statistics`InverseStatisticalFunctions`Private`"]

InverseErf[p_] := InverseErf[0,p];
InverseErfc[p_] := InverseErf[Infinity,-p];
InverseErf[x_, 0] := x;

InverseErf[x_,p_?NumberQ] :=
    Module[{e},
	e /; (e = inverf[x,p]) =!= $Failed
	] /; (Im[p]==0 && ((NumberQ[x] && Im[x]==0) ||
		SameQ[x,Infinity] || SameQ[x,-Infinity]))

InverseGammaRegularized[a_?NumberQ,x_,p_?NumberQ] :=
    Module[{e},
	e /; (e = invgam[a,x,p]) =!= $Failed
	] /; (Im[p]==0 && Im[a]==0 && a>0 && (SameQ[x,Infinity] ||
			(NumberQ[x] && Im[x]==0 && x >= 0)))

InverseBetaRegularized[x_?NumberQ,p_?NumberQ,a_?NumberQ,b_?NumberQ] :=
    Module[{e},
	e /; (e = invbet[x,p,a,b]) =!= $Failed
	] /; (Im[x]==0 && x >= 0 && x <= 1 && Im[p]==0 &&
			Im[a]==0 && a>0 && b>0 && Im[b]==0) 

mp[j_] := Max[j, $MachinePrecision]
precis[x_] := If[MachineNumberQ[x], $MachinePrecision, Precision[x]]

invgam[a_,x_,pp_] :=
    Module[{prec=Min[precis[a],precis[x],precis[pp]],p,x0,
		gamx,gama,cflag,j,pa,temp,app,dy,x1,x2},
		If[prec == $MachinePrecision, prec--];
		If[SameQ[prec,Infinity], Return[$Failed]];
		gama = N[Gamma[a],prec];
		If[SameQ[Head[x],DirectedInfinity],
			cflag = True;
			gamx = 0;
			p = pp+1,
		    (* else *)
			If[cflag = (x>a),
				gamx = Gamma[a,N[x,prec],Infinity]/gama;
				p = pp+1 - gamx,
			    (* else *)
				gamx = Gamma[a,0,N[x,prec]]/gama;
				p = pp + gamx
				]
			];
		(* p == GammaRegularized[a,0,x0] where x0 is still unknown
		gamx is the value of the incomplete gamma function from
		0 to x or x to Infinity depending on cflag. *)
		If[p==0, Return[p]];
		If[p<0 || p>1, Return[$Failed]];
		If[a == 1,
			p = If[cflag, gamx - pp, 1 - p];
			Return[-N[Log[p],prec]]];
		app = N[a];
		pa = (p a gama)^(1/a);
		temp = pa/(app+1);
		x0 = pa(1+temp(1+temp((5+3app)/(2(app+2))+
			temp((31+app(33+8app))/(3(2+app)(3+app))+temp(
			(2888+app(5661+app(3971+app(1179+125app))))/
				(24(2+app)^2(3+app)(4+app))+
			temp((9074+app(20997+app(18375+app(7575+app(
				1471+108app)))))/
			(10(2+app)^2(3+app)(4+app)(5+app))))))));
		If[x0 > 5, (* if x0 is large it may be very wrong *)
			If[cflag,
				temp = inverf[Infinity,2(pp - gamx)],
			    (* else *)
				temp = inverf[-Infinity,2 N[p]]];
			x1 = app (1 - 1/(9app) + temp/Sqrt[4.5 app])^3;
			(* if x1 is small it may be very wrong *)
			If[x1 > 5, x0 = x1]];
		If[x0 > app,
			p = If[cflag, gamx - pp, cflag = True; 1 - p],
		    (* else *)
			cflag = False];
		If[p<0, Return[$Failed]];
	(* now use Taylor series to quartically converge to x0. *)
		p *= gama;
		j = 0;
		While[j < prec,
			If[j > 0 && NumberQ[x0],
				j = Max[j,Min[4j,
					Max[1,If[x1 == 0,
						Max[4j,4],
					    (* else *)
						Round[-1.7 Log[Abs[N[x1]]]]
					]]]];
				If[j > prec, j = If[j > 2prec, prec, prec-1]];
				x0 = SetPrecision[x0,mp[j]],
			    (* else *)
				j = 1
				];
			dy = p-If[cflag,Gamma[a,x0,Infinity],Gamma[a,0,x0]];
			While[dy == 0 && j < prec,
			    (* adding a fake zero would corrupt everything *)
				If[j < $MachinePrecision,
					j = 2 $MachinePrecision, j *= 2];
				If[j > prec, j = prec];
				x0 = SetPrecision[x0,mp[j]];
				dy = p-If[cflag,
					Gamma[a,x0,Infinity],
					Gamma[a,0,x0]]
				];
			If[dy == 0, Return[x0]];
			dy = SetPrecision[If[cflag, -dy, dy],mp[j]];
			x1 = Exp[x0] x0^-a dy;
			x2 = x1((1-a+x0)/2+x1((1+a(-3+2a)+(4-4a+2x0)x0)/6+
				x1(1+a(-6+a(11-6a))+x0(11+a(-29+18a)+
				x0(18-18a+6x0)))/24));
			x1 *= x0 (1+SetPrecision[x2,mp[j]]);
			x0 += SetPrecision[x1,mp[prec]]
			];
		x0
		];
	
invbet[x_,pp_,a_,b_] :=
    Module[{prec=Min[precis[a],precis[b],precis[x],precis[pp]],
		p,pc,x0,x1,x2,betx,betab,c=a/(a+b),cflag,ar,br,it,xa,xb},
		If[prec == $MachinePrecision, prec--];
		If[SameQ[prec,Infinity], Return[$Failed]];
		betab = N[Beta[a,b],prec];
		If[cflag = (x>c),
			betx = Beta[N[1-x,prec],b,a]/betab;
			p = pp+1 - betx;
			pc = betx - pp,
		    (* else *)
			betx = Beta[N[x,prec],a,b]/betab;
			p = pp + betx;
			pc = 1-betx-pp,
			];
	(* p == BetaRegularized[0,x0,a,b] where x0 is still unknown
	  pc == BetaRegularized[x0,1,a,b] *)
		If[p==0, Return[p]];
		If[pc==0, Return[1-pc]];
		If[p<0 || pc<0, Return[$Failed]];
		If[a==1, Return[1-pc^(1/b)]];
		If[b==1, Return[p^(1/a)]];

		xa = N[a p betab]^(1/a);
		x1 = (b-1.)/(a+1.);
		x2 = x1(-4+a(-1+a)+b(5+3a))/(2(a+1)(a+2));
		xa *= 1+xa(x1+xa(x2+xa x1(18+a(4+a(-10+a(-1+a)))+
			b(-47+a(-30+a(11+6a))+b(31+a(33+8a))))/
			(3(a+1)^2(a+2)(a+3))));
		If[xa < 0, xa = 10	(* just make it invalid. *)];

		xb = N[b pc betab]^(1/b);
		x1 = (a-1.)/(b+1.);
		x2 = x1(-4+b(-1+b)+a(5+3b))/(2(b+1)(b+2));
		xb *= 1+xb(x1+xb(x2+xb x1(18+b(4+b(-10+b(-1+b)))+
			a(-47+b(-30+b(11+6b))+a(31+b(33+8b))))/
			(3(b+1)^2(b+2)(b+3))));
		If[xb < 0, xb = 10	(* just make it invalid *)];

		If[a<1 || b<1 || Min[xa,xb]<0.1,
			If[xa < xb,	(* choose the better initial value. *)
				x0 = xa; ar = a; br = b; cflag = False,
			    (* else *)
				x0 = xb; ar = b; br = a; p = pc; cflag = True],
		    (* else *)
			(* try AS 26.5.22 since both xa and xb may be wrong. *)
			x1 = Sqrt[2.] If[p < pc, inverf[Infinity,-2p],
						inverf[-Infinity,2pc]];
			x2 = 2/(1/(2a-1)+1/(2b-1));
			ar = (x1^2-3)/6;
			br = (1/(2a-1)-1/(2b-1))(ar+5/6-2/(3x2));
			br += x1 Sqrt[x2+ar]/x2;
			If[br > 0,
				x0 = a/(a+b Exp[2br]);
				ar = a; br = b; cflag = False,
			    (* else *)
				x0 = b/(b+a Exp[-2br]);
				ar = b; br = a; p = pc; cflag = True]
				];
		p *= betab;

		(* use binary search to improve it. *)
		betx = p-Beta[x0,ar,br];
		If[betx > 0,
			xb = If[x0 < .25, 2x0, (1+2x0)/3];
			x2 = p-Beta[xb,ar,br];
			While[x2 > 0,
				x0 = xb;
				betx = x2;
				xb = If[x0 < .25, 2x0, (1+2x0)/3];
				x2 = p-Beta[xb,ar,br]
				];
			xa = x0,
		    (* else *)
			xa = If[x0 < .75, 2x0/3, 2x0-1];
			x2 = p-Beta[xa,ar,br];
			While[x2 < 0,
				x0 = xa;
				betx = x2;
				xa = If[x0 < .75, 2x0/3, 2x0-1];
				x2 = p-Beta[xa,ar,br]
				];
			xb = x0
			];

		j = 0;
		While[j < prec,
			If[j > 0 && NumberQ[x0],
				j = Max[j,Min[3j,
					Max[1,If[x1 == 0,
						Max[j*3,4],
					    (* else *)
						Round[-1.3 Log[Abs[N[x1]]]]
						]
					]]];
				If[j > prec, j = If[j > 2prec, prec, prec-1]];
				x0 = SetPrecision[x0,mp[j]],
			    (* else *)
				j = 1
				];
			dy = p-Beta[x0,ar,br];
			While[dy == 0 && j < prec,
			    (* adding a fake zero would corrupt everything. *)
				If[j < $MachinePrecision,
					j = 2 $MachinePrecision, j *= 2];
				If[j > prec, j = prec];
				x0 = SetPrecision[x0,mp[j]];
				dy = p-Beta[x0,ar,br]
				];
			If[dy == 0, Return[If[cflag, 1-x0, x0]]];
			dy = SetPrecision[dy,mp[j]];
			x1 = x0^-ar (1-x0)^-br dy;
			x2 = x1((1-ar+x0(ar+br-2))/2+x1(1+ar(-3+2ar)+
				x0(-6+ar(10-4ar)+4br(1-ar)+
				x0(6+ar(-7+2ar)+br(-7+4ar+2br))))/6);
			x1 *= x0 (1-x0) (1+x2);
			x0 += x1
			];

		If[cflag, 1-x0, x0]
		];

inverf[xx_,pp_] :=
    Module[{prec=Min[precis[xx],precis[pp]],
		p,x=Abs[xx],x0,sgnxx,erfxx,cflag,rflag,i=16/3},
		If[prec == $MachinePrecision, prec--];
		If[SameQ[prec,Infinity], Return[$Failed]];
		If[SameQ[Head[xx],DirectedInfinity],
			cflag = True;
			erfxx = 0;
			sgnxx = Sign[xx[[1]]];
			p = pp+sgnxx,
		    (* else *)
			sgnxx = Sign[xx];
			If[cflag = (x>1/2),
				erfxx = Erf[N[Abs[xx],prec],Infinity];
				p = pp+sgnxx(1-erfxx),
			    (* else *)
				erfxx = Erf[N[Abs[xx],prec]];
				p = pp+sgnxx erfxx
				]
			];
		If[p==0, Return[p]];
		If[rflag = (p<0), p = -p];
		If[p<.9,
			x0 = (0.8862244286 + (-1.220051315 +
				(0.2194688168 + 0.1329123619p)p)p)p/
				(1 - (1.37690286 + (0.01110212605 +
				(-0.4939796583 + 0.0964477892p)p)p)p);
			If[cflag = (x0 > 1/2), p = 1-p],
		    (* else *)
			(* p = 1-p , with roundoff error minimized. *)
			If[cflag,
				If[rflag,
					(* p = 1 + (pp+sgnxx(1-erfxx)) *)
					p = pp-sgnxx erfxx;
					If[sgnxx==1, p += 2],
				    (* else *)
					(* p = 1 - (pp+sgnxx(1-erfxx)) *)
					p = sgnxx erfxx-pp;
					If[sgnxx==-1, p += 2]
					],
			    (* else *)
				p = 1-p
				];
			cflag = True;
			If[p<0, Return[$Failed]];
			If[p==0, Return[If[rflag,-Infinity,Infinity]]];
			x0 = N[Sqrt[(-0.555075007823429663  + p(
					-22943097299079.457 + p(
					-3.90911954269625125*10^21 + p(
					-4.20106086583717546*10^26 + p(
					-2.42091291012776922*10^29 + p(
					 7.8542986097983114*10^29))))))/
				(1 + p(41894651750846.3689 + p(
					 7.23863012895805264*10^21 + p(
					 7.9066722226489883*10^26 + p(
					 4.63843915091736701*10^29 + p(
					-1.56919477544184394*10^30))))))
				-Log[p]-Log[-Log[p]]/2]]
			];

		erfxx = p;
		p = N[Sqrt[Pi],prec]/2;
		i=0;
		While[i<prec,
			If[i > 0 && NumberQ[x0],
				i = Min[3i,Max[1,If[x == 0,
						Max[i*3,10],
					    (* else *)
						Round[-8 Log[Abs[N[x]]]]
						]
					]];
				If[i > prec, i = If[i > 2prec, prec, prec-1]];
				x0 = SetPrecision[x0,mp[i]],
			    (* else *)
				i = 1
				];
			x = If[cflag,
				(erfxx-Erf[x0,Infinity])p Exp[x0^2],
			    (* else *)
				(Erf[x0]-erfxx)p Exp[x0^2]
				];
			x = SetPrecision[x,mp[i]];
			x *= 1 - x x0/(1 + 2x x0);
			x0 -= x
			];

		If[rflag, -x0, x0]
		];

End[]  (* Statistics`InverseStatisticalFunctions`Private *)

SetAttributes[ InverseErf ,ReadProtected];
SetAttributes[ InverseErfc ,ReadProtected];
SetAttributes[ InverseGammaRegularized ,ReadProtected];
SetAttributes[ InverseBetaRegularized ,ReadProtected];

Protect[InverseErf, InverseErfc, InverseGammaRegularized, 
	InverseBetaRegularized];

EndPackage[ ] (* Statistics`InverseStatisticalFunctions` *)

(*:Limitations:
	The inverse functions defined only work for real arguments within
	ranges which give real results.
*)


