/*********  A Markov Chain Monte Carlo (MCMC) Algorithm:
            Random-Walk Metropolis
Author: Fridtjof Thomas <fth@cts-fs1.du.se>
Date: 18 Jun 1998

Provided without performance guarantees for free, public,
non-commercial use.


There are many MCMC algorithms, not just one.
Here is code for "Random-Walk Metropolis".
The code will generate three figures, meant to
replicate Figure 1.1 (p. 6) in
  Gilks, Richardson, and Spiegelhalter (editors):
  Markov Chain Monte Carlo in Practice, Chapman&Hall, 1996.


first part:
proposal distribution N(x,0.5) - the chain mixes nicely.

second part:
proposal distribution N(x,0.1) - almost all steps are
accepted, but since the variance of the proposal
distribution is low, the chain mixes slowly.

third part: proposal distribution N(x,10) -
many steps are rejected and the chain mixes purely as a result.


A central question in MCMC is to determine c
onvergence of the Markov chain.
There is a site under construction
devoted to these issues
(by Mengerson and Robert):
 MCMC Diagnostics
 http://www.ensae.fr/crest/statistique/robert/McDiag/


Many texts about MCMC are available
from the MCMC Preprint Service:
 http://www.stats.bris.ac.uk/MCMC/


Have fun!


***********************************************************/

/***********************************************************
*                                                          *
*  Iterations from random-walk Metropolis algorithm.       *
*                                                          *
*  See: Gilks, Richardson, and Spiegelhalter:              *
*       MCMC in Practice, p.6.                             *
*                                                          *
*     (Fridtjof Thomas, 98-02-04)                          *
*                                                          *
*                                                          *
***********************************************************/



NEW;

/* defining a procedure for the calculations
   return syntax: (iteration number, value)    */


PROC(1)=RandomMetro(startValue,stdDev,numIter);
        LOCAL t,x,ser,y,u,rat,alpha,ap;
        t=0;
                x=startValue;
                ser=t~x~2~-2;
                DO WHILE t<numIter;
                        y=x+RNDN(1,1)*stdDev;
                        u=RNDU(1,1);
                        rat=PDFN(y)/PDFN(x);
                        alpha=MINC(1|rat);

                        IF u<=alpha; ap=y;
                                ELSE; ap=x;
                        ENDIF;

                        x=ap;

                        t=t+1;
                        ser=ser|(t~ap~2~-2);
                ENDO;
        RETP(ser);
ENDP;


/* proc generating data */

PROC(3)=Figure1_1(startVal1,startVal2,startVal3,numIter);
        LOCAL data1,data2,data3;
                {data1}=RandomMetro(startVal1,0.5,numIter);
                {data2}=RandomMetro(startVal2,0.1,numIter);
                {data3}=RandomMetro(startVal3,10,numIter);
        RETP(data1,data2,data3);
ENDP;


/* generating data */

        /* change start values and number of iterations here!!! */
{data1,data2,data3}=Figure1_1(-5,-5,-5,1000);




/* graphic output */

LIBRARY PGRAPH;
GRAPHSET;

BEGWIND;

        WINDOW(3,0,0);

        _pcolor={2,4,4};
        _pltype={6,3,3};

        SETWIND(1);
                TITLE("Proposal distribution N(X,0.5)");
                XLABEL("iteration t");
                YLABEL("X");
                XY(data1[.,1],data1[.,2:4]);
        NEXTWIND;
                TITLE("Proposal distribution N(X,0.1)");
                XLABEL("iteration t");
                YLABEL("X");
                XY(data2[.,1],data2[.,2:4]);
        NEXTWIND;
                TITLE("Proposal distribution N(X,10)");
                XLABEL("iteration t");
                YLABEL("X");
                XY(data3[.,1],data3[.,2:4]);

ENDWIND;


END;