rlp_math.cpp

Linux/windows 环境下跨平台开发程序

	else {
		gcf(&gammcf, a, x, &gln);			return gammcf;
		}
	return 0.0;
}

//continued fraction for incomplete beta function, used by betai()
double betacf(double a, double b, double x)
{
	double qap, qam, qab, em, tem, d, bz, bm = 1.0, bp, bpp, az = 1.0, am = 1.0, ap, app, aold;
	int m;

	qab = a+b;		qap = a+1.0;		qam = a-1.0;	bz = 1.0-qab*x/qap;
	for(m = 1; m <= 100; m++) {
		em = (double)m;			tem = em+em;
		d = em*(b-em)*x/((qam+tem)*(a+tem));
		ap = az + d * am;		bp = bz + d *bm;
		d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
		app = ap + d * az;		bpp = bp + d * bz;
		aold = az;				am = ap/bpp;
		bm = bp/bpp;			az = app/bpp;
		bz = 1.0;
		if(fabs(az-aold) <= (_PREC * fabs(az))) return az;	//success: return
		}
	return az;												//fail: iterations exceeded
}

//The incomplete beta function Ix(a,b) for 0 <= x <= 1
double betai(double a, double b, double x)
{
	double bt;

	if(x < 0.0 || x > 1.0) return 0.0;		//range !
	if(x == 0.0 || x == 1.0) bt = 0.0;
	else
		bt = exp(gammln(a+b)-gammln(a)-gammln(b)+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;
}

//the binomial coefficient
double bincof(double n, double k)
{
	if(n<0 || k<0 || k > n) return 0.0;
	return exp(gammln(n+1.0) - gammln(k+1.0) - gammln(n-k+1.0));
}

//the cumulative binomial distribution
double binomdistf(double k, double n, double p)
{
	if(k > n || n < 0.0 || p < 0.0 || p >1.0) return 0.0;
	return betai(n-k, k+1, p);
}

//the beta function
double betaf(double z, double w)
{
	return exp(gammln(z)+gammln(w)-gammln(z+w));
}

//the error function: not all compilers have a built in erf()
double errf(double x)
{
	return x < 0.0 ? -gammp(0.5, x*x) : gammp(0.5, x*x);
}

//the complementary error function
double  errfc(double x)
{
	double t, z, ans;

	z = fabs(x);
	t = 1.0/(1.0+0.5*z);
	ans = t * exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
		t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
		t*(-0.82215223+t*0.17087277)))))))));
	return x >= 0.0 ? ans : 2.0-ans;
}

//cumulative normal distribution
double norm_dist(double x, double m, double s)
{
	return 0.5 + errf((x - m)/(s * _SQRT2))/2.0;
}

//chi square distribution
double chi_dist(double x, double df, double)
{
	return gammq(df/2.0, x/2);
}

//t-distribution
double t_dist(double t, double df, double)
{
	return betai(df/2.0, 0.5, (df/(df+t*t)));
}

//poisson distribution
double pois_dist(double x, double m, double)
{
	return gammq(x+1.0, m);
}

//f-distribution
double f_dist(double f, double df1, double df2)
{
	return betai(df2/2.0, df1/2.0, df2/(df2+df1*f));
}

//---------------------------------------------------------------------------
// Inverse of statitistical functions:
//    Use a combination of the Newton-Raphson method and bisection
// Ref.: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Rcipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 273 ff.

// funcd supplies the function value fn and the derivative df of the function sf at x
void funcd(double x, double *fn, double *df, double (*sf)(double, double, double), 
		   double df1, double df2, double p)
{
	double y1, y2;

	*fn = (sf)(x, df1, df2);
	y1 = (sf)(x * 0.995, df1, df2);
	y2 = (sf)(x * 1.005, df1, df2);
	*df = (y2-y1)*100.0/x;
	*fn = *fn - p;
}

//distinv does actual bisection and Newton-Raphson root finding
double distinv(double (*sf)(double, double, double), double df1, double df2, double p, double x0)
{
	int j;
	double df, dx, dxold, f, fh, fl;
	double swap, temp, xh, xl, rts; 
	double x1 = 0.0001, x2 = 10000;
	char info[80];

	funcd(x1, &fl, &df, sf, df1, df2, p);
	funcd(x2, &fh, &df, sf, df1, df2, p);
	for(j = 0; fl*fh >= 0 && j < 10; j++) {
		x1 /= 2.0;		x2 *= 2.0;
		funcd(x1, &fl, &df, sf, df1, df2, p);
		funcd(x2, &fh, &df, sf, df1, df2, p);
		}
	if(fl*fh >= 0) {
		sprintf(info, "Value for inverse distribution\nmust be between %g and %g!", x1, x2);
		InfoBox(info);
		return 0.0;
		}
	if(fl < 0.0) {
		xl = x1;		xh = x2;
		}
	else {
		xh = x1;		xl = x2;
		swap = fl;		fl = fh;	fh = swap;
		}
	rts = x0;	dxold = fabs(x2-x1);	dx = dxold;
	funcd(rts, &f, &df, sf, df1, df2, p);
	for(j = 1; j <= 100; j++) {
        if((((rts-xh)*df-f)*((rts-xl)*df-f) >= 0.0) || (fabs(2.0*f) > fabs(dxold * df))) {
			dxold = dx;		dx = 0.5 * (xh-xl);		rts = xl + dx;
			if(xl == rts) return rts;
			}
		else {
			dxold = dx;		dx = f/df;		temp = rts;		rts -= dx;
			if(temp == rts) return rts;
			}
		if(fabs(dx) < _PREC) return rts;
		funcd(rts, &f, &df, sf, df1, df2, p);
		if(f < 0.0) {
			xl = rts;	fl = f;
			}
		else {
			xh = rts;	fh = f;
			}
		}
	return 0.0;
}

//---------------------------------------------------------------------------
//some statistical basics
//do quartiles, median of data
void d_quartile(int n, double *v, double *q1, double *q2, double *q3)
{
	int n2, n3;
	double f1, f2;

	if(!v || n<2) return;
	SortArray(n, v);			n2 = n >> 1;
	if(q1) {
		n3 = n2 >> 1;
		switch(n%4) {
		case 3:		n3 ++;		f1 = 2.0;		f2 = 2.0;		break;
		case 2:		n3 ++;		f1 = 3.0;		f2 = 1.0;		break;
		case 1:		n3 ++;		f1 = 4.0;		f2 = 0.0;		break;
		default:	f1 = 1.0;	f2 = 3.0;						break;
			}
		*q1 = (f1*v[n3-1] + f2*v[n3])/4.0;
		}
	if(q2) {
		if(n & 1) *q2 = v[n2];
		else *q2 = (v[n2-1] + v[n2])/2.0;
		}
	if(q3) {
		n3 = n2 >> 1;
		switch(n%4) {
		case 3:		n3++;		f1 = 2.0;		f2 = 2.0;	break;
		case 2:		f1 = 3.0;	f2 = 1.0;					break;
		case 1:		f1 = 4.0;	f2 = 0.0;					break;
		default:	f1 = 1.0;	f2 = 3.0;					break;
			}
		n3 += n2;
		*q3 = (f2*v[n3-1] + f1*v[n3])/4.0;
		}
}

//do arithmethic mean
double d_amean(int n, double *v)
{
	int i;
	double sum;

	for(i = 0, sum = 0.0; i < n; i++) {
		sum += (v[i]);
		}
	return (sum/n);
}


//do geometric mean
double d_gmean(int n, double *v)
{
	int i;
	double sum;

	for(i = 0, sum = 0.0; i < n; i++) {
		if(v[i] <= 0.0) return 0.0;
		sum += log(v[i]);
		}
	return exp(sum/n);
}

//do harmonic mean
double d_hmean(int n, double *v)
{
	int i;

	double sum;

	for(i = 0, sum = 0.0; i < n; i++) {
		if(v[i] == 0.0) return 0.0;
		sum += 1.0/(v[i]);
		}
	return (n/sum);
}

//---------------------------------------------------------------------------
// Pearsons linear correlation
// Ref.: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Rcipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 503 ff.
double d_pearson(double *x, double *y, int n, char *dest, DataObj *data)
{
	int j, r, c;
	double yt, xt, t, df, res[3];
	double syy=0.0, sxy=0.0, sxx=0.0, ay=0.0, ax=0.0;
	AccRange *rD;

	for(j = 0;	j < n; j++) {				// find means
		ax += x[j];			ay += y[j];
		}
	ax /= n;			ay /= n;
	for(j = 0; j < n; j++) {				// correlation
		xt = x[j] - ax;		yt = y[j] - ay;
		sxx += xt*xt;		syy += yt*yt;		sxy += xt * yt;
		}
	res[0] = sxy/sqrt(sxx*syy);		//pearsons r
	if(dest) {
		res[1] = 0.5 * log((1.0+res[0]+_PREC)/(1.0-res[0]+_PREC));	//Fishers z-transform
		df = n-2;
		t = res[0]*sqrt(df/((1.0-res[0]+_PREC)*(1.0+res[0]+_PREC)));	//Student's t
		res[2] = betai(0.5*df, 0.5, df/(df+t*t));					//probability
		}
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(j = 0; j < 3 && rD->GetNext(&c, &r); j++) {
			data->SetValue(r, c, res[j]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	return res[0];
}

//---------------------------------------------------------------------------
// Spearman rank-order correlation
// Ref.: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Recipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 507 ff.

//Given a sorted array w, crank replaces the elements by their rank
void crank(int n, double *w0, double *s)
{
	int j=1, ji, jt;
	double t, rank, *w = w0-1;

	*s = 0.0;
	while (j < n) {
		if(w[j+1] != w[j]) {
			w[j] = j;		++j;
			}
		else {
			for(jt = j+1; jt <= n; jt++)
				if(w[jt] != w[j]) break;
			rank = 0.5 * (j+jt-1);
			for(ji = j; ji <= (jt-1); ji++) w[ji] = rank;
			t = jt -j;
			*s += t*t*t -t;
			j = jt;
			}
		}
	if(j == n) w[n] = n;
}

//the actual rank correlation
double d_spearman(double *x, double *y, int n, char *dest, DataObj *data)
{
	int j, r, c;
	double vard, t, sg, sf, fac, en3n, en, df, aved, tmp;
	double res[5];
	AccRange *rD;

	SortArray2(n, x, y);		crank(n, x, &sf);
	SortArray2(n, y, x);		crank(n, y, &sg);
	for(j = 0, res[0] = 0.0; j < n; j++) res[0] += ((tmp = (x[j]-y[j]))*tmp);
	en = n;						en3n = en*en*en -en;
	aved = en3n/6.0 - (sf+sg)/12.0;
	fac = (1.0-sf/en3n)*(1.0-sg/en3n);
	vard = ((en-1.0)*en*en*((en+1.0)*(en+1.0))/36.0)*fac;
	vard = ((en-1.0)*en*en*((tmp = (en+1.0))*tmp)/36.0)*fac;
	res[1] = (res[0]-aved)/sqrt(vard);
	res[2] = errfc(fabs(res[1])/_SQRT2);
	res[3] = (1.0-(6.0/en3n)*(res[0]+0.5*(sf+sg)))/fac;
	t = res[3]*sqrt((en-2.0)/((res[3]+1.0)*(1.0-res[3])));
	df = en-2.0;
    res[4] = betai(0.5*df, 0.5, df/(df+t*t));
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(j = 0; j < 5 && rD->GetNext(&c, &r); j++) {
			data->SetValue(r, c, res[j]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	return res[3];
}

//linear regression
double d_regression(double *x, double *y, int n, char *dest, DataObj *data)
{
	double sx, sy, dx, dy, sxy, sxx, syy, sdy;
	double res[10];		// slope, intercept, mean x, mean y, SE of slope, 
						//   variance(x), variance(y), variance(fit), F of regression 
	int i, j, r, c;
	AccRange *rD;

	if(n < 2) return 0.0;
	for(i = 0, 	sx = sy = 0.0; i < n; i++) {
		sx += x[i];			sy += y[i];
		}
	res[2] = sx /n;			res[3] = sy/n;
	sxy = sxx = syy = 0.0;
	for(i = 0; i < n; i++) {
		dx = x[i]-res[2];	dy = y[i]-res[3];
		sxx += (dx*dx);		syy += (dy*dy);		sxy += (dx*dy);
		}
	res[0] = sxy / sxx;		res[1] = res[3] - res[0] * res[2];
	for(i = 0, sdy = 0.0; i < n; i++) {
		dy = y[i] - (res[1] + x[i] *res[0]);
		sdy += (dy * dy);
		}
	sdy = sdy/(n-2);		res[4] = sqrt(sdy/sxx);
	res[5] = sxx/(n-1);		res[6] = syy/(n-1);			res[7] = sdy;
	res[8] = sxy/sdy*sxy/sxx;
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(j = 0; j < 9 && rD->GetNext(&c, &r); j++) {
			data->SetValue(r, c, res[j]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	return n;
}

//t-test
double d_ttest(double *x, double *y, int n1, int n2, char *dest, DataObj *data)
{
	int i, r, c;
	double sx, sy, mx, my, d, df, p;
	double res[9];			// mean1, SD1, n1, mean2, SD2, n2, p if variances equal,
	AccRange *rD;			//    corrected df, corrected p

	for(i=0, sx = 0.0; i < n1; sx += x[i], i++);				mx = sx/n1;
	for(i=0, sy = 0.0; i < n2; sy += y[i], i++);				my = sy/n2;
	for(i=0, sx = 0.0; i < n1; sx += ((d=(x[i]-mx))*d), i++);
	for(i=0, sy = 0.0; i < n2; sy += ((d=(y[i]-my))*d), i++);
    d = ((sx+sy)/(n1+n2-2)) * ((double)(n1+n2)/(double)(n1*n2));
	d = (mx-my)/sqrt(d);										//Student's t

	//Welch's correction for differences in variance
	df = (sx/(double)n1)*(sx/(double)n1)/(double)(n1+1)+(sy/(double)n2)*(sy/(double)n2)/(double)(n2+1);
	df = (sx/(double)n1+sy/(double)n2)*(sx/(double)n1+sy/(double)n2)/df;
	df -= 2.0;
	p = betai(df/2.0, 0.5, (df/(df+d*d)));
	if((dest) && (data) && (rD = new AccRange(dest))) {
		res[0] = mx;	res[1] = sqrt(sx/(double)(n1-1));	res[2] = n1;
		res[3] = my;	res[4] = sqrt(sy/(double)(n2-1));	res[5] = n2;
		res[7] = df;	df = (n1-1) + (n2-1);	res[6] = betai(df/2.0, 0.5, (df/(df+d*d)));
		res[8] = p;
		rD->GetFirst(&c, &r);
		for(i = 0; i < 9 && rD->GetNext(&c, &r); i++) {
			data->SetValue(r, c, res[i]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	return p;
}

//f-test
double d_ftest(double *x, double *y, int n1, int n2, char *dest, DataObj *data)
{
	int i, r, c;
	double sx, sy, mx, my, d, df1, df2, p;
	double res[6];			// mean1, SD1, n1, mean2, SD2, n2
	AccRange *rD;

	for(i=0, sx = 0.0; i < n1; sx += x[i], i++);				mx = sx/n1;
	for(i=0, sy = 0.0; i < n2; sy += y[i], i++);				my = sy/n2;
	for(i=0, sx = 0.0; i < n1; sx += ((d=(x[i]-mx))*d), i++);	sx /= (n1-1);
	for(i=0, sy = 0.0; i < n2; sy += ((d=(y[i]-my))*d), i++);	sy /= (n2-1);
	d = sx/sy;		df1 = n1-1;		df2 = n2-1;
	p= 2.0 * betai(df2/2.0, df1/2.0, df2/(df2+df1*d));
	if((dest) && (data) && (rD = new AccRange(dest))) {
		res[0] = mx;	res[1] = sqrt(sx);	res[2] = n1;
		res[3] = my;	res[4] = sqrt(sy);	res[5] = n2;
		rD->GetFirst(&c, &r);
		for(i = 0; i < 6 && rD->GetNext(&c, &r); i++) {
			data->SetValue(r, c, res[i]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	return p;
}