/*
metrop.c
metropolis for N(0,1) based on Uniform candidates

gcc -shared -lgsl -lgslblasnative -o metrop.so metrop.c

dyn.load("metrop.so")
norm_function(n,alpha)
{
  .C("metrop",as.integer(n),as.double(alpha),as.double(1:n))[[3]]
}

*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>


double min(double a, double b)
{
  if (a < b)
    return(a);
  else
    return(b);   
}

void metrop(long *np, double *alphap, double *vec)
{
  long n,i;
  double alpha,x,can,u,aprob;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
  n=*np;alpha=*alphap;
  x=0.0;
  vec[0]=x;
  for (i=1;i<n;i++)
    {
      can=x+gsl_ran_flat(r,-alpha,alpha);
      aprob=min(1.0, gsl_ran_ugaussian_pdf(can)/gsl_ran_ugaussian_pdf(x));
      u=gsl_ran_flat(r,0.0,1.0);
      if (u < aprob)
	x=can;
      vec[i]=x;
    }
}


/* eof  */