statfunc.cpp

可进行多种分布函数的计算

#include <cmath>

#include "statfunc.h"

#define NO_DATA -1000
#define MAX_FISHER 500

double CalcZ(double Q)
{
	//(0<Q <1.0)
	//Hastings偺嬤帡幃傪巊梡丅
	//妋棪Q偵懳墳偡傞惓婯暘晍偺z傪曉偡丅

	double z;

	const double a0=2.30753;
	const double a1=0.27061;
	const double b1=0.99229;
	const double b2=0.04481;

	if((Q<=0)||(Q>=1.0)){
		return NO_DATA;
	}
	else if(Q<=0.5){
		z=sqrt(log(1/Q/Q));
		return z-(a0+a1*z)/(1+b1*z+b2*z*z);
	}
	else{
		Q=1-Q;
		z=sqrt(log(1/Q/Q));
		return -(z-(a0+a1*z)/(1+b1*z+b2*z*z));
	}
}

double CalcPValue_Z(double u)
{
	//惓婯暘晍昞偺u偵懳墳偡傞妋棪傪曉偡丅
	//弌揟丗僷僜僐儞摑寁夝愅僴儞僪僽僢僋(婎慴摑寁曇) p.285丅

	double temp;

	const double d1=0.0498673470;
	const double d2=0.0211410061;
	const double d3=0.0032776263;
	const double d4=0.0000380036;
	const double d5=0.0000488906;
	const double d6=0.0000053830;

	if(u>0){
		temp=1+d1*u+d2*pow(u,2)+d3*pow(u,3)+d4*pow(u,4)+d5*pow(u,5)+d6*pow(u,6);
		return pow(temp,-16)/2;
	}
	else{
		u=-u;
		temp=1+d1*u+d2*pow(u,2)+d3*pow(u,3)+d4*pow(u,4)+d5*pow(u,5)+d6*pow(u,6);
		return 1-pow(temp,-16)/2;
	}

}

double CalcPValue_Kai(int f,double t)
{
	//t偵懳墳偡傞帺桼搙f偺僇僀俀忔暘晍偺妋棪傪寁嶼偡傞丅
	//弌揟丗僷僜僐儞摑寁夝愅僴儞僪僽僢僋(婎慴摑寁曇) p.286丅

	double ans,kai,temp,sum;
	int i,j,n;

	if(t<0){
		t=0;
	}

	if(f>40){
		temp=(pow(t/(double)f,1/3)-(1-2/9/(double)f))/sqrt(2/9/(double)f);
		ans=CalcPValue_Z(temp);
	}
	else if(f<=40){
		kai=sqrt(t);
		n=f%2;
		if(n==0){
			for(i=2,sum=1;i<=f-2;i+=2){
				for(j=2,temp=1;j<=i;j+=2){
					temp=temp*j;
				}
				sum=sum+pow(kai,i)/temp;
			}

			ans=exp(-t/2)*sum;
		}
		else{
			for(i=1,sum=0;i<=f-2;i+=2){
				for(j=1,temp=1;j<=i;j+=2){
					temp=temp*j;
				}
				sum=sum+pow(kai,i)/temp;
			}

			ans=2*CalcPValue_Z(kai)+sqrt(2/3.141592)*exp(-t/2)*sum;
		}
	}

	return ans;
}

double CalcKai(int f,double Q)
{
	//妋棪Q偵懳墳偡傞帺桼搙f偺僇僀俀忔暘晍偺抣傪寁嶼偡傞丅
	//弌揟丗僷僜僐儞摑寁夝愅僴儞僪僽僢僋(婎慴摑寁曇) p.286丅

	double temp,alpha1,alpha2,zL,zU,z0;
	int i;

	if(f==1){
		return CalcZ(Q/2)*CalcZ(Q/2);
	}
	else if(f==2){
		return -2*log(Q);
	}
	else if(f>=3){
		zL=0;
		zU=1;
		z0=0.5;
		alpha1=Q;
		temp=1/z0-1;
		alpha2=CalcPValue_Kai(f,temp);
		for(i=0;i<500;i++){
			if(alpha2<alpha1){
				zL=z0;
			}
			else{
				zU=z0;
			}
			z0=(zL+zU)/2;
			temp=1/z0-1;
			alpha2=CalcPValue_Kai(f,temp);
		}

		return 1/z0-1;
	}
	else{
		return NO_DATA;
	}

}


double CalcPValue_t(int f,double t)
{
	//t偵懳墳偡傞帺桼搙f偺倲暘晍昞偺妋棪傪寁嶼偡傞丅
	//弌揟丗僷僜僐儞摑寁夝愅僴儞僪僽僢僋(婎慴摑寁曇) p.289丅

	if(t>0){
		return 0.5*CalcPValue_F(1,f,t*t);
	}
	else{
		return 1-0.5*CalcPValue_F(1,f,t*t);
	}
}

double CalcPValue_F(int m,int n,double F)
{
	//F偵懳墳偡傞帺桼搙m,n偺俥暘晍昞偺妋棪傪寁嶼偡傞丅
	//弌揟丗僷僜僐儞摑寁夝愅僴儞僪僽僢僋(婎慴摑寁曇) p.287丅

	double x,Ix,U,f1,f2;
	int ma,na;

	if(F<0){
		F=0;
	}

	x=n/(n+m*F);

	ma=m%2;
	na=n%2;

	if((na==1)&&(ma==1)){
		Ix=1-2/3.141592*atan(sqrt(1/x-1));
		U=sqrt(x*(1-x))/3.141592;
	}
	else if((na==1)&&(ma==0)){
		Ix=sqrt(x);
		U=sqrt(x)*(1-x)/2;
	}
	else if((na==0)&&(ma==1)){
		Ix=1-sqrt(1-x);
		U=x*sqrt(1-x)/2;
	}
	else if((na==0)&&(ma==0)){
		Ix=x;
		U=x*(1-x);
	}

	if(ma==0){
		ma=2;
	}

	if(na==0){
		na=2;
	}

	for(f2=na;f2<n;f2+=2){
		f1=ma;
		Ix=Ix-2*U/f2;
		U=(f1+f2)*x*U/f2;
	}
	for(f1=ma;f1<m;f1+=2){
		Ix=Ix+2*U/f1;
		U=(f1+f2)*(1-x)*U/f1;
	}

	return Ix;
}

double CalcTValue(int f,double Q)
{
	//Paulson-Takeuchi偺嬤帡幃傪巊梡丅
	//妋棪Q偵懳墳偡傞帺桼搙f偺倲暘晍昞偺抣傪寁嶼偡傞丅
	//弌揟丗僷僜僐儞摑寁夝愅僴儞僪僽僢僋(婎慴摑寁曇) p.289丅

	double u,Y1,Y2,Y3,Y4,Y5;

	if(f>0){
		u=CalcZ(Q);
		Y1=(pow(u,3)+u)/4;
		Y2=(5*pow(u,5)+16*pow(u,3)+3*u)/96;
		Y3=(3*pow(u,7)+19*pow(u,5)+17*pow(u,3)-15*u)/384;
		Y4=(79*pow(u,9)+776*pow(u,7)+1482*pow(u,5)-1920*pow(u,3)-945*u)/92160L;
		Y5=(27*pow(u,11)+339*pow(u,9)+930*pow(u,7)-1782*pow(u,5)-765*pow(u,3)+17955*u)/368640L;
		return u+Y1/f+Y2/f/f+Y3/f/f/f+Y4/f/f/f/f+Y5/f/f/f/f/f;
	}
	else{
		return NO_DATA;
	}
}

double CalcFValue(int m,int n,double Q)
{
	//妋棪Q偵懳墳偡傞帺桼搙m,n偺俥暘晍昞偺抣傪寁嶼偡傞丅
	//弌揟丗僷僜僐儞摑寁夝愅僴儞僪僽僢僋(婎慴摑寁曇) p.288-289丅

	double u,temp,a,b,free1,free2,xL,x0,xU,alpha1,alpha2;
	int i;

	if((Q<0)||(Q>1)){
		return NO_DATA;
	}
	else if((m>30)&&(n>30)){
		u=CalcZ(Q);
		free1=m;
		free2=n;
		a=2/(9*free1);
		b=2/(9*free2);
		temp=((1-a)*(1-b)+u*sqrt((1-a)*(1-a)*b+(1-b)*(1-b)*a-a*b*u*u))/((1-b)*(1-b)-b*u*u);
		return temp*temp*temp;
	}
	else{
		xL=0;
		x0=0.5;
		xU=1;
		alpha1=Q;
		temp=(1/x0-1)*m/n;
		alpha2=CalcPValue_F(m,n,temp);
		for(i=0;i<500;i++){
			if(alpha1>alpha2){
				xL=x0;
			}
			else{
				xU=x0;
			}
			x0=(xL+xU)/2;
			temp=(1/x0-1)*m/n;
			alpha2=CalcPValue_F(m,n,temp);
		}

		return (1/x0-1)*m/n;
	}

}

double CalcFisher(int a,int b,int c,int d)
{
	//Fisher専掕偺妋棪傪寁嶼偡傞丅
	//悢抣偑僆乕僶乕僼儘乕偟側偄傛偆偵懳悢傪偲傞丅
	//孮攏戝妛偺惵栘巵偺傾儖僑儕僘儉傪巊梡偟偰偄傞丅
	//
	//	a	|	b
	//--------------
	//	c	|	d
	//
	
	int n,i;
	double logdata[MAX_FISHER],temp;

	n=a+b+c+d;
	logdata[0]=0;

	for(i=1;i<MAX_FISHER;i++){
		logdata[i]=logdata[i-1]+log(i);
	}

	temp=logdata[a+b]+logdata[a+c]+logdata[b+d]+logdata[c+d]
			-logdata[a]-logdata[b]-logdata[c]-logdata[d]-logdata[n];

	return exp(temp);
}

double CalcPValue_WSR(int count,int T)
{
	//Wilcoxon signed rank専掕偺椉懁妋棪傪寁嶼偡傞丅
	//孮攏戝妛偺惵栘巵偺傾儖僑儕僘儉傪巊梡偟偰偄傞丅

	int i,n,sum,limit;
	double* table,denom,total=0;
	double p=0;
	
	limit=(count+1)*count/2+1;
	table=new double[limit];
	
	for(i=0;i<limit;i++){
		table[i]=0;
	}
	
	table[0]=1;
	
	for(n=1,sum=0;n<=count;sum+=n++){
		for(i=sum;i>=0;i--){
			table[i+n]+=table[i];
		}
	}
	
	denom=pow(2.0,-count);
	
	if(T<0){
		T=0;
	}
	else if(T>=limit){
		T=limit-1;
	}
	
	for(i=0;i<limit;i++){
		total+=(double)table[i]*denom;
		if(T==i){
			p=total*2;
		}
	}

	delete table;
	return p;
}

double* g_table;

double CalcPValue_MW(int n1,int n2,int U)
{
	//Mann-Whitney偺U専掕偺椉懁妋棪傪寁嶼偡傞丅
	//U偵偼摑寁検U1丄U2偺偆偪彫偝偄曽傪戙擖偡傞丅
	//孮攏戝妛偺惵栘巵偺傾儖僑儕僘儉傪巊梡偟偰偄傞丅

	int i;
	double total=0.0;
	double denom;
	double p=0;
	
	g_table=new double[n1*n2+1];

	for(i=0;i<n1*n2+1;i++){
		g_table[i]=0;
	}

	u(n1,n2,0,n1*n2);
	denom=1.0/combination(n1+n2,n2);
	for(i=0;i<=n1*n2;i++){
		total+=g_table[i];
		if(U==i){
			p=total*denom*2;
		}
	}
	
	delete g_table;
	return p;
}

void increment(double* q,int n)
{
	//孮攏戝妛偺惵栘巵偺傾儖僑儕僘儉傪巊梡偟偰偄傞丅

	while(n-->=0){
		(*q++)++;
	}
}

void u(int n,int m,int begin,int end)
{
	//孮攏戝妛偺惵栘巵偺傾儖僑儕僘儉傪巊梡偟偰偄傞丅

	if((n==1)||(m==1)){
		increment(g_table+begin,end-begin);
	}
	else{
		u(n-1,m,begin,begin+(n-1)*m);
		u(n,m-1,end-n*(m-1),end);
	}
}

double combination(int n,int i)
{
	//孮攏戝妛偺惵栘巵偺傾儖僑儕僘儉傪巊梡偟偰偄傞丅

	int j;
	double retv;
	
	if((i<0)||(i>n)||(n<0)){
		return 0;
	}
	
	retv=1.0;
	
	if(i>n-i){
		i=n-i;
	}
	
	for(j=0;j<i;j++){
		retv*=(double)(n-j)/(double)(j+1);
	}
	
	return retv;
}