- [[Markov chain]]: memory-less; to get the next sample, we only need to consider the current sample - [[Monte Carlo]]: the algorithm randomly generates samples that we can then analyze - [[Python]] # Idea The Metropolis Hastings MCMC algorithm relies on a simple trick to determine whether to accept a proposal. ```python probability_proposal = 0.5 probability_current = 0.5 probability_accept = probability_proposal / probability_current accept = np.random.rand() < p_accept # randomly determine whether to accept if accept: current_position = proposed # memoryless Markov chain! ``` If `probability_proposal` is equal or greater than `probability_current`, `probability_accept` is always equal or greater than 1, and thus, will always be accepted (i.e., `accept == 1` and `current_position` will be replaced by the `proposed` value). But if `probability_proposal` is less than `probability_current`, `probability_accept` will be a value between 0 and 1. So if `probability_proposal` is half the size of `probability_current`, `probability_accept == 0.5`, and `accept` has 0.5 probability of taking on the value 1 (accept). Function to implement test the above algorithm: ```python def decide(p_proposal, p_current, n): p_accept = p_proposal / p_current print(f"p_accept: {p_accept}") prop = np.mean([np.random.rand() < p_accept for i in range(n)]) print(f"proportion accepted: {prop}") decide(0.01, 1.0, 5000) ``` ## Basic Metropolis Hasting algorithm ```r n_samples <- 1000 p <- rep(NA , n_samples) # posterior samples p[1] <- 0.5 # initial value H <- 6 # heads T <- 3 # tails N <- H + T # total flips for (i in 2:n_samples) { # generate new proposal p_new <- rnorm(n = 1, mean = p[i - 1], sd = 0.1) # p_new <- runif(n = 1, 0, 1) # alternatively, proposal generated from uniform distribution # ensure proposal p has to be within 0 and 1 p_new <- min(c(p_new, 1)) p_new <- max(c(p_new, 0)) # MCMC algorithm # calculate probability of data given preivous step/proposal and current step proposal q0 <- dbinom(H, N, p[i - 1]) # previous step q1 <- dbinom(H, N, p_new) # current step # should we reject or acccept q1 p[i] <- ifelse(runif(1) < q1 / q0, p_new, p[i-1]) } ``` # References - https://twiecki.io/blog/2015/11/10/mcmc-sampling/ - https://wiseodd.github.io/techblog/2015/10/17/metropolis-hastings/ - https://stackoverflow.com/questions/54853017/why-is-my-python-implementation-of-metropolis-algorithm-mcmc-so-slow