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

/* Use the Metropolis-Hastings algorithm to sample from a
   d-dimensional multivariate normal N(0,I) distribution.
   
   A random walk proposal distribution is used with a range of
   different step sizes.

   The goal is to assess the performance of different step sizes for
   estimating E X(1) = 0 based on the sample mean from a chain of a
   given length.

   The process is replicated to get an estimate of the MSE for each
   step length.

   According to theory, an acceptance rate around 0.23 should be
   optimal.   
*/

int main(int argc, char** argv)
{
  /* The standard deviation of the step. */
  double f;
  
  /* The current point in the Markov Chain. */
  double* X;

  /* The current sum of the entire history of the Markov Chain. */
  double M;
  
  /* The current MSE estimate for X(1)^2 * ... * X(d)^2. */
  double MSE;

  /* The proposal point. */
  double* x;

  /* The dimension. */
  int d = 10;

  /* The length of the Markov chain realization. */
  int chainlen = 1000;

  /* The number of replications to estimate the MSE. */
  int nrep = 1000;

  /* Keep track of the number of acceptances. */
  int acc;

  /* Temporary variables. */
  int r, it, j;
  double R;
  
  /* Allocations. */
  X = (double*)malloc(d*sizeof(double));
  x = (double*)malloc(d*sizeof(double));

  /* Consider a range of different step lengths. */
  for(f=0.3; f<=1.5; f+=0.1)
    {
      /* Replication to estimate the MSE. */ 
      MSE = 0;
      acc = 0;
      for (r=0; r<nrep; ++r)
	{
	  /* Start at the mean (no burn-in). */
	  memset(X, 0, d*sizeof(double));
    
	  /* Generate one realization from the Markov chain. */
	  M = 0;
	  for (it=0; it<chainlen; ++it)
	    {
	      /* The proposal point. */
	      for (j=0; j<d; ++j)
		{
		  x[j] = X[j] + f*sqrt(-2*log(drand48())) * 
		    cos(2*M_PI*drand48());
		}

	      /* The Metropolis-Hastings ratio. */
	      R = 0;
	      for (j=0; j<d; ++j) R += 0.5*(X[j]*X[j] - x[j]*x[j]);
	      
	      /* Make a random decision to stay or move. */
	      if (log(drand48()) < R) 
		{
		  memcpy(X, x, d*sizeof(double));
		  ++acc;
		}

	      /* Update the average. */
	      M += X[0];
	    }
    
	  M /= chainlen;

	  /* Update the squared error for the sample mean. */
	  MSE += M*M;
	}

      /* Print out the step size, the MSE, and the acceptance rate. */
      MSE /= nrep;
      printf("%f\t%f\t%f\n", f, MSE, (double)acc/(double)(nrep*chainlen));
    }

  return 0;
}