/* Author: Murali Haran, Statistics Department, Penn State University */ /* Version: November, 2006 */ /* A C function for computing consistent batch means estimate of standard error from: */ /* Citation: Galin L. Jones, Murali Haran, Brian S. Caffo, and Ronald Neath, "Fixed-Width Output Analysis for Markov Chain Monte Carlo" (2006), Journal of the American Statistical Association, 101:1537--1547 */ /* To compile: gcc -o cbm cbm.c -lm */ /* To run: ./cbm mcmclen numcol */ /* file is assumed to contain a matrix of output from MCMC */ /* with mcmclen number of rows(chain length) and numcol number of columns (parameters) */ #include #include #include #include #define VERBOSE 0 /* 1=> print descriptions of output*/ /* assume that the vector values are g(x_1),..,g(x_n) */ /* and expectation of interest is E(g(x)) */ /* n is the length of the vector, and vec is the vector */ /* bs is the batchsize. Default is sqrt(n), if bs==0, use cuberoot of n */ int cbm(int n,double *vec,int bs,double *bmres) { double est,serr,sigmahatsq; int a,b,Nbm; /* batch number and batch size, product */ int i,k; double *Y; if (bs) /* if sqroot (default) */ { b=(int ) sqrt(n); a=n/b; /* use 'div' operator */ } else /* cuberoot */ { b=(int ) pow(n,1/3); a=n/b; /* use 'div' operator */ } Y=(double *)malloc(a*sizeof(double)); /* memory allocation for batch means */ /* computing batch means according to batch sizes */ est=0.0; for (k=0;k mcmclen numcol\n",argv[0]); exit(1); } else { strlength = strlen(argv[1]); inputpath = malloc(strlength*sizeof(char)+1); /* CHECK this added 1 to avoid strange errors*/ sprintf(inputpath,"%s",argv[1]); mcmclen = atoi(argv[2]); numcol = atoi(argv[3]); } /* allocate memory according to length of MCMC run and number of columns */ /* X[mcmclen] */ X=(double *) malloc(mcmclen*sizeof(double)); /* mcmc[numcol][mcmclen] */ mcmc = (double **) malloc(numcol*sizeof(double)); for (i=0;i