myutil.c

MCMC（马尔可夫－盟特卡罗方法）实现的程序

static double posarg(double);
static double negarg(double);
static double asform(double);

double lgamma(double arg)
{
	signgamR1 = 1.0;
	if (arg <= 0.0)
		return (negarg(arg));
	if (arg > 8.0)
		return (asform(arg));
	return (log(posarg(arg)));
}

/* Equation 6.1.41 Abramowitz and Stegun */
/* See also ACM algorithm 291 */

static double asform(arg)
double arg;
{
	double log();
	double n, argsq;
	int i;

	argsq = 1. / (arg * arg);
	for (n = 0, i = M - 1; i >= 0; i--) {
		n = n * argsq + p1[i];
	}
	return ((arg - .5) * log(arg) - arg + hl2pi + n / arg);
}

static double negarg(arg)
double arg;
{
	double temp;
	double log(), sin(), posarg();

	arg = -arg;
	temp = sin(xpi * arg);
	if (temp == 0.0)
		/* removed DOMAIN_ERROR */ printerr("negarg: temp == 0.0");
	if (temp < 0.0)
		temp = -temp;
	else
		signgamR1 = -1;
	return (-log(arg * posarg(arg) * temp / xpi));
}

static double posarg(arg)
double arg;
{
	double n, d, s;
	register i;

	if (arg < 2.0)
		return (posarg(arg + 1.0) / arg);
	if (arg > 3.0)
		return ((arg - 1.0) * posarg(arg - 1.0));

	s = arg - 2.;
	for (n = 0, d = 0, i = N - 1; i >= 0; i--) {
		n = n * s + p2[i];
		d = d * s + q2[i];
	}
	return (n / d);
}
 
/* end of lgamma insertion */


double lfactl(int n)
{
	static double lookup[1000];
	if(n < 0)printerr("lfactl: n < 0");
	if(n <= 1)return 0.0;
/*	return lgamma(n+1.0); */
	if(n >= 1000)return lgamma(n+1.0);
 	if(lookup[n] == 0.0)lookup[n] = lgamma(n+1.0);
	return lookup[n];  
}


double expdev(void)
{

      int k;
      ++jindic;
      if(jindic>97)jindic=0;
      k = jindic+27;
      if(k>97)k=k-98;
      rand_table[jindic] = rand_table[jindic]^rand_table[k];
	return(-log(((unsigned)rand_table[jindic] + 1.0)/4294967298.0));
}

void opengfsr(void)
{
	
	FILE 	*rt;
	int 	j;
	
	rt = fopen("INTFILE","r");
	if(rt==NULL) {
		printf("I need INTFILE! Where is INTFILE?\n");
		exit(1);
	}
	for(j=0;j<98;++j) {fscanf(rt,"%d",&rand_table[j]);
	//printf("ASDF%d\n", rand_table[j]);/*TEST*/
	}
	fscanf(rt,"%d",&jindic); 
	fclose(rt); 
}

void closegfsr(void)
{
	FILE 	*rt;
	int 	j;
	
	rt = fopen("INTFILE","w");
	for(j=0;j<98;++j)fprintf(rt,"%d\n",rand_table[j]);
	fprintf(rt,"%d\n",jindic);
	fclose(rt);
}

float norm4(void)
{
	static int iset=0;
	static float gset;
	float fac,rsq,v1,v2;

	if  (iset == 0) {
		do {
			v1=2.0*gfsr4()-1.0;
			v2=2.0*gfsr4()-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}


double norm8(void)
{
	static int iset=0;
	static double gset;
	double fac,rsq,v1,v2;

	if  (iset == 0) {
		do {
			v1=2.0*gfsr8()-1.0;
			v2=2.0*gfsr8()-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}

void mom(float x[],int n,float *x1,float *x2,float *x3,float *x4,
		float *min,float *max)
{
	int i;
	double s1,s2,s3,s4,an,an1,dx,dx2,xi,var;
	
	s1 = x[0];
	s2 = 0.0;
	s3 = 0.0;
	s4 = 0.0;
	*min = s1;
	*max = s1;
	for(i=1;i<n;++i){
		xi = x[i];
		an = i+1;
		an1 = i;
		dx = (xi-s1)/an;
		dx2 = dx*dx;
		s4 -= dx*(4.0*s3-dx*(6.0*s2+an1*(1.0+pow(an1,3.0))*dx2));
		s3 -= dx*(3.0*s2-an*an1*(an-2.0)*dx2);
		s2 += an*an1*dx2;
		s1 += dx;
		if(xi<*min)*min=xi;
		if(xi>*max)*max=xi;
	}
	*x1 = s1;
	var = n>1 ? s2/(n-1) : 0.0;
	*x2 = sqrt(var);
	*x3 = var>0.0 ? 1.0/(n-1)*s3/pow(var,1.5) : 0.0;
	*x4 = var>0.0 ? 1.0/(n-1)*s4/pow(var,2.0)-3.0 : 0.0;
	return;
}

void isort(char dir,int n,int * x)  /* This is adapted from R 0.16 */
{
	int i, j, h, asc;
	int xtmp;

	if(dir == 'a' || dir == 'A')asc = 1;
	else asc = 0;
	h = 1;
	do {
		h = 3 * h + 1;
	}
	while (h <= n);

	do {
		h = h / 3;
		for (i = h; i < n; i++) {
			xtmp = x[i];
			j = i;
			if(asc){
				while (x[j - h] > xtmp) {
					x[j] = x[j - h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			else{
				while (x[j - h] < xtmp) {
					x[j] = x[j - h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			
		end:	x[j] = xtmp;
		}
	} while (h != 1);
}

void fsort(char dir,int n,float * x)  /* This is adapted from R 0.16 */
{
	int i, j, h, asc;
	float xtmp;

	if(dir == 'a' || dir == 'A')asc = 1;
	else asc = 0;
	h = 1;
	do {
		h = 3 * h + 1;
	}
	while (h <= n);

	do {
		h = h / 3;
		for (i = h; i < n; i++) {
			xtmp = x[i];
			j = i;
			if(asc){
				while (x[j - h] > xtmp) {
					x[j] = x[j - h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			else{
				while (x[j - h] < xtmp) {
					x[j] = x[j - h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			
		end:	x[j] = xtmp;
		}
	} while (h != 1);
}

void dsort(char dir,int n,double * x)  /* This is adapted from R 0.16 */
{
	int i, j, h, asc;
	double xtmp;

	if(dir == 'a' || dir == 'A')asc = 1;
	else asc = 0;
	h = 1;
	do {
		h = 3 * h + 1;
	}
	while (h <= n);

	do {
		h = h / 3;
		for (i = h; i < n; i++) {
			xtmp = x[i];
			j = i;
			if(asc){
				while (x[j - h] > xtmp) {
					x[j] = x[j - h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			else{
				while (x[j - h] < xtmp) {
					x[j] = x[j - h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			
		end:	x[j] = xtmp;
		}
	} while (h != 1);
}


void isorti(char dir, int n, int * x, int *indx) 
{
	int i, j, h, asc,indtmp;
	int xtmp,*priv;

	priv = (int *)malloc(n*sizeof(int));
	for(j=0;j<n;++j)priv[j] = x[j];
	if(dir == 'a' || dir == 'A')asc = 1;
	else asc = 0;
	for(j=0;j<n;++j)indx[j] = j;
	h = 1;
	do {
		h = 3 * h + 1;
	}
	while (h <= n);

	do {
		h = h / 3;
		for (i = h; i < n; i++) {
			xtmp = priv[i];
			indtmp = indx[i];
			j = i;
			if(asc){
				while (priv[j - h] > xtmp) {
					priv[j] = priv[j - h];
					indx[j] = indx[j-h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			else{
				while (priv[j - h] < xtmp) {
					priv[j] = priv[j - h];
					indx[j] = indx[j-h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			
		end:	priv[j] = xtmp;indx[j] = indtmp;
		}
	} while (h != 1);
	free(priv);
}

void isorti2(char dir, int n, int *x, int *priv, int *indx) 
{
	int i, j, h, asc,indtmp;
	int xtmp;

	for(j=0;j<n;++j)priv[j] = x[j];
	if(dir == 'a' || dir == 'A')asc = 1;
	else asc = 0;
	for(j=0;j<n;++j)indx[j] = j;
	h = 1;
	do {
		h = 3 * h + 1;
	}
	while (h <= n);

	do {
		h = h / 3;
		for (i = h; i < n; i++) {
			xtmp = priv[i];
			indtmp = indx[i];
			j = i;
			if(asc){
				while (priv[j - h] > xtmp) {
					priv[j] = priv[j - h];
					indx[j] = indx[j-h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			else{
				while (priv[j - h] < xtmp) {
					priv[j] = priv[j - h];
					indx[j] = indx[j-h];
					j = j - h;
					if (j < h)
						goto end;
				}
			}
			
		end:	priv[j] = xtmp;indx[j] = indtmp;
		}
	} while (h != 1);
}