################################################################################
#
# 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()