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