################################################################################ # # RevBayes Example: Writing your on MCMC algorithm for a simple Binomial model # # authors: Mike May and Sebastian Hoehna # ################################################################################ # Make up some coin flips! # Feel free to change these numbers n <- 100 # the number of flips x <- 63 # the number of heads # Initialize the chain with starting values alpha <- 1 beta <- 1 p <- rbeta(n=1,alpha,beta)[1] # specify the likelihood function function likelihood(p) { if(p < 0 || p > 1) return 0 l = dbinomial(x,p,n,log=false) return l } # specify the prior function function prior(p) { if(p < 0 || p > 1) return 0 pp = dbeta(p,alpha,beta,log=false) return pp } # specify the uniform move function move_uniform( Natural weight) { for (i in 1:weight) { # Propose a new value of p p_prime <- runif(n=1,0.0,1.0)[1] # Compute the acceptance probability R <- ( likelihood(p_prime) / likelihood(p) ) * ( prior(p_prime) / prior(p) ) # Accept or reject the proposal u <- runif(1,0,1)[1] if (u < R){ # Accept the proposal p <- p_prime } else { # Reject the proposal # (we don't have to do anything here) } } } # specify the sliding-window move function move_slide( RealPos delta, Natural weight) { for (i in 1:weight) { # Propose a new value of p p_prime <- p + runif(n=1,-delta,delta)[1] # Compute the acceptance probability R <- ( likelihood(p_prime) / likelihood(p) ) * ( prior(p_prime) / prior(p) ) # Accept or reject the proposal u <- runif(1,0,1)[1] if (u < R){ # Accept the proposal p <- p_prime } else { # Reject the proposal # (we don't have to do anything here) } } } # specify the scaling move function move_scale( RealPos lambda, Natural weight) { for (i in 1:weight) { # Propose a new value of p sf <- exp( lambda * ( runif(n=1,0,1)[1] - 0.5 ) ) p_prime <- p * sf # Compute the acceptance probability R <- ( likelihood(p_prime) / likelihood(p) ) * ( prior(p_prime) / prior(p) ) * sf # Accept or reject the proposal u <- runif(1,0,1)[1] if (u < R){ # Accept the proposal p <- p_prime } else { # Reject the proposal # (we don't have to do anything here) } } } # Prepare a file to log our samples write("iteration","p","\n",file="binomial_MH_moves.log") write(0,p,"\n",file="binomial_MH_moves.log",append=TRUE) # Print the initial values to the screen print("iteration","p") print(0,p) # Write the MH algorithm reps = 10000 printgen = 10 for (rep in 1:reps) { move_uniform(1) move_slide(1,1) move_scale(1,1) if ( (rep % printgen) == 0 ) { # Write the samples to a file write(rep,p,"\n",file="binomial_MH_moves.log",append=TRUE) # Print the samples to the screen print(rep,p) } } # end MCMC q()