mymcmc {Bhat} | R Documentation |

This function generates MCMC-based samples from a (posterior) density f (not necessarily normalized). It uses a Metropolis algorithm in conjunction with a multivariate normal proposal distribution which is updated adaptively by monitoring the correlations of succesive increments of at least 2 pilot chains. The method is described in De Gunst, Dewanji and Luebeck (submitted). The adaptive method is similar to the one proposed in Gelfand and Sahu (JCGS 3:261–276, 1994).

```
mymcmc(
x,
nlogf,
m1,
m2 = m1,
m3,
scl1 = 0.5,
scl2 = 2,
skip = 1,
covm = 0,
nfcn = 0,
plot = FALSE,
plot.range = 0
)
```

`x` |
a list with components 'label' (of mode character), 'est' (the parameter vector with the initial guess), 'low' (vector with lower bounds), and 'upp' (vector with upper bounds) |

`nlogf` |
negative log of the density function (not necessarily normalized) for which samples are to be obtained |

`m1` |
length of first pilot run (not used when covm supplied) |

`m2` |
length of second pilot run (not used when covm supplied ) |

`m3` |
length of final run |

`scl1` |
scale for covariance of mv normal proposal (second pilot run) |

`scl2` |
scale for covariance of mv normal proposal (final run) |

`skip` |
number of cycles skipped for graphical output |

`covm` |
covariance matrix for multivariate normal proposal distribution. If supplied, all pilot runs will be skipped and a run of length m3 will be produced. Useful to continue a simulation from a given point with specified covm |

`nfcn` |
number of function calls |

`plot` |
logical variable. If TRUE the chain and the negative log density (nlogf) is plotted. The first m1+m2 cycles are shown in green, other cycles in red |

`plot.range` |
[Not documented. Leave as default] |

standard output reports a summary of the acceptance fraction, the current values of nlogf and the parameters for every (100*skip) th cycle. Plotted chains show values only for every (skip) th cycle.

list with the following components:

`f` |
values of nlogf for the samples obtained |

`mcmc` |
the chain (samples obtained) |

`covm` |
current covariance matrix for mv normal proposal distribution |

This function is part of the Bhat exploration tool

E. Georg Luebeck (FHCRC)

too numerous to be listed here

```
# generate some Poisson counts on the fly
dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))
# neg. log-likelihood of Poisson model with 'linear-quadratic' mean:
nlogf <- function (x) {
ds <- data[, 1]
y <- data[, 2]
g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds))
return(sum(g - y * log(g)))
}
# start MCMC near mle
x <- list(label=c("a","b","c"), est=c(20, 0.5, 0.05), low=c(0,0,0), upp=c(100,10,.1))
# samples from posterior density (~exp(-nlogf))) with non-informative
# (random uniform) priors for "a", "b" and "c".
out <- mymcmc(x, nlogf, m1=2000, m2=2000, m3=10000, scl1=0.5, scl2=2, skip=10, plot=TRUE)
# start MCMC from some other point: e.g. try x$est <- c(16,.2,.02)
```

[Package *Bhat* version 0.9-12 Index]