d13r6b.cpp

使用VC++编写的大量数学算法的源代码

#include "iostream.h"
#include "math.h"
#include "stdlib.h"

double gammln(double xx)
{
	double cof[7];
	cof[1] = 76.18009173;
	cof[2] = -86.50532033;
	cof[3] = 24.01409822;
	cof[4] = -1.231739516;
	cof[5] = 0.00120858003;
	cof[6] = -0.00000536382;
	double stp = 2.50662827465;
	double half = 0.5;
	double one = 1.0;
	double x,tmp,ser,fpf = 5.5;
	x = xx - one;
	tmp = x + fpf;
	tmp = (x + half) * log(tmp) - tmp;
	ser = one;
	for (int j = 1; j<=6; j++)
	{
		x = x + one;
		ser = ser + cof[j] / x;
	}
	return tmp + log(stp * ser);
}

double betacf(double a, double b, double x)
{
	int em,tem,m,itmax = 100;
	double d,ap,bp,aap,bpp,aold,eps = 0.0000003;
	double am = 1.0;
	double bm = 1.0;
	double az = 1.0;
	double qab = a + b;
	double qap = a + 1.0;
	double qam = a - 1.0;
	double bz = 1.0 - qab * x / qap;
	for (m = 1; m<=itmax; m++)
	{
		em = m;
		tem = em + em;
		d = em * (b - m) * x / ((qam + tem) * (a + tem));
		ap = az + d * am;
		bp = bz + d * bm;
		d = -(a + em) * (qab + em) * x / ((a + tem) * (qap + tem));
		aap = ap + d * az;
		bpp = bp + d * bz;
		aold = az;
		am = ap / bpp;
		bm = bp / bpp;
		az = aap / bpp;
		bz = 1.0;
		if (fabs(az - aold) < eps * fabs(az))
		{
			return az;
		}
	}
	cout<< "a or b too big, or itmax too small"<<endl;
	exit(1);
}

void avevar(double data[], int n, double& ave, double& var)
{
    ave = 0.0;
    var = 0.0;
	int j;
	double s;
    for (j = 1; j<=n; j++)
	{
        ave = ave + data[j];
    }
    ave = ave / n;
    for (j = 1; j<=n; j++)
	{
        s = data[j] - ave;
        var = var + s * s;
    }
    var = var / (n - 1);
}

double betai(double a, double b, double x)
{
	double bt,aaa;
	if (x < 0.0 || x > 1.0)
	{
		cout<<"bad argument x in betai"<<endl;
	}
	if (x == 0.0 || x == 1.0)
	{
		bt = 0.0;
	}
	else
	{
		aaa = gammln(a + b) - gammln(a) - gammln(b);
		bt = exp(aaa + a * log(x) + b * log(1.0 - x));
	}
	if (x < (a + 1.0) / (a + b + 2.0))
	{
		return bt * betacf(a, b, x) / a;
	}
	else
	{
		return 1.0 - bt * betacf(b, a, 1.0 - x) / b;
	}
}

void tptest(double data1[], double data2[], int n, double& t, double& prob)
{
	double ave1,ave2,var1,var2;
    avevar(data1, n, ave1, var1);
    avevar(data2, n, ave2, var2);
    double cov = 0.0;
    for (int j = 1; j<=n; j++)
	{
        cov = cov + (data1[j] - ave1) * (data2[j] - ave2);
    }
    int df = n - 1;
    cov = cov / df;
    double sd = sqrt((var1 + var2 - 2.0 * cov) / n);
    t = (ave1 - ave2) / sd;
    prob = betai(0.5 * df, 0.5, df / (df + t * t));
}

double ran1(long& idum)
{
    int j,iff=-1;
	static long ix1,ix2,ix3;
	static double r[98];
    long m1 = 259200; long m2 = 134456; long m3 = 243000;
	long ia1 = 7141; long ia2 = 8121; long ia3 = 4561;
	long ic1 = 54773; long ic2 = 28411; long ic3 = 51349;
	double rm1 = 0.0000038580247; double rm2 = 0.0000074373773;    
    if (idum < 0 || iff == 0)
	{
        iff = 1;
        ix1 = (ic1 - idum) % m1;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = ix1 % m2;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix3 = ix1 % m3;
        for (j = 1; j<=97; j++)
		{
            ix1 = (ia1 * ix1 + ic1) % m1;
            ix2 = (ia2 * ix2 + ic2) % m2;
            r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
        }
        idum = 1;
    }
    ix1 = (ia1 * ix1 + ic1) % m1;
    ix2 = (ia2 * ix2 + ic2) % m2;
    ix3 = (ia3 * ix3 + ic3) % m3;
    j = 1 + int((97 * ix3) / m3);
    if (j > 97 || j < 1)
	{
		cout<<"abnormal exit in ran1"<<endl;
		exit(1);
	}
	double temp=r[j];
    r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
	return temp;
}

double gasdev(long& idum)
{
    static int iset;
	static double gset;
	double v1,v2,r,fac;
    if (iset == 0)
	{
		do
		{
			v1 = 2.0 * ran1(idum) - 1.0;
			v2 = 2.0 * ran1(idum) - 1.0;
			r = v1 * v1 + v2 * v2;
		}while (r >= 1.0 || r == 0);
		fac = sqrt(-2.0 * log(r) / r);
		gset = v1 * fac;
		iset = 1;
		return v2 * fac;
	}
    else
	{
		iset = 0;
		return gset;
    }
}

void main()
{
    //program d13r6b
    //driver for routine tptest
    //compare two correlated distributions vs. two
    //uncorrelated distributions
    int i,j,npts = 500;
    double eps = 0.01;
    int nshft = 10;
    double anoise = 0.3;
    double data1[501], data2[501], data3[501];
    long idum = -5;
    cout<<endl;
    cout<<"                 Correlated:        Uncorrelated:"<<endl;
    cout<<"     shift       t    probability    t    probability"<<endl;
    double offset = (nshft / 2) * eps;
    for (j = 1; j<=npts; j++)
	{
        double gauss = gasdev(idum);
        data1[j] = gauss;
        data2[j] = gauss + anoise * gasdev(idum);
        data3[j] = gasdev(idum) + anoise * gasdev(idum);
    }
	double ave1,ave2,ave3,var1,var2,var3;
    avevar(data1, npts, ave1, var1);
    avevar(data2, npts, ave2, var2);
    avevar(data3, npts, ave3, var3);
    for (j = 1; j<=npts; j++)
	{
        data1[j] = data1[j] - ave1 + offset;
        data2[j] = data2[j] - ave2;
        data3[j] = data3[j] - ave3;
    }
	double t1,t2,prob1,prob2;
	cout.setf(ios::fixed|ios::right);
	cout.precision(4);
	for (i = 1; i<=nshft + 1; i++)
	{
        double shift = i * eps;
        for (j = 1; j<=npts; j++)
		{
            data2[j] = data2[j] + eps;
            data3[j] = data3[j] + eps;
        }
        tptest(data1, data2, npts, t1, prob1);
        tptest(data1, data3, npts, t2, prob2);
		cout.width(10);
        cout<<shift;
		cout.width(10);
        cout<<t1;
		cout.width(10);
        cout<<prob1;
		cout.width(10);
        cout<<t2;
		cout.width(10);
        cout<<prob2<<endl;
	}
}