MCMCmetrop1R            package:MCMCpack            R Documentation

_M_e_t_r_o_p_o_l_i_s _S_a_m_p_l_i_n_g _f_r_o_m _U_s_e_r-_W_r_i_t_t_e_n _R _f_u_n_c_t_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     This function allows a user to construct a sample from a
     user-defined R function using a random walk Metropolis algorithm.
     It assumes the parameters to be sampled are in a single block.

_U_s_a_g_e:

     MCMCmetrop1R(fun, theta.init, burnin = 500, mcmc = 20000, thin = 1,
                  tune = 1, verbose = 0, seed=NA, logfun = TRUE,
                  optim.trace = 0, optim.REPORT = 10, optim.maxit = 500, ...)

_A_r_g_u_m_e_n_t_s:

     fun: The (log)density from which to take a sample. This should be
          a function defined in R and it should take a single argument.
          Additional arguments can be passed implicitly by either
          putting them in the global environment or by passing them as
          additional arguments to 'MCMCmetrop1R()'. See the Details
          section  and the examples below. 

theta.init: Starting values for the sampling. Must be of the
          appropriate dimension.

  burnin: The number of burn-in iterations for the sampler.

    mcmc: The number of MCMC iterations after burnin.

    thin: The thinning interval used in the simulation.  The number of
          MCMC iterations must be divisible by this value.

    tune: Can be either a positive scalar or a k-vector, where k is the
          length of theta.

 verbose: A switch which determines whether or not the progress of the
          sampler is printed to the screen.  If 'verbose' is greater
          than 0 the iteration number, the theta vector, the function
          value, and the Metropolis acceptance rate are sent to the
          screen every 'verbose'th iteration.

    seed: The seed for the random number generator.  If NA, the
          Mersenne Twister generator is used with default seed 12345;
          if an integer is  passed it is used to seed the Mersenne
          twister.  The user can also pass a list of length two to use
          the L'Ecuyer random number generator, which is suitable for
          parallel computation.  The first element of the list is the
          L'Ecuyer seed, which is a vector of length six or NA (if NA 
          a default seed of 'rep(12345,6)' is used).  The second
          element of  list is a positive substream number. See the
          MCMCpack  specification for more details.

  logfun: Logical indicating whether 'fun' returns the natural log of a
          density function (TRUE) or a density (FALSE).

optim.trace: The value of the 'trace' parameter sent to 'optim' during
          an initial maximization of 'fun'. 

optim.REPORT: The value of the 'REPORT' parameter sent to 'optim'
          during an initial maximization of 'fun'.

optim.maxit: The value of the 'maxit' parameter sent to 'optim' during
          an initial maximization of 'fun'.

     ...: Additional arguments.

_D_e_t_a_i_l_s:

     MCMCmetrop1R produces a sample from a user-defined (log)density
     using a random walk Metropolis algorithm with multivariate normal
     proposal distribution. See Gelman et al. (2003) and Robert &
     Casella (2004) for details of the random walk Metropolis
     algorithm.

     The proposal distribution is centered at the current value of
     theta and has variance-covariance V = T (-1*H)^{-1} T, where T is
     a the diagonal positive definite matrix formed from the 'tune' and
     H is the approximate Hessian of 'fun' evaluated at it's mode. This
     last calculation is done via an initial call to 'optim'.

     Passing data into a (log)density function does not conform to
     standard R  practice.  The (log)density function should have only
     one argument - the  parameter vector.  The data used by the
     (log)density need to be either in the  global environment or
     defined in the call to 'MCMCmetrop1R()', which will  create an
     environment that contains the data and then call the (log)density
     function in that environment.  The names of the data objects  are
     not checked, so be careful that these match the names used with
     the  (log)density function.  _NOTE:_ This interface will likely be
     changed in  a future implementation of MCMCpack.

_V_a_l_u_e:

     An mcmc object that contains the posterior density sample.  This 
     object can be summarized by functions provided by the coda
     package.

_R_e_f_e_r_e_n_c_e_s:

     Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald B. Rubin.
     2003. _Bayesian Data Analysis_. 2nd Edition. Boca Raton: Chapman &
     Hall/CRC.

     Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein.  2004.  
     _Scythe Statistical Library 1.0._ <URL: http://scythe.wustl.edu>.

     Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
     _Output Analysis and Diagnostics for MCMC (CODA)_. <URL:
     http://www-fis.iarc.fr/coda/>.

     Christian P. Robert and George Casella. 2004. _Monte Carlo
     Statistical Methods_. 2nd Edition. New York: Springer.

_S_e_e _A_l_s_o:

     'plot.mcmc', 'summary.mcmc', 'optim'

_E_x_a_m_p_l_e_s:

       ## Not run: 
         
         ## logistic regression with an improper uniform prior
         ## X and y are passed as args to MCMCmetrop1R

         logitfun <- function(beta){
           eta <- X %*% beta
           p <- 1.0/(1.0+exp(-eta))
           sum( y * log(p) + (1-y)*log(1-p) )
         }
        
         x1 <- rnorm(1000)
         x2 <- rnorm(1000)
         Xdata <- cbind(1,x1,x2)
         p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2))
         yvector <- rbinom(1000, 1, p)
         
         post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0),
                                   X=Xdata, y=yvector,
                                   thin=1, mcmc=40000, burnin=500,
                                   tune=c(1.5, 1.5, 1.5),
                                   verbose=500, logfun=TRUE, optim.maxit=100)

         raftery.diag(post.samp)
         plot(post.samp)
         summary(post.samp)
         ## ##################################################
                                   
         
         ## logistic regression with an improper uniform prior
         ## X and y are now in the global environment

         logitfun <- function(beta){
           eta <- X %*% beta
           p <- 1.0/(1.0+exp(-eta))
           sum( y * log(p) + (1-y)*log(1-p) )
         }
         
         x1 <- rnorm(1000)
         x2 <- rnorm(1000)
         X <- cbind(1,x1,x2)
         p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2))
         y <- rbinom(1000, 1, p)
         
         post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0),
                                   thin=1, mcmc=40000, burnin=500,
                                   tune=c(1.5, 1.5, 1.5),
                                   verbose=500, logfun=TRUE, optim.maxit=100)

         raftery.diag(post.samp)
         plot(post.samp)
         summary(post.samp)
         ## ##################################################

         ##  negative binomial regression with an improper unform prior
         ## X and y are passed as args to MCMCmetrop1R
         negbinfun <- function(theta){
           k <- length(theta)
           beta <- theta[1:(k-1)]
           alpha <- exp(theta[k])
           mu <- exp(X %*% beta)
           log.like <- sum(
                           lgamma(y+alpha) - lfactorial(y) - lgamma(alpha) +
                           alpha * log(alpha/(alpha+mu)) +
                           y * log(mu/(alpha+mu))
                          )
         }

         n <- 1000
         x1 <- rnorm(n)
         x2 <- rnorm(n)
         XX <- cbind(1,x1,x2)
         mu <- exp(1.5+x1+2*x2)*rgamma(n,1)
         yy <- rpois(n, mu)

         post.samp <- MCMCmetrop1R(negbinfun, theta.init=c(0,0,0,0), y=yy, X=XX,
                                   thin=1, mcmc=35000, burnin=1000,
                                   tune=1.5, verbose=500, logfun=TRUE,
                                   optim.maxit=500, seed=list(NA,1))
         raftery.diag(post.samp)
         plot(post.samp)
         summary(post.samp)
         ## ##################################################
      
       ## End(Not run)

