MultipleR2.m
BeginPackage["MultipleR2`"]
Needs["Statistics`ContinuousDistributions`"]
probr::usage=
  "probr[n,p,rho2,x]computes the probability that \!\(R^2\) is less 
	than x ,under the random regression model."
confidence::usage=
  "confidence[n,p,r2,plevel]computes the lower or upper bound for a plevel \
confidence interval for 
		\[Rho]^2 given an observed r^2 under the random regression model."
confidencefix::usage=
  "confidencefix[n,p,r2,plevel]computes the lower or upper bound for a plevel \
confidence interval for \[Rho]^2 given an observed  r^2 under the fixed \
model."			
probf::usage="probf[n,p,fval]computes  the probability that
	the F is less than or equal to fval."
powerR2::usage=
  "powerR2[n,p,rho2,alpha] computes the probability  of rejecting  the \
hypothesis that  \!\(\[Rho]^2\)= 0  when  \!\(\[Rho]^2\)= rho2."		
powerR2fix::usage=
  "powerR2fix[n,p,rho2,alpha] computes the probability  of rejecting  the \
hypothesis that  \!\(\[Rho]^2\)= 0  when  \!\(\[Rho]^2\)= rho2."		
invprobf::usage="invprobf[n,p,upperp]is the inverse function for probf
	returning the fval."
stdnormal::usage="stdnormal[z] computes the probability a standard normal 
		variable is less than z."
findsample::usage=
  "findsample[p,rho2,power,alpha] for the given power and alpha it computes \
the sample size necessary to reject the hypothesis that  \!\(\[Rho]^2\)= 0  \
when  \!\(\[Rho]^2\)= rho2, under the random regression model."
findsamplefix::usage=
  "findsamplefix[p,rho2,power,alpha] for the given power and alpha it \
computes the sample size necessary to reject the hypothesis that  \
\!\(\[Rho]^2\)= 0  when  \!\(\[Rho]^2\)= rho2, under the random fixed model."	\
		
invnormal::usage="invnormal[plevel]is the inverse of stdnormal."
r2cross::usage=
  "r2cross[n,p,r2,plevel,pcross] using a plevel one-sided confidence for \!\(\
\[Rho]^2\) and the observed r2, it establishes a  one-sided confidence \
interval for the cross-validated  R^2 at the pcross level."	
crossR2confidence::usage=
  "crossR2confidence[n,p,rho2,pcross] for a given  \!\(\[Rho]^2\) (rho2) it \
returns a one-sided confidence interval for the cross-validated \!\(R^2\) at \
the pcross level."
		crossR2sample::usage=
  "crossR2sample[p,rho2,pcross,diff] it gives the necessary sample to obtain \
a cross-validated  \!\(R^2\) that does not differ from rho2 by more than diff \
with probability pcross."
		probfnoncentral::usage=
  "probfnoncentral[x, dfn, dfd, non] it integrates the noncentral \
f-distribution with degrees of freedom dfn and dfd from 0 to x.  The non is \
the noncentrality parameter." 
Begin["`Private`"]

	stdnormal[z_]:=Block[{f,prob,zabs},f=Exp[-.5 x^2]/Sqrt[2*Pi];
	zabs=Abs[z];prob=NIntegrate[f,{x,0,zabs}];
		If[z\[GreaterEqual]0,prob+=.5];If[z<0,prob=.5-prob];
	Return[prob]]

invnormal[plevel_]:=Block[{z,zupper,zlower,k,i,pcounter},
		If[plevel<.5000,pcounter=N[1.0-plevel],
			pcounter=plevel];i=0;z=0.0;zupper=3.5;
		zlower=0.0;k=.5000;
		While[Abs[pcounter-stdnormal[z]]>.00001,
			z=z+k (zupper-zlower);
			If[pcounter-stdnormal[z]<0,{k=-.5000,zupper=z},{k=.5000,zlower=z}];
					i+=1;
			If[i>150,Break[]]];
		If[plevel<.5000,z=-z];
		Return[z]]
confidence[n_,p_,r2_,plevel_]:=
		Block[{rho2,a,b,f,fval,int,RR2,lower,upper,k,ssum,ck,cbeta,prob,x,kint,e},
		f=FRatioDistribution[p,n-p-1];
		fval=Quantile[f,plevel];
		int=(p*fval)/(n-p-1);
		RR2=int/(1+int);
		If[RR2>=r2,Return[0.0]];
		rho2=r2; x=r2;
		a=(n-1)/2; b=(n-p-1)/2;
		lower=0.000; upper=1.000;prob=plevel-.1;
		While[Abs[prob-plevel]>.001,
						cbeta=(1-rho2)^a/Beta[1,p/2,b];
				k=0;ck=3;ssum=0.0;
		While[ck>=.0000000000001 ,
					k+=1;
				e=rho2^k;
						ck=((((Gamma[a+k]/Gamma[a])^2)*e)/((Gamma[p/2+k]/Gamma[p/2])*
						Gamma[k+1]))*Beta[x,p/2+k,b];
					cp=ck;
			If[cp <.0000000000000001, ck=.00000000000001];
			If[k > 2*p, ck=cp];
						ssum=ssum+ck;];
		prob=(Beta[x,p/2,b]+ssum)*cbeta;
		If[prob>plevel, {kint=.5,lower=rho2}, {kint=-.5,upper=rho2}];
			rho2+=kint*(upper-lower)];
		Return[rho2-=kint*(upper-lower)]]		

probr[n_,p_,rho2_,x_]:=0 /; n-p-1<3
probr[n_,p_,rho2_,x_]:=
	Block[{a,b,cbeta,k,ck,ssum,prob},
					a=(n-1)/2; b=(n-p-1)/2;
				cbeta=(1-rho2)^a/Beta[1,p/2,b];
				k=0;ck=3;ssum=0.0;
		While[ck>=.0000000000001 ,
					k+=1;
						ck=(((Gamma[a+k]/Gamma[a])^2*rho2^k)/((Gamma[p/2+k]/Gamma[p/2])*
						Gamma[k+1]))*Beta[x,p/2+k,b];
			cp=ck;
			If[cp <.0000000000000001, ck=.00000000000001];
			If[k > 2*p, ck=cp];
				ssum=ssum+ck;];
				prob=(Beta[x,p/2,b]+ssum)*cbeta;
		Return[prob]
		 ]  /; n-p-1>=3

probf[n_,p_,fval_]:= 
		Block[{v1,v2,k1,prob},v1=p;
		v2=n-p-1;
		k1=N[v1*fval/(v2+v1*fval)];
		prob=BetaRegularized[k1,v1/2,v2/2];
		Return[prob]]

invprobf[n_,p_,upperp_]:= 
		Block[{v1,v2,bval,fval},v1=p;
		v2=n-p-1;
		bval=InverseBetaRegularized[upperp,v1/2,v2/2];
		fval=bval*v2/(v1-bval*v1);
			Return[fval] ]
powerR2[n_,p_,rho2_,alpha_:.05]:=
	Block[{f,pp,R2,fval,int,pval},
		f=FRatioDistribution[p,n-p-1];
			pp=N[1-alpha];
				fval=Quantile[f,pp];
		int=(p*fval)/(n-p-1);
		R2=int/(1+int);
		pval=1-probr[n,p,rho2,R2];
Return[pval] ]	
findsamplefix[p_,rho2_,power_,alpha_:.05]:=
	Block[{pval,n,probx,lambda,f,fcrit,dfe,nupper},
		pval=1-alpha;
		n=p+5;
		probx=0;
		While[power-probx>0,
			dfe=n-p-1;
			f=FRatioDistribution[p,dfe];
			fcrit=Quantile[f,pval];
			lambda=(rho2/(1-rho2))*(n-1);
			probx=1-probfnoncentral[fcrit,p,dfe,lambda];
			n+=50];
		nupper=n-50;
		n-=100;
		While[n<nupper,
			dfe=n-p-1;
			f=FRatioDistribution[p,dfe];
			fcrit=Quantile[f,pval];
			lambda=(rho2/(1-rho2))*(n-1);
			fnon=NoncentralFRatioDistribution[p,dfe,lambda];
			probx=1-CDF[fnon,fcrit];
			++n;
		If[probx-power>-.001, Break[]]
					];
		Return[n-1] ]
findsample[p_,rho2_,power_,alpha_:.05]:=
	Block[{pval,n,probx,f,fcrit,dfe,nupper,rr,r2},
		pval=1-alpha;
		n=p+5;
		probx=0;
		While[power-probx>0,
			dfe=n-p-1;
			f=FRatioDistribution[p,dfe];
			fcrit=Quantile[f,pval];
			rr=(p*fcrit)/(n-p-1);
			r2=rr/(1+rr);
			probx=1-probr[n,p,rho2,r2];
			n+=50];
		nupper=n-50;
		n-=100;
		While[n<nupper,
			dfe=n-p-1;
			f=FRatioDistribution[p,dfe];
			fcrit=Quantile[f,pval];
			rr=(p*fcrit)/(n-p-1);
			r2=rr/(1+rr);
			probx=1-probr[n,p,rho2,r2];
			++n;
		If[probx-power>-.001, Break[]]
					];
		Return[n-1] ]
powerR2fix[n_,p_,rho2_,alpha_:.05]:=
	Block[{f,pp,fval,dfe,lambda,pval},
		f=FRatioDistribution[p,n-p-1];
		pp=1-alpha;
		fval=Quantile[f,pp];
		dfe=n-p-1;
		lambda=(rho2/(1-rho2))*(n-1);
		pval=1-probfnoncentral[fval,p,dfe,lambda];
		Return[pval] ]	
confidencefix[n_,p_,r2_,plevel_]:=
	Block[{rho2,lower,upper,prob,fcrit,kint,lambda,f,fval,int,RR2},
			f=FRatioDistribution[p,n-p-1];
		fval=Quantile[f,plevel];
		int=(p*fval)/(n-p-1);
		RR2=int/(1+int);
		If[RR2>=r2,Return[0.0]];
		rho2=r2;
		lower=0;
		upper=1.000;
		prob=plevel-.1;
		fcrit=(r2/(1-r2))*((n-p-1)/p);
		While[Abs[prob-plevel]>.001,
		lambda=(rho2/(1-rho2))*(n-1);
		prob=probfnoncentral[fcrit,p,n-p-1,lambda];
	  If[prob>plevel,{kint=.5,lower=rho2},{kint=-.5,upper=rho2}];
				rho2+=kint*(upper-lower)
										];
				Return[rho2-=kint*(upper-lower)]]
r2cross[n_,p_,r2_,plevel_,pcross_]:=
	Block[{lambda,c,f,rho2low,rho2,a,b,lower,upper,prob,probf,rho2cross,k,ck,
      ssum,x,kint},
		rho2=r2/2; x=r2;
		a=(n-1)/2; b=(n-p-1)/2;
		lower=0.000; upper=r2;prob=plevel-.1;
		While[Abs[prob-plevel]>.001,
						cbeta=(1-rho2)^a/Beta[1,p/2,b];
				k=0;ck=3;ssum=0.0;
		While[ck>=.0000000000001 ,
					k+=1;
						ck=(((Gamma[a+k]/Gamma[a])^2*rho2^k)/((Gamma[p/2+k]/Gamma[p/2])*
						Gamma[k+1]))*Beta[x,p/2+k,b];
				ssum=ssum+ck;];
		prob=(Beta[x,p/2,b]+ssum)*cbeta;
		If[prob>plevel, {kint=.5,lower=rho2}, {kint=-.5,upper=rho2}];
			rho2+=kint*(upper-lower)];
	rho2low=rho2-kint*(upper-lower);	
		(* computation of the cross interval *)
		lambda=((n-p-1)*rho2low)/(1-rho2low);
		k=.5; upper=rho2low;lower=0;rho2cross=0;probf=0;
		While[Abs[pcross-probf]>.001,
			rho2cross+=k*(upper-lower);
			f=NoncentralFRatioDistribution[1,p-1,lambda];
		c=(p-1)/(rho2low/rho2cross-1);
		probf=1-CDF[f,c];
			If[pcross-probf>0,{k=-.5,upper=rho2cross},{k=.5,lower=rho2cross}]
					
	];
	  {rho2low,prob,rho2cross,probf} ]
crossR2confidence[n_,p_,rho2_,pcross_]:=
	Block[{lambda,k,upper,lower,rho2cross,probf,f,c  },
		lambda=((n-p-1)*rho2)/(1-rho2);
		k=.5; upper=rho2;lower=0;rho2cross=0;probf=0;
		While[Abs[pcross-probf]>.001,
			rho2cross+=k*(upper-lower);
			f=NoncentralFRatioDistribution[1,p-1,lambda];
		c=(p-1)/(rho2/rho2cross-1);
		probf=1-CDF[f,c];
			If[pcross-probf>0,{k=-.5,upper=rho2cross},{k=.5,lower=rho2cross}]
					
	];
	  {rho2cross,probf} ]
crossR2sample[p_,rho2_,pcross_,diff_]:=
	Block[{pval,n,probxf,rho2cross,lambda,nupper},
		rho2cross=rho2-diff;
		n=p+5;
		probx=0;
		While[pcross-probx>0,
			lambda=((n-p-1)*rho2)/(1-rho2);
			f=NoncentralFRatioDistribution[1,p-1,lambda];
			c=(p-1)/(rho2/rho2cross-1);
			probx=1-CDF[f,c];
			n+=50];
		nupper=n-50;
		n-=100;
		While[n<nupper,
			lambda=((n-p-1)*rho2)/(1-rho2);
			f=NoncentralFRatioDistribution[1,p-1,lambda];
			c=(p-1)/(rho2/rho2cross-1);
			probx=1-CDF[f,c];
			++n;
			If[probx-pcross>-.001, Break[]]
		];
		{probx,n-1}]
\!\(probfnoncentral[x_Real, dfn_, dfd_, non_Real] := 
    Block[{j, ssum, v1, xx, v2, c, b}, \nj = 0; ssum = 0.0; v1 = dfn\/2; 
      v2 = dfd\/2; \nxx = \(x\ dfn\)\/\(x\ dfn + dfd\); c = 1; b = 1; \n
      While[c\ b >  .0000000000001, \n\t\t\tj += 1; 
        c = \((non\/2)\)\^j\/\(j!\); \n\t\t\t
        b = BetaRegularized[xx, v1 + j, v2]; \n\t\t\tssum = ssum + c\ b; ]; \n
      \t\tssum = 
        Exp[\(-\(non\/2\)\)]\ \((ssum + BetaRegularized[xx, v1, v2])\); \n\t\t
      Return[ssum]]\)
End[]
Protect[probr,confidence, confidencefix, probf,powerR2,powerR2fix, invprobf,
  stdnormal, findsample, findsamplefix,invnormal,crossR2sample,
  crossR2confidence,probfnoncentral]
EndPackage[]
