# Simulate some data (i.e. shoot some arrows)
# First we need the number of arrows to shoot
n = 10
# Then we need some true mean distance
mu_true = 1
# Simulate the observed mean distance of the arrows we shot
arrow_mean = rgamma(1, n, n/mu_true)[1]


# Initialize the chain with some starting value
alpha = 1.0
mu = rexp(1, 1/alpha)[1]


# Define the likelihood function on the mean 
function likelihood(mu){
    if(mu < 0.0)
        return 0.0

    return dgamma(arrow_mean, n, n/mu, log=false)
}

# Define the prior function on the mean 
function prior(mu){
    if(mu < 0.0)
        return 0.0

    return dexp(mu, 1/alpha, log=false)
}


# Prepare a file to log our samples
write("iteration","mu","\n",file="archery_MH.log")
write(0,mu,"\n",file="archery_MH.log",append=TRUE)

# Print the initial values to the screen 
print("iteration","mu")
print(0,mu)

# Write the MH algorithm
printgen = 10
reps = 10000
delta = 1

for(rep in 1:reps){
    # Propose a new value of mu
    mu_prime <- mu + runif(n=1,-delta,delta)[1]

    # Compute the acceptance probability
    R = ( likelihood(mu_prime)/likelihood(mu) ) * ( prior(mu_prime)/prior(mu) ) 
    
    # Accept or reject the proposal
    u = runif(1,0,1)[1] 
    if(u < R){
        # Accept the proposal
        mu = mu_prime 
    }

    if ( (rep % printgen) == 0 ) {
        # Write the samples to a file
        write(rep,mu,"\n",file="archery_MH.log",append=TRUE)
        # Print the samples to the screen
        print(rep,mu)
    }
} # end MCMC



q()