14stdt0.c

特殊函数T分布函数的计算,有兴趣的可以试试

 
  #include "math.h"
  #include "stdio.h"
  double gam1(x)
  double x;
  { int i;
    double y,t,s,u;
    static double a[11]={ 0.0000677106,-0.0003442342,
           0.0015397681,-0.0024467480,0.0109736958,
           -0.0002109075,0.0742379071,0.0815782188,
           0.4118402518,0.4227843370,1.0};
    if (x<=0.0)
      { printf("err**x<=0!\n"); return(-1.0);}
    y=x;
    if (y<=1.0)
      { t=1.0/(y*(y+1.0)); y=y+2.0;}
    else if (y<=2.0)
      { t=1.0/y; y=y+1.0;}
    else if (y<=3.0) t=1.0;
    else
      { t=1.0;
        while (y>3.0)
          { y=y-1.0; t=t*y;}
      }
    s=a[0]; u=y-2.0;
    for (i=1; i<=10; i++)
      s=s*u+a[i];
    s=s*t;
    return(s);
  }


  double beta(a,b,x)
  double a,b,x;
  { double y;
    double bt(double,double,double);
    if (a<=0.0)
      { printf("err**a<=0!"); return(-1.0);}
    if (b<=0.0)
      { printf("err**b<=0!"); return(-1.0);}
    if ((x<0.0)||(x>1.0))
      { printf("err**x<0 or x>1 !");
        return(1.0e+70);
      }
    if ((x==0.0)||(x==1.0)) y=0.0;
    else
      { y=a*log(x)+b*log(1.0-x);
        y=exp(y);
        y=y*gam1(a+b)/(gam1(a)*gam1(b));
      }
    if (x<(a+1.0)/(a+b+2.0))
      y=y*bt(a,b,x)/a;
    else
      y=1.0-y*bt(b,a,1.0-x)/b;
    return(y);
  }

  static double bt(a,b,x)
  double a,b,x;
  { int k;
    double d,p0,q0,p1,q1,s0,s1;
    p0=0.0; q0=1.0; p1=1.0; q1=1.0;
    for (k=1; k<=100; k++)
      { d=(a+k)*(a+b+k)*x;
        d=-d/((a+k+k)*(a+k+k+1.0));
        p0=p1+d*p0; q0=q1+d*q0; s0=p0/q0;
        d=k*(b-k)*x;
        d=d/((a+k+k-1.0)*(a+k+k));
        p1=p0+d*p1; q1=q0+d*q1; s1=p1/q1;
        if (fabs(s1-s0)<fabs(s1)*1.0e-07)
          return(s1);
      }
    printf("a or b too big !");
    return(s1);
  }


  double stdt(t,n)
  int n;
  double t;
  { double y;
    if (t<0.0) t=-t;
    y=1.0-beta(n/2.0,0.5,n/(n+t*t));
    return(y);
  }


  main()
  { int n;
    double t,y;
    printf("\n");
    for (n=1; n<=5; n++)
      { t=0.5; y=stdt(t,n);
        printf("P(%4.2f, %d)=%e\n",t,n,y);
        t=5.0; y=stdt(t,n);
        printf("P(%4.2f, %d)=%e\n",t,n,y);
      }
    printf("\n");
  }