分布函数分位点及p-value计算及程序实现

各分布函数分位点及p-value计算及程序实现

一、各分布函数计算

1.标准正态分布函数的计算

方法（1）：利用泰勒展示求解

由于分布函数

2

2

()

2

x

t

F t e dx

π

-

-∞

=?不能被求解成为简单的初等函数，所以不能

直接用定义初等复合函数的方法计算其值。通过把正态密度函数

2

2

f()

2

x

x e

π

-

=

展开为泰勒级数,对泰勒级数的各项求定积分，然后相加各项数值即得分布函数F（t）的值。

22462

2

123

()=...+(-1)...

1!2!3!!

22

x n

n

n

x x x x

f x e

n

ππ

-

=-+-++

?2?2?2?2

（1）（1）

22

22

()0.5

22

x x

t t

F t e dx e dx

ππ

--

-∞

==+

??（2）

由（1），（2）得出

35721

123

()0.5(...+(-1)...)

1!2!3!!n

2

n

n

n

t t t t

F t t

n

π

+

=+-+-++

?2?3?2?5?2?7?2?(2+1)

（3）因为（3）式是一个无穷级数，计算时只能对前面有限项相加。选择越多，结果越准确。但数值精度并不要求无限高，可选取50项就足够，因而这里n=0,1,2…50。

流程图如下：

代开C++软件，输入以下程序代码：

#include

#include

#define PI 3.141592653

main( )

{

float z,t,y=0,jsx,F,b=1 ; //jsx为每个级数项的值，F为所要求分布函数的值

int n=1,a=1,m=1 ;

printf("Please enter t:");

scanf("%f",&t); //输入分位点t

while(n<101) //n<101,得出的级数项值共有50个

{

z=pow(t,n);

jsx=a*z/(b*n);

y+=jsx;

a=-1*a ; //正负符号标志 b=b*2*m ; m++; n=n+2; }

F=0.5+y/sqrt(2*PI);

printf("The result you want is:\n"); printf("%f",F); }

输入数据：比如t=0.02,x=-0.01,x=0 输出:

方法（2）：积分的近似算法： 正态分布函数2

2

()2x t

F t dx π

-

=?

的计算归为积分计算。

设将区间[a,b]分为n 等份，共有n+1个分点。分点x ,,0,1,...k b a

a kh h k n

-=+=

=，在每个子区间1[x ,x ]k k +（k=0,1,….,n-1）上采用梯形公式()[()()]2

b a

b a f x dx f a f b -≈+?则可得复合梯形公式11

101

[()()][()2()()]22n n n k k k k k h h

T f x f x f a f x f b --+===+=++∑∑。如果将求积区间再二分

一次，则可提高求积精度。此时分点增至2n+1个，则由原来每个子区间1[x ,x ]

k k +经过二分只增加了一个分点112

1

x (x x )2k k k ++=+，再用复合梯形公式求得该子区

间上的积分值为112

[()2()()]4k k k h

f x f x f x ++++（其中b a h n -=代表二分前的步

长）。将每个子区间上的积分值相加得：11

21100

2[()()]()42n n n k k k k k h h T f x f x f x --++===++∑∑

C++程序算法： #include #include double f(double x) {

return exp(-x*x/2);

}

double F(double a,double b,double ep=1e-6)

{

double h,s1=0,s2=(b-a)*(f(a)+f(b))/2;

int n,k;

for(int n=1;fabs(s1-s2)>ep;n*=2)

{

h=(b-a)/n;

s1 = s2;

s2 = 0;

for(int k=0;k

{

s2 += h*f(a+(k+0.5)*h);

}

s2 = (s1+s2)/2;

}

return s2*sqrt(1/(8*atan(1.0)));

}

int main()

{

double c=0,x;

printf("please enter x:");

scanf("%lf",&x);

printf("%lf",0.5+F(c,x));

return 0;

}

结果显示：（输入x的值为2）

输出：

方法（3）：连分式展开（算法程序见下面最后部分的总程序编写）

Φ（x）为标准正态分布函数，其是对称的，只要求出x>0时Φ（x）的值也就可以求得x<0时Φ（x）=1-Φ（-x）。由于Φ（x）的近似式：

Φ（x）≈

222

1()2(1)

......,03 2135(21)

()12

1......,3

n

f x x x x nx

x

n

f x n

x

x x x x

?--

+++++≤≤

?+

?

?

?

-++++>

??

当时

当时

其中

2

2

()

x

f x-

=采用递推式表达，

当0≤x ≤3时，令12

1a 0

(1),,1,......1(21)n k k k kx

a k n n k a

++=??-?==-?++?

则Φ（x ）= 11()21f x x a ++ 当x>3时，令11a 0

,,1,.....1n k k k

a k n n x a

++=??

?==-?+?

则Φ（x ）= 1()1f x x a -+ 一般当n>28时，精度可达10-12

2.Beta 分布函数的计算

Beta 分布的分布函数I (a,b) x 递推如下：

1I (a+1,b)=I (a,b)-(,)1I (a,b+1)=I (a,b)+(,)

(1,)(,)

(,1)(1)(,)

x

x x x x x

x x

x x U a b a

U a b b a b U a b xU a b a a b U a b x U a b a ?

??

???+?+=??+?+=-?

其中1(,)(1)(,)a b x

U a b x x B a b =- 由于利用Beta 分布的分布函数计算t 分布，F 分布，二项分布时，参数a,b 的值要么是正整数，要么就是1/2的倍数。所以考虑参数a,b 的值是正整数或者是1/2倍数情况下I (a,b) x 的计算。此时递推公式初值选取有以下4种：

1111112,(,)(,)=1-222222x x a b U π===当时，

111,1(,1)),I (222

x x a b U x ===-当时，

1111

1,(1,)I (1,2222

x x a b U ===当时，

1,1(1,1)(1),I (1,1)=x x x a b U x x ===-当时，

3.t 分布函数的计算 I (a,b) x 与 t 分布的关系

2

1

11(,),0

222(|),11,0(,)222

x x

n I t n T t n x n t n t I ?-?>?==?≤+???其中。由于n 为整数，所以n/2可以递归到?，即可以使用Beta 分布的递推公式。

4.F 分布函数计算

I (a,b) x 与 F 分布的关系

(|,)(,),22y m n mx

F x m n I y n mx

==+其中。由于m,n 均为整数，所以m/2,n/2均可以

递归到?，?，即可以使用Beta 分布的递推公式。

5．二项分布函数的计算

由二项分布函数的定义：()[]00,0x |,)(1)

,01,x n k n k

k k x B n p p p x n x n

-=

6．卡方分布函数的计算

卡方分布函数H （x|n ）的递推计算公式如下： (|)(|2)2(|),(3,4,.....)(|)(|2)2H x n H x n f x n n x f x n f x n n =--??

=?=-?-?其中1221(|)()22()2

n x

x f x n e n --=Γ

则递推初值为2

22(|2)1,(|1)1(|1)21,(|2)2

x x x H x e f x H x f x e ---?=-=

????

=Φ-=??

7．Poisson 分布函数的计算

[]P x |1H 2|2x 1λλ=-+（）（（））,x ≥0 二、各分布函数分位数的计算

1.标准正态分布分位数的计算 采用精度达到 1.2x10

-8

的Toda 近似公式：110

2

()

i i u y b y α≈∑,其中

[]l n 4(1)y α

α=-- b i 的各数值用数组B[11]表示如下：

B[11]=(1.570796288,0.0370*******,-0.8364353589e-3,-0.2250947176e-3, 0.6841218299e-5,0.5824238515e-5,-0.1045274970e-5,0.8360937017e-7, -0.3231081277e-8,0.3657763036e-10,0.6936233982e-12)

2. t 分布函数的计算

由T 分布函数与I (a,b) x 的关系可得：

2

1110t 0,t t | n (,)1,x 2222t x n n

T I n αααααα<<>=-=

+当时，且满足（）=1-其中由此得1(,)22x n I

，则122

1(,),t ()22t p

n n

n n ααβ-==+1-1

t 0t ()=-t ()2

n n αααα><当时,，

3. F 分布函数分位数的计算

F 分布的α分位数(m,)F n α满足(|m,)1F F n αα=-。由F 分布函数与Beta 分布函

数的关系得F α满足(|m ,)

(,)1,y 22y mF m n

F F n I n mF ααα

α==-=+其中，则

(,)22

mF m n

n mF αααβ=+，

所以(,)

22(m,)(1(,))

22

m n n F n m n

m αααββ=-

4. 卡方分布函数α分位数2()n αχ的计算

①当n=1时，由H （x|1）=2Φ

-1=1-α,得Φ

=1-α/2，

u α2

=，22(1)(u )ααχ2

=

②当n=2时，2

(2)χ分布就是12λ=的指数分布1

()2

e ,由2(x |2)11x

H e α-=-=-得

2(2)2ln αχα=-

③当n ≥3

时，23

2()[19n n u n ααχ2

≈-

+，

其中u α是标准正态分布的α分位数。 程序算法如下：

#include #include #include #include

#define Pi 3.141592653

#define MaxTime 500 //定义最大迭代步数 #define eps 1.0e-10 //定义精度

double KaFangFx(float x,int Freedom); double PoissonFx(float x,double p); double KaFangFx1(float x,int Freedom); double GaossFx1(float x);

double KaFangUa0(float af,int Freedom); double KaFangUa(float af,int Freedom); double Gama(int n);

double KaFangPx(float x,int Freedom); double GaossFx(float x); double GaossUa(float af); double GaossUa1(float af); double GaossPx(float x);

double BetaFx(float x,double a,double b); double TdistFx(float x,int Freedom);

double FdistFx(float x,int Freedom_m,int Freedom_n); double BetaUa(float af,double a,double b); double TdistUa(float af,int Freedom);

double FdistUa(float af,int Freedom_m,int Freedom_n); double BinominalFx(float x,float p,int n);

double BetaFx(float x,double a,double b)//贝塔分布函数{

int m,n;

double I,U;

double ta,tb;

m=(int)(2*a);

n=(int)(2*b);

if(m%2==1&&n%2==1)

{

ta=0.5;

tb=0.5;

U=sqrt(x*(1.0-x))/Pi;

I=1.0-2.0/Pi*atan(sqrt((1.0-x)/x));

}

else if(m%2==1&&n%2==0)

{

ta=0.5;

tb=1;

U=0.5*sqrt(x)*(1.0-x);

I=sqrt(x);

}

else if(m%2==0&&n%2==1)

{

ta=1;

tb=0.5;

U=0.5*x*sqrt(1.0-x);

I=1.0-sqrt(1.0-x);

}

else if(m%2==0&&n%2==0)

{

ta=1;

tb=1;

U=x*(1.0-x);

I=x;

}

while(ta

{

I=I-U/ta;

U=(ta+tb)/ta*x*U;

ta+=1;

}

while(tb

{

I=I+U/tb;

U=(ta+tb)/tb*(1.0-x)*U;

tb+=1;

}

return I;

}

double TdistFx(float x,int Freedom_t) //T分布t(n)在x处的分布函数计算

{

double t,prob;

printf("输入t分布的自由度Freedom_t:\n");

scanf("%d",&Freedom_t);

printf("输入t分布的x:\n");

scanf("%f",&x);

t=Freedom_t/(Freedom_t+x*x);

if(x>0)

prob=1.0-0.5*BetaFx(t,Freedom_t/2.0,0.5);

else

prob=0.5*BetaFx(t,Freedom_t/2.0,0.5);

printf("T分布T（%f,%d）=%f\n",x,Freedom_t,prob);

}

double FdistFx(float x,int Freedom_m,int Freedom_n) //F分布F（m,n）在x处分布函数计算

{

double y,prob;

printf("F分布的自由度Freedom_m,Freedom_n:\n");

scanf("%d%d",&Freedom_m,&Freedom_n);

printf("输入P(X<=x)分布值P中的x:\n");

scanf("%f",&x);

if(x<=0) return 0.0;

y=Freedom_m*x/(Freedom_n+Freedom_m*x);

prob=BetaFx(y,Freedom_m/2.0,Freedom_n/2.0);

printf("F分布F（%f,%d,%d)=%f\n",x,Freedom_m,Freedom_n,prob);

}

double BinominalFx(float x,float p,int n)//二项分布的分布函数计算{double prob;

printf("输入二项分布参数p:\n");

scanf("%f",&p);

printf("输入二项分布x:\n");

scanf("%f",&x);

printf("输入二项分布n:\n");

scanf("%d",&n);

if(x<0) prob=0.0;

else if(x>=n) prob=1.0;

else prob=BetaFx(1.0-p,n-int(x),int(x)+1);

printf("二项分布B(%f，%f，%d)=%f\n",x,p,n,prob);

}

double BetaUa(float af,double a,double b)//贝塔分布的上侧分位数计算{

int times=0;

double x1,x2,xn;

double f1,f2,fn,ua;

x1=0.0;

x2=1.0;

f1=BetaFx(x1,a,b)-(1.0-af);

f2=BetaFx(x2,a,b)-(1.0-af);

while(fabs((x2-x1)/2.0)>eps)

{

xn=(x1+x2)/2.0;

fn=BetaFx(xn,a,b)-(1.0-af);

if(f1*fn<0) x2=xn;

else if(fn*f2<0) x1=xn;

f1=BetaFx(x1,a,b)-(1.0-af);

f2=BetaFx(x2,a,b)-(1.0-af);

times++;

if(times>MaxTime) break;

}

printf("times=%5d\n",times);

ua=xn;

return ua;

}

double TdistUa(float af,int Freedom)//T分布的上侧分位数计算

{

double ua,tbp,bf;

printf("T分布的上侧分位概率af:\n");

scanf("%f",&af);

printf("T分布的自由度Freedom\n");

scanf("%d",&Freedom);

bf=1-af;

printf("F分布的下侧分位数bf=%f\n",bf);

if(af<=0.5)

{

tbp=BetaUa(1-2*af,Freedom/2.0,0.5);

ua=sqrt(Freedom/tbp-Freedom);

}

else if(af>0.5)

{

tbp=BetaUa(1-2*(1.0-af),Freedom/2.0,0.5);

ua=-sqrt(Freedom/tbp-Freedom);

}

printf("T分布下侧分位数T(%f,%d)=%f\n",bf,Freedom,ua);

}

double FdistUa(float af,int Freedom_m,int Freedom_n)//F分布的上侧分位数计算

{

float ua,tbp,bf;

printf("F分布的上侧概率分位数af:\n");

scanf("%f",&af);

printf("F分布的自由度Freedom_m,Freedom_n:\n");

scanf("%d%d",&Freedom_m,&Freedom_n);

bf=1-af;

printf("F分布的下侧分位数bf=%f\n",bf);

tbp=BetaUa(af,Freedom_m/2.0,Freedom_n/2.0);

ua=Freedom_n*tbp/(Freedom_m*(1.0-tbp));

printf("F分布F(%f,%d,%d)=%f\n",bf,Freedom_m,Freedom_n,ua);

}

double KaFangFx1(float x,int Freedom)

{

int k,n;

double f,h,prob;

k=Freedom%2;

if(k==1)

{

f=exp(-x/2.0)/sqrt(2*Pi*x);

h=2.0*GaossFx1(sqrt(x))-1.0;

n=1;

while(n

{

n=n+2;

f=x/(n-2.0)*f;

h=h-2.0*f;

}

}

else

{

f=exp(-x/2.0)/2.0;

h=1.0-exp(-x/2.0);

n=2;

while(n

{

n=n+2;

f=x/(n-2.0)*f;

h=h-2.0*f;

}

}

prob=h;

return prob;

}

double KaFangFx(float x,int Freedom) //卡方分布函数计算{

int k,n;

double f,h,prob;

printf("卡方分布的自由度Freedom,\n");

scanf("%d",&Freedom);

printf("输入P(X<=x)分布值P中的x:\n");

scanf("%f",&x);

k=Freedom%2;

if(k==1)

{

f=exp(-x/2.0)/sqrt(2*Pi*x);

h=2.0*GaossFx1(sqrt(x))-1.0;

n=1;

while(n

{

n=n+2;

f=x/(n-2.0)*f;

h=h-2.0*f;

}

}

else

{

f=exp(-x/2.0)/2.0;

h=1.0-exp(-x/2.0);

n=2;

while(n

{

n=n+2;

f=x/(n-2.0)*f;

h=h-2.0*f;

}

}

prob=h;

printf("卡方分布K(%f,%d)=%f\n",x,Freedom,prob);

}

double PoissonFx(float x,float p) //Poisson分布函数的计算

{

double prob;

printf("输入Poisson分布参数p:\n");

scanf("%f",&p);

printf("输入Poisson分布x:\n");

scanf("%f",&x);

prob=1.0-KaFangFx1(2*p,2*(int(x))+1);

printf("Poisson分布P（%f,%f）=%f\n",x,p,prob);

}

double KaFangUa0(float af,int Freedom)

{

double ua,p,temp;

if(Freedom==1)

{

p=1.0-(1.0-af+1.0)/2.0;

temp=GaossUa(p);

ua=temp*temp;

}

else if(Freedom==2)

{

ua=-2.0*log(af);

}

else

{

temp=1.0-2.0/(9.0*Freedom)+sqrt(2.0/(9.0*Freedom))*GaossUa(af); ua=Freedom*(temp*temp*temp);

}

return ua;

}

double KaFangUa(float af,int Freedom)//卡方分布的上侧分位数的计算{

int times;

double ua,x0,xn,bf;

printf("输入卡方分布的自由度Freedom:\n");

scanf("%d",&Freedom);

printf("输入卡方分布af:\n");

scanf("%f",&af);

bf=1-af;

printf("F分布的下侧分位数bf=%f\n",bf);

x0=KaFangUa0(af,Freedom);

printf("x0=%12.8f\n",x0);

if(Freedom==1||Freedom==2)

{

ua=x0;

}

else

{//采用牛顿迭代法

times=1;

xn=x0-(KaFangFx1(x0,Freedom)-1+af)/KaFangPx(x0,Freedom);

while(fabs(xn-x0)>eps)

{

x0=xn;

xn=x0-(KaFangFx1(x0,Freedom)-1+af)/KaFangPx(x0,Freedom);

times++;

if(times>MaxTime) break;

}

printf("times=%5d\n",times);

ua=xn;

}

printf("卡方分布下侧分位数K(%f，%d)=%f\n",bf,Freedom,ua);

}

double KaFangPx(float x,int Freedom)//卡方分布的密度函数

{

double p,g;

if(x<=0) return 0.0;

g=Gama(Freedom);

p=1.0/pow(2.0,Freedom/2.0)/g*exp(-x/2.0)*pow(x,Freedom/2.0-1.0); return p;

}

double Gama(int n)//伽马分布函数Gama(n/2)

{

double g;

int i,k;

k=n/2;

if(n%2==1)

{

g=sqrt(Pi)*0.5;

for(i=1;i

g*=(i+0.5);

}

else

{

g=1.0;

for(i=1;i

g*=i;

}

return g;

}

double GaossFx1(float x)

{

double prob,t,temp;

int i,n,symbol;

temp=x;

if(x<0)

x=-x;

n=28;

if(x>=0&&x<=3.0)

{

t=0.0;

for(i=n;i>=1;i--)

{

if(i%2==1) symbol=-1;

else symbol=1;

t=symbol*i*x*x/(2.0*i+1.0+t); }

prob=0.5+GaossPx(x)*x/(1.0+t); }

else if(x>3.0)

{t=0.0;

for(i=n;i>=1;i--)

t=1.0*i/(x+t);

prob=1-GaossPx(x)/(x+t);

}

x=temp;

if(x<0)

prob=1.0-prob;

return prob;

}

double GaossFx(float x)//正态分布函数的计算{

double prob,t,temp;

int i,n,symbol;

printf("输入P(X<=x)分布值P中的x:\n"); scanf("%f",&x);

temp=x;

if(x<0)

x=-x;

n=28;//连分式展开的阶数

if(x>=0&&x<=3.0)//连分式展开方法

{

t=0.0;

for(i=n;i>=1;i--)

{

if(i%2==1) symbol=-1;

else symbol=1;

t=symbol*i*x*x/(2.0*i+1.0+t);

}

prob=0.5+GaossPx(x)*x/(1.0+t);

}

else if(x>3.0)

{t=0.0;

for(i=n;i>=1;i--)

t=1.0*i/(x+t);

prob=1-GaossPx(x)/(x+t);

}

x=temp;

if(x<0)

prob=1.0-prob;

printf("正态分布F(%f)=%f\n",x,prob);

}

double GaossUa(float af)

{double b[11]={

0.1570796288*10,

0.3706987906*1.0e-1,

-0.8364353589*1.0e-3,

-0.2250947176*1.0e-3,

0.6841218299*1.0e-5,

0.5824238515*1.0e-5,

-0.1045274970*1.0e-5,

0.8360937017*1.0e-7,

-0.3231081277*1.0e-8,

0.3657763036*1.0e-10,

0.6936233982*1.0e-12

};

double y,ua;

int i;

y=-log(4*af*(1.0-af));

ua=0.0;

for(i=0;i<=10;i++)

ua+=b[i]*pow(y,i);

ua=sqrt(ua*y);

if(af>0.5) ua=-ua;

return ua;

}

double GaossUa1(float af)//正态分布的上侧分位数计算{double b[11]={

0.1570796288*10,

0.3706987906*1.0e-1,

-0.8364353589*1.0e-3,

-0.2250947176*1.0e-3,

0.6841218299*1.0e-5,

0.5824238515*1.0e-5,

-0.1045274970*1.0e-5,

0.8360937017*1.0e-7,

-0.3231081277*1.0e-8,

0.3657763036*1.0e-10,

0.6936233982*1.0e-12

};

double y,ua,bf;

int i;

printf("正态分布上侧分位数af:\n");

scanf("%f",&af);

bf=1-af;

y=-log(4*af*(1.0-af));

ua=0.0;

for(i=0;i<=10;i++)

ua+=b[i]*pow(y,i);

ua=sqrt(ua*y);

if(af>0.5) ua=-ua;

printf("正态分布下侧分位数bf=%f,u=%f\n",bf,ua);

}

double GaossPx(float x)//正态分布的密度函数

{

double f;

f=1.0/sqrt(2.0*Pi)*exp(-x*x/2.0);

return f;

}

void main()

{int ch,Freedom_m,Freedom_n,Freedom,n;

float x, p,af;

while(1)

{

printf("\n §※§※§各分布函数值的计算§※§※§\t\n\n"); printf("\t 1. F分布F(m,n)在x处分布函数值\t\n");

printf("\t 2. T分布t(n)在x处分布函数值\t\n");

printf("\t 3. 二项分布的分布函数\t\n");

printf("\t 4. 卡方分布分布函数\t\n");

printf("\t 5. Poisson分布的分布函数\t\n");

printf("\t 6. 正态分布分布函数值\t\n");

printf("\n §※§※§各上侧分位数的计算§※§※§\t\n\n"); printf("\t 7. F分布F(m,n)的上侧分位数\t\n");

printf("\t 8. T分布t(n)的上侧分位数\t\n");

printf("\t 9. 卡方分布的上侧分位数\t\n");

printf("\t 10.正态分布的上侧分位数\t\n");

printf("\t 11. 退出系统.\t\n");

while(1)

{

printf(" 请选择: ");

scanf("%d",&ch);

if(ch>=1&&ch<=11)break;

else printf("\n 输入有误，请重新选择： 1~11: ");

}

switch(ch)

{

case 1:FdistFx(x,Freedom_m,Freedom_n);break;

case 2:TdistFx(x,Freedom);break;

case 3:BinominalFx(x,p,n);break;

case 4:KaFangFx(x,Freedom);break;

case 5:PoissonFx(x,p);break;

case 6:GaossFx(x);break;

case 7:FdistUa(af,Freedom_m,Freedom_n);break;

case 8:TdistUa(af,Freedom);break;

case 9:KaFangUa(af,Freedom);break; case 10:GaossUa1(af);break;

case 11:exit(0);

default: break;

}

}

}

结果显示截图：

参考书籍：《概率论与数理统计教程》高等教育出版《数值分析》第五版清华大学出版

《统计计算》高惠旋编著北京大学出版