mog1.c

he AutoMix package is a C program for Unix-like systems, implementing the automatic reversible jum

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#include "mog1.h"

#define tpi 6.283185307179586477

double my_infinity(void) {
  double zero = 0;
  return 1.0/zero;
}
double my_nan(void) {
  double zero = 0;
  return zero/zero;
}
#define INFINITY my_infinity()
#define NAN my_nan()

#define KMAX 3

const double hyper_mean = 0.0;
const double hyper_prec = 1.0;
const double hyper_a = 2.0;
const double hyper_b = 1.0;
const double hyper_delta = 1.0;
double hyper_w[KMAX];

double normrnd(double mean, double prec);
double gammarnd(double a, double b);
void dirichletrnd(int dim, double *param, double *sample);

double normprob(double mean, double var, double x);

double lognormprob(double mean, double var, double x);
double loggammaprob( double a, double b, double x );
double logdirichletprob(int dim, double *param, double *data);

double boxm();

extern double loggamma(double s);
extern double sdrand();
extern double rgamma(double s);

/* Function to return number of models */
void getkmax(int *kmax){
	*kmax = KMAX;
	return;
}

/* Function to return the dimension of each model */ 
/* each model has nc componenets, each having a mean and variance */
/* and nc-1 mixing coefficients (sum-to-1 constraint) */
void getnk(int kmax,int *nk){

	int k;
	int nc; // # mixture componenents
	for(k=0;k<kmax;k++){
		nc = k+1;
		nk[k] = 2*nc + nc-1;
	}
	return;

}

/* Function to return initial conditions for RWM runs */
/* sample them from the prior */
/*
	int i;
	FILE* fp = fopen("blah.txt","wt");

	for( i=0; i<10000; i++ ) 
		fprintf(fp, "%f,", normrnd(0,1.0/2.0) );

	fprintf(fp, "%f", normrnd(0,1.0/2.0) );

	fclose(fp);

	exit(1);
*/

void getic(int k, int nkk, double *theta){
	int ci;
	int nc;
	double pmean, pprec;
	double pweight[KMAX];

	nc = k+1;

	for( ci=0; ci<nc; ci++ ) 
		hyper_w[ci] = hyper_delta;

	dirichletrnd(nc, hyper_w, pweight);

	for( ci=0; ci<nc; ci++ ) {
		pprec = gammarnd(hyper_a, 1.0/hyper_b);
		pmean = normrnd(hyper_mean, 1.0/(hyper_prec*pprec) );

		theta[ci*3+0] = pmean;
		theta[ci*3+1] = pprec;

		if( ci<nc-1 )
			theta[ci*3+2] = pweight[ci];
	}
}

/* Function to return log target distribution up to additive const at (k,theta)
   value also returned in llh */

double weight[KMAX], mean[KMAX], prec[KMAX];

double lpost(int k,int nkk,double *theta, double *llh){ 	 

	int di, nd;
	int ci, nc;

	double sum;
	double llikelihood, lprior;

	nc = k+1;

	// log prior contribution
	lprior = 0;

	sum = 0;
	// load the current parameters
	for( ci=0; ci<nc; ci++ ) {
		hyper_w[ci] = hyper_delta;

		mean[ci] = theta[ci*3+0];
		prec[ci] = theta[ci*3+1];

		if( ci<nc-1 ) {
			weight[ci] = theta[ci*3+2];
			sum += weight[ci];
		} else {
			weight[ci] = 1.0-sum;
		}

		if( nc>1 && ((weight[ci]<=0.0)||(weight[ci]>=1.0) ) ){
			return -INFINITY; // RWM proposed impossible weights
		}

		if( prec[ci]<=0 ) {
			return -INFINITY; // RWM proposed impossible precision
		}
	}
	// prior on mixing coefficients
	lprior += logdirichletprob(nc, hyper_w, weight );
	for( ci=0; ci<nc; ci++ ) {
		// prior on precisions
		lprior += loggammaprob( hyper_a, 1.0/hyper_b, prec[ci] );

		// prior on means
		lprior += lognormprob( hyper_mean, 1.0/(hyper_prec*prec[ci]), mean[ci] );
	}

	// log likelihood contribution
	llikelihood = 0;

	nd = nData;
	for( di=0; di<nd; di++ ) {
		double ldata = 0;
		double coeff[KMAX];
		double exponent[KMAX];
		double maxExponent = -INFINITY;
		double arg;

		for( ci=0; ci<nc; ci++ ) {
			arg = (data[di]-mean[ci]);
			exponent[ci] = -arg*arg/(2.0)*prec[ci];
			coeff[ci] = weight[ci]/sqrt(tpi/prec[ci]);

			if(exponent[ci]>maxExponent) maxExponent = exponent[ci];
		}
		for( ci=0; ci<nc; ci++ ) {
			ldata += coeff[ci]*exp( exponent[ci] - maxExponent );
		}
		ldata = log(ldata) + maxExponent;

		llikelihood += ldata;
	}

	return llikelihood + lprior;
} 

double normrnd(double mean, double var) {
	return mean + boxm()*sqrt(var);
}

double gammarnd(double a, double b) {
	return rgamma(a)*b;
}

void dirichletrnd(int dim, double *param, double *sample) {

	int di;
	double sum = 0.0;
	for(di=0; di<dim; di++) {
		sample[di] = rgamma(param[di]);
		sum += sample[di];
	}
	for(di=0; di<dim; di++) {
		sample[di] /= sum;
	}
}

double normprob(double mean, double var, double x) {
	double arg = x-mean;
	return exp(-arg*arg/(2.0*var)) / sqrt(tpi*var);
}

double lognormprob(double mean, double var, double x) {
	double arg = x-mean;
	return -arg*arg/(2.0*var) - log(sqrt(tpi*var));
}

double loggammaprob( double a, double b, double x ) {
	 return -( a*log(b) + loggamma(a) ) + (a-1.0)*log(x) -x/b;
}

double logdirichletprob(int dim, double *param, double *data) {
	
	int di;
	double pr = 0;
	double paramSum = 0;
	double gammaSum = 0;
	for( di=0; di<dim; di++ ) {
		pr += (param[di]-1.0)*log( data[di] );

		paramSum += param[di];
		gammaSum += loggamma(param[di]);
	}
	pr += loggamma(paramSum) - gammaSum;
	
	return pr;
}

double boxm(){

  /*Function for returning single N(0,1) variable */

  double u1,u2,out;

  u1=sdrand();
  u2=sdrand();
  u1=tpi*u1;
  u2=sqrt(-2.0*log(u2));
  out=u2*cos(u1);

  return out;

}