[R-gregmisc-users] SF.net SVN: r-gregmisc: [1057] trunk/RMCMC/R/mh.R

[<wa...@us...>]

Revision: 1057
          http://svn.sourceforge.net/r-gregmisc/?rev=1057&view=rev
Author:   warnes
Date:     2007-02-16 15:11:53 -0800 (Fri, 16 Feb 2007)

Log Message:
-----------
Add initial code for mh() to generate a Metropolis-Hastings sampler

Added Paths:
-----------
    trunk/RMCMC/R/mh.R

Added: trunk/RMCMC/R/mh.R
===================================================================
--- trunk/RMCMC/R/mh.R	                        (rev 0)
+++ trunk/RMCMC/R/mh.R	2007-02-16 23:11:53 UTC (rev 1057)
@@ -0,0 +1,96 @@
+## return a function that generates a Metropolis-Hastings MCMC sampler
+## using the provided posterior and proposal.  Extra arguments will be
+## appended to the call to the proposal.
+
+mh <- function( posterior, qgenerator, qdensity, x0, ..., log=FALSE)
+  {
+    npar <- unlist(length(x0))
+    # ensure everything works, so we show errors early on
+    cat("Testing that the posterior density function works properly on x0...")
+    posterior(unlist(x0))
+    cat("Success.\n")
+
+    cat("Testing that the proposal generator function works properly on x0...")
+    qgenerator(unlist(x0))
+    cat("Success.\n")
+
+    cat("Testing that the proposal generator function works properly on x0...")
+    qdensity(unlist(x0), unlist(x0))
+    cat("Success.\n")
+
+    constructor <- match.call()
+    
+    fun <- function(x0, niter)
+      {
+        ## test parameters
+        if(niter <1) niter <- 0
+        
+        ## step 0: setup
+           X <- matrix( NA, nrow=niter+1, ncol=npar ) # sampler state
+           X[1,] <- unlist(x0)                        # initial state
+           # note that X is 1 element longer than the number of iterations
+
+           Y <- matrix( NA, nrow=niter+1, ncol=npar ) # proposed state
+
+           p.X <- numeric(length=npar)  # p(X): posterior density at X
+           p.Y <- numeric(length=npar)  # p(Y): posterior density at Y
+
+           q.Y.X <- numeric(length=npar) # q(X,Y):proposal density from X to Y 
+           q.X.Y <- numeric(length=npar) # q(Y,X):proposal density from Y to X
+
+           alpha <- numeric(length=npar) # alpha(X,Y): trans. prob form X to Y
+
+
+        for(t in 1:niter+1)
+          {
+            ## step 1: generate Y
+            Y[t] <- proposal(X[t], ...)
+
+            ## step 2: calculate transition probability
+            if(log)
+              {
+                p.X[t] <- posterior(X[t], log=TRUE)
+                p.Y[t] <- posterior(Y[t], log=TRUE)
+            
+                q.X.Y[y] <- qdensity( X[t], Y[t], log=TRUE)
+                q.Y.X[y] <- qdensity( Y[t], X[t], log=TRUE)
+
+                alpha <- max(1, exp(p.X[t] - p.Y[t] + q.Y.X[t] - q.X.Y[t]))
+              }
+            else
+              {
+                p.X[t] <- posterior(X[t])
+                p.Y[t] <- posterior(Y[t])
+            
+                q.X.Y[t] <- qdensity( X[t], Y[t])
+                q.Y.X[t] <- qdensity( Y[t], X[t])
+
+                alpha[t] <- max(1, p.X[t] / p.Y[t] * q.Y.X[t] / q.X.Y[t])
+              }
+
+            ## step 3: flip a coin to decide whether X_t or Y is the next state
+            u <- runif()
+            if( u < alpha[t] )
+              X[t+1] <- Y
+            else
+              X[t+1] <- X[t]
+          }
+              
+        retval <- list(
+                       niter=niter,
+                       X = X,
+                       Y = Y,
+                       p.X = p.X,
+                       p.Y = p.Y,
+                       q.Y.X = q.Y.X,
+                       q.X.Y = q.X.Y,
+                       alpha = alpha,
+                       call = match.call(),
+                       constructor=constructor,
+                       ) 
+        
+        
+      }
+
+    
+  }


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.