MCMC sampling

Introduction

If the new solution costs less than the old solution then it is guaranteed to be accepted and a poorer solution is unlikely to be accepted. it is to calculate the integral, proportion, and area.

Data simulation

# Plot the exponential distribution:
x <- seq(from = 0, to = 10, length=20)
y <- exp(-x)  #y is only density function
plot(x,y)
lines(x,y)

mean(x )

## [1] 5

Set up a target function

# set up a target function, we want to sample from in this case, is exponent function 
target <- function(x){
  if (x < 0) {
    return(0)
  } else {
    return(exp(-x))   #target function
  }
}

target(3)

## [1] 0.04978707

target(-1)

## [1] 0

Sampling from a proposal distribution proportional =normal to the target =exponent

# program a Metropolis–Hastings scheme to sample 
x <- rep(0, 1000)
x[1] <- 30     # this is just a starting value, which I’ve set arbitrarily to 3

for (i in 2:1000) {
  # i=2
  currentx <- x[i - 1]
  proposedx <-   rnorm(1, mean = currentx, sd = 1)
  A <- target(proposedx) / target(currentx)
  if (runif(1) < A) {
    x[i] <- proposedx
  } else {
    x[i] <- currentx
  }
}

plot(x)

Wrap this up in a function

easyMCMC <- function(niter, startval, proposalsd) {
  x <- rep(0, niter)
  x[1] <- startval
  for (i in 2:niter) {
    currentx <- x[i - 1]  #status 1
    proposedx <-  rnorm(1, mean = currentx, sd = proposalsd)  #status 2
    #transform matrix, that is the matter, to sample target sampling
    
    A <- target(proposedx) / target(currentx)  
    if (runif(1) < A) {
      x[i] <- proposedx
    } else {
      x[i] <- currentx
    }
  }
  return(x)
}

Run the MCMC scheme 3 times with different start values

It will be stable even if the start value is large.

# 
par(mfrow = c(1, 1))
z1 <- easyMCMC(1000,50,1)
z2 <- easyMCMC(1000,5,1)
z3 <- easyMCMC(1000,0.5,1)
plot(z1, type = "l",col = "red")
lines(z2, col = "black")
lines(z3, col = "green")

# posterior mean
mean(z1)

## [1] 5.572372

How to calculate pie?