######
# Code for running MCMC algorithm for radio-tagging data
# Uses MH random walk updates. Created Sept 08 (RK)
######

# Define the function: with input parameters:
# nt = number of iterations
# nburn  = burn-in

radiocode <- function(nt,nburn){

# Define the parameter values:
# ni = number of release years
# nj = number of recapture years

    ni = 1
    nj = 6

# Read in the data:

data <- matrix(c(3,2,3,3,2,1,22),nrow=ni,byrow=T)

# Read in the priors:

# Beta priors on the parameters

    alpha <- 1
    beta <- 1

# Parameters for MH updates (standard deviation for Normal random walk):

    delta <- 0.5

# Set initial parameter value for phi:

    param <- 0.9

# sample is an array in which we put the sample from the posterior distribution.

    sample <- array(0,dim=c(nt,1))

    sample <- data.frame(sample)

# Label the column of the array "sample"

    names(sample) <- c("phi")

# Calculate log(likelihood) for initial state using a separate function "likhood":

    likhood <- calclikhood(ni, nj, data, param)

# Set up a vector for parameter values and associated log(likelihood):

    output <- dim(2)

# MCMC updates - MH algorithm:

# Cycle through each iteration:

    for (t in 1:nt){

# Update the parameters in the model using function "updateparam":

        output <- updateparam(param, ni, nj, data, likhood, alpha, beta, delta)

# Set parameter values and log(likelihood) value of current state to be the output from
# the MH step:

        param <- output[1]
        likhood <- output[2]

# Record the set of parameter values:

        sample[t,1] <- param

    }

# Calculate the mean and standard deviation of the parameters
# following pre-specified burn-in:

mn <- mean(sample[(nburn+1):nt,1])
std <- sqrt(var(sample[(nburn+1):nt,1]))

# Output the posterior mean and standard deviation of the parameter
# following burn-in to the screen:

cat("Posterior summary estimates for parameter:  ", "\n")
cat("\n")
cat("mean  (SD)", "\n")
cat("phi: ", "\n")
cat(mn, "   (", std, ")", "\n")

# Output the sample from the posterior distribution:

sample

}


######
# This function calculates the log likelihood of capture-recapture data
######

calclikhood <- function(ni, nj, data, param){

# First we set the survival and recapture probs:

# Set up the set of cell probabilities (q), and initially set them to be all equal to zero:

    q <- array(0,dim=c(ni,nj+1))

# Set the survival rate to be the current parameter value

    phi <- param

# Set initial value of the log(likelihood) to be zero

    likhood <- 0

# Calculate the Multinomial cell probabilities and corresponding contribution to the
# log(likelihood):

    for (i in 1:ni){
        for (j in i:nj) {

            q[i,j] <- (1-phi)*(phi^(j-i))
            likhood <- likhood + data[i,j]*log(q[i,j])

        }

# Probability of an animal never being seen again

        q[i,nj+1] <- 1 - sum(q[i,i:nj])
        likhood <- likhood + data[i,nj+1]*log(q[i,nj+1])

    }

# Output the log(likelihood) value:

likhood

}

######
# Function for updating the parameters values:
######

updateparam <- function(param, ni, nj, data, likhood, alpha, beta, delta){

# Propose to update parameter phi using
# random walk MH with Normal proposal density:

# Keep a record of the current parameter value being updated

    oldparam <- param

# Propose a new value for the parameter using a random walk with
# Normal proposal density

    param <- rnorm(1, param, delta)

# Automatically reject any moves that are propose values outside [0,1]

    if (param >= 0 & param <= 1) {

# Calculate the log(acceptance probability):

# Calculate the new likelihood value for the proposed move:
# Calculate the numerator (num) and denominator (den) in turn:

        newlikhood <- calclikhood(ni, nj, data, param)

# Include the likelihood term in the acceptance probability

        num <- newlikhood
        den <- likhood

# Include the likelihood and prior (Beta) terms in the acceptance probability

        num <- newlikhood + log(dbeta(param,alpha,beta))
        den <- likhood + log(dbeta(oldparam,alpha,beta))

# All other prior terms (for other parameters) cancel in the acceptance probability.

# Acceptance probability of MH step:

       A <- min(1,exp(num-den))

    }

    else {

        A <- 0

    }

# To do the accept/reject step of the algorithm:
# Simulate a random number in [0,1]:

    u <- runif(1)

# Accept the move with probability A:

    if (u <= A) {

# Accept the proposed move:
# Update the log(likelihood) value:

        likhood <- newlikhood

    }

    else {

# Reject proposed move so parameter stays at current value:

        param <- oldparam

    }


# Set the values to be outputted from the function to be the
# set of parameter values and log(likelihood) value:

        output <- c(param, likhood)

# Output the parameter values:

    c(output,A)
}