################################################################################ # # 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 } # Prepare a file to log our samples write("iteration","p","\n",file="binomial_MH.log") write(0,p,"\n",file="binomial_MH.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){ # 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) } if ( (rep % printgen) == 0 ) { # Write the samples to a file write(rep,p,"\n",file="binomial_MH.log",append=TRUE) # Print the samples to the screen print(rep,p) } } # end MCMC q()