# MetropolisHastings: Metropolis-Hastings sampler In bbricks: Bayesian Methods and Graphical Model Structures for Statistical Modeling

## Description

sample xhat from the target distribution p(xhat), with given proposal distribution q(xhat|x), the acceptance rate is:
min( 1 , dp(xhat)/dp(x) * dq(x|xhat)dq(xhat|x) )
where dp() is the density function of the target distribution, dq() the the density function of the proposal distribution. A new sample xhat is drawn from the sampler of the proposal distribution rq(). See examples.

## Usage

 1 MetropolisHastings(nsamples, xini, dp, dq, rq)

## Arguments

 nsamples integer, number of samples to draw xini initial sample, the chain of samples starts from here. xini must be a matrix of one row, or a numeric vector that will be converted to a matrix of one row. dp function(x), the LOG density function of the target distribution log dp(x), DON'T FORGET THE LOG. dq function(xhat,x), the LOG density function of the proposal distribution log dq(xhat | x), DON'T FORGET THE LOG. rq function(x), the generator of the proposal distribution rq(xhat | x).

## Value

a matrix of nsamples rows.

## Examples

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 ## example1: independent Metropolis-Hastings algorithm, get 5000 samples from Beta(2.7,6.3) ## with independent uniform proposal U(0,1), and independent normal proposal N(0.5,1). ## step1: define p() and q() dp <- function(x) if(x>0&x<1) dbeta(x,2.7,6.3,log = TRUE) else -Inf dq1 <- function(xnew,x) 0 #uniform proposal log density dq2 <- function(xnew,x) dnorm(xnew,0.5,1) #normal proposal log density rq1 <- function(x) runif(1,0,1) #uniform proposal sampler rq2 <- function(x) rnorm(1,0.5,1) #normal proposal sampler ## step2: get 5000 samples, with two different proposals X1 <- MetropolisHastings(nsamples = 5000,xini = runif(1,0,1),dp=dp,dq=dq1,rq=rq1) X2 <- MetropolisHastings(nsamples = 5000,xini = runif(1,0,1),dp=dp,dq=dq2,rq=rq2) ## step3: plot the result, calculate acceptance rate sum(diff(X1)!=0)/nrow(X1) #the acceptance rate of uniform proposal sum(diff(X2)!=0)/nrow(X2) #the acceptance rate of normal proposal ## Clearly Uniform, compare to Normal, can better resemble Beta, so the acceptance rate is higher ## plot the results hist(X1) hist(X2) hist(rbeta(5000,2.7,6.3)) ## example2: independent Metropolis-Hastings algorithm, sample from an improper distribution ## p(x) = -|x|+1, where -1

bbricks documentation built on July 8, 2020, 7:29 p.m.