Some more changes getting it ready to send to Josh

NMML · Dec 31, 2014 · c10ec08 · c10ec08
1 parent eb35431
commit c10ec08
Show file tree

Hide file tree

Showing 8 changed files with 193 additions and 163 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,7 +1,8 @@
 # Generated by roxygen2 (4.0.2): do not edit by hand
 
 export(make_ikfdnm_likelihood)
-export(sim_group)
+export(make_ikfdnm_sampler)
+export(sim_data)
 import(Rcpp)
 import(RcppArmadillo)
 useDynLib(kfdnm)
diff --git a/R/r_wrappers.R b/R/r_wrappers.R
@@ -3,9 +3,9 @@
 #' @description This function assimilates the data and desirted model to create a likelihood function that is a function with only one argument, \code{par},
 #' This likelihood function can teither be optimized via \code{\link{optim}} or placed within an MCMC routine.
 #' 
-#' @param dnm_survival A model statement for the survival parameters of the N-mixture portion of the model
-#' @param dnm_recruit A model statement for the recruitment portion of the N-mixture model
-#' @param dnm_det A model statement for the the detection portion of the N-mixture model
+#' @param survival A model statement for the survival parameters of the N-mixture portion of the model
+#' @param recruit A model statement for the recruitment portion of the N-mixture model
+#' @param detection A model statement for the the detection portion of the N-mixture model
 #' @param kf_survival_effects A model statement describing the difference between the known-fate survival and the survival of the N-mixture described in 
 #' \code{dnm_survival}. 
 #' @param fixed_list A list object with named entries: \code{survival}, \code{recruit}, and \code{det}. Each of the entries is a vector of length 
@@ -17,10 +17,8 @@
 #' then time within group. The data must contain the following \code{colnames}: 
 #' \itemize{
 #' \item \code{group} A factor variable indicating the survey group.
-#' \item \code{num_released} The number of known-fate individuals released into the population at the time of the survey
-#' \item \code{num_returns} The number of known-fate individuals observed at each survey time. These are a binomial sample of the \code{num_release} of the 
-#' previous row.
 #' \item \code{n} The observed abundance count at each survey time.
+#' \item \code{R} The number of new KF individuals sampled from the population at each time.
 #' }
 #' 
 #' The \code{kf_survival_effects} model state provides a method for specifying the \emph{difference} in survival between the N-mixtue individuals and the 
@@ -178,81 +176,81 @@ make_ikfdnm_sampler = function(dnm_survival=~1, dnm_recruit=~1, dnm_det=~1, kf_s
 }
 
 
-#' @title Perform MCMC sampling for integrated known-fate, dynamic N-mixture model
-#' @description This function used a KFDNM likelihood function, as well as, a user defined prior distribution
-#' for the parameters and uses it to perform MCMC sampling of the parameter vector. 
-#' 
-#' @param kfdnm_lik A likelihood function created with the \code{\link{make_ikfdnm_likelihood}} function
-#' @param prior_dist A user defined function that evaluates to the log-prior distribution of the parameter vector
-#' @param sample_N Logical indicating that the abundance process should be sampled.
-#' @param sample_SG Logical indicating wheather the survival (S) and recruitement process should be sampled.
-#' @param burn The number of iterations for burnin of the MCMC algorithm
-#' @param iter The number of MCMC iterations
-#' @param block Block size for adaptation of proposal variance. 
-#' @param par Initial parameter value. If not provided a optimization is executed first to obtain 
-#' realistic initial values. 
-#' @return A list containg the MCMC sample of parameters, abundance, and dynamic abundance components.   
-
-ikfdnm_mcmc = function(dnm_survival=~1, dnm_recruit=~1, dnm_det=~1, kf_survival_effects=NULL, 
-                       fixed_list=NULL, data, N_max, ln_prior, sample_abundance=FALSE,
-                       burn, iter, block, par){
-  n2ll=make_ikfdnm_likelihood(
-    dnm_survival=dnm_survival, dnm_recruit=dnm_recruit, dnm_det=dnm_det, 
-    kf_survival_effects=kf_survival_effects, fixed_list=fixed_list, 
-    data=data, N_max=N_max
-    )
-  NSG_sampler = make_ikfdnm_sampler(
-    dnm_survival=dnm_survival, dnm_recruit=dnm_recruit, dnm_det=dnm_det, 
-    kf_survival_effects=kf_survival_effects, fixed_list=fixed_list, 
-    data=data, N_max=N_max
-    )
-  d = attr(kfdnm_lik, "npar")
-  ln_post = function(par){-n2ll(par)/2 + prior_dist(par)}
-  if(missing(par)){
-    message("Computing starting values...")
-    par_curr=optim(rep(0,d), ln_post, control=list(fnscale=-0.5))$par
-  } else{par_curr=par}
-  ln_post_curr = ln_post(par_curr)
-  tune = 2.4^2/d
-  Sigma = 0.01*diag(d)
-  r_opt = 0.234
-  r_hist = rep(0,burn+iter)
-  par_store = matrix(NA, burn+iter, d)
-  ### Update par
-  message("Begininng MCMC...")
-  flush.console()
-  pb = txtProgressBar(min = 1, max = burn+iter, style = 3)
-  for(i in 1:(burn+iter)){
-    par_prop = par_curr + sqrt(tune)*crossprod(chol(Sigma),rnorm(d,0,1))
-    ln_post_prop = ln_post(par_prop)
-    mh = exp(ln_post_prop - ln_post_curr)
-    if(runif(1,0,1)<=mh){
-      par_curr = par_prop
-      r_hist[i] = 1
-      ln_post_curr =ln_post_prop
-    }
-    par_store[i,]=par_curr
-    ### Adapt par update proposal
-    if(i%%block==0){
-      r_hat = sum(r_hist[(i-block+1):i])/block
-      Sigma_hat = var(par_store[(i-block+1):i,])
-      tune=exp(log(tune) + 2*((i/block)^-0.8)*(r_hat-r_opt))
-      Sigma = Sigma + ((i/block)^-0.8)*(Sigma_hat-Sigma)
-    }
-    ### Update abundance
-    if(sample_abundance){
-      abund=NSG_sampler(par_curr)
-    }
-
-    setTxtProgressBar(pb,i)
-  }
-  close(pb)
-  message("Processing posterior sample for output...")
-  return(list(
-    par_sample=par_store,
-    N=NULL,
-    S=NULL,
-    G=NULL,
-    accept=r_hist
-    ))
-}
+# #' @title Perform MCMC sampling for integrated known-fate, dynamic N-mixture model
+# #' @description This function used a KFDNM likelihood function, as well as, a user defined prior distribution
+# #' for the parameters and uses it to perform MCMC sampling of the parameter vector. 
+# #' 
+# #' @param kfdnm_lik A likelihood function created with the \code{\link{make_ikfdnm_likelihood}} function
+# #' @param prior_dist A user defined function that evaluates to the log-prior distribution of the parameter vector
+# #' @param sample_N Logical indicating that the abundance process should be sampled.
+# #' @param sample_SG Logical indicating wheather the survival (S) and recruitement process should be sampled.
+# #' @param burn The number of iterations for burnin of the MCMC algorithm
+# #' @param iter The number of MCMC iterations
+# #' @param block Block size for adaptation of proposal variance. 
+# #' @param par Initial parameter value. If not provided a optimization is executed first to obtain 
+# #' realistic initial values. 
+# #' @return A list containg the MCMC sample of parameters, abundance, and dynamic abundance components.   
+# 
+# ikfdnm_mcmc = function(dnm_survival=~1, dnm_recruit=~1, dnm_det=~1, kf_survival_effects=NULL, 
+#                        fixed_list=NULL, data, N_max, ln_prior, sample_abundance=FALSE,
+#                        burn, iter, block, par){
+#   n2ll=make_ikfdnm_likelihood(
+#     dnm_survival=dnm_survival, dnm_recruit=dnm_recruit, dnm_det=dnm_det, 
+#     kf_survival_effects=kf_survival_effects, fixed_list=fixed_list, 
+#     data=data, N_max=N_max
+#     )
+#   NSG_sampler = make_ikfdnm_sampler(
+#     dnm_survival=dnm_survival, dnm_recruit=dnm_recruit, dnm_det=dnm_det, 
+#     kf_survival_effects=kf_survival_effects, fixed_list=fixed_list, 
+#     data=data, N_max=N_max
+#     )
+#   d = attr(kfdnm_lik, "npar")
+#   ln_post = function(par){-n2ll(par)/2 + prior_dist(par)}
+#   if(missing(par)){
+#     message("Computing starting values...")
+#     par_curr=optim(rep(0,d), ln_post, control=list(fnscale=-0.5))$par
+#   } else{par_curr=par}
+#   ln_post_curr = ln_post(par_curr)
+#   tune = 2.4^2/d
+#   Sigma = 0.01*diag(d)
+#   r_opt = 0.234
+#   r_hist = rep(0,burn+iter)
+#   par_store = matrix(NA, burn+iter, d)
+#   ### Update par
+#   message("Begininng MCMC...")
+#   flush.console()
+#   pb = txtProgressBar(min = 1, max = burn+iter, style = 3)
+#   for(i in 1:(burn+iter)){
+#     par_prop = par_curr + sqrt(tune)*crossprod(chol(Sigma),rnorm(d,0,1))
+#     ln_post_prop = ln_post(par_prop)
+#     mh = exp(ln_post_prop - ln_post_curr)
+#     if(runif(1,0,1)<=mh){
+#       par_curr = par_prop
+#       r_hist[i] = 1
+#       ln_post_curr =ln_post_prop
+#     }
+#     par_store[i,]=par_curr
+#     ### Adapt par update proposal
+#     if(i%%block==0){
+#       r_hat = sum(r_hist[(i-block+1):i])/block
+#       Sigma_hat = var(par_store[(i-block+1):i,])
+#       tune=exp(log(tune) + 2*((i/block)^-0.8)*(r_hat-r_opt))
+#       Sigma = Sigma + ((i/block)^-0.8)*(Sigma_hat-Sigma)
+#     }
+#     ### Update abundance
+#     if(sample_abundance){
+#       abund=NSG_sampler(par_curr)
+#     }
+#     
+#     setTxtProgressBar(pb,i)
+#   }
+#   close(pb)
+#   message("Processing posterior sample for output...")
+#   return(list(
+#     par_sample=par_store,
+#     N=NULL,
+#     S=NULL,
+#     G=NULL,
+#     accept=r_hist
+#     ))
+# }
diff --git a/R/sim_group.R b/R/sim_group.R
@@ -15,34 +15,46 @@
 #' 
 #' @export
 #' 
-sim_group = function(num_kf, num_years, num_surveys, recruit_rate, init_rate, 
+sim_data = function(num_kf, num_years, num_surveys, recruit_rate, init_rate, 
                     survival_dnm, survival_kf, perfect_survey_rate, detection){
   out = data.frame(year=rep(1:num_years, each=num_surveys))
   out$survey=rep(1:num_surveys, num_years)
   out$perfect_survey=rbinom(num_years*num_surveys, 1, perfect_survey_rate)
   out$Y = 0
   out$R = 0
-  out$R[1]=num_kf
   N = rep(NA,num_years*num_surveys) 
   N[1] = rpois(1,init_rate)+num_kf
   S = rep(NA,num_years*num_surveys)
   G = rep(0,num_years*num_surveys)
   for(r in 2:(num_years*num_surveys)){
-    out$Y[r]=rbinom(1, out$Y[r-1]+out$R[r-1], survival_kf)
     S[r]=rbinom(1, N[r-1]-out$R[r-1], survival_dnm)
     if(out$survey[r]==1) {
       G[r]=rpois(1,recruit_rate)
     }
     N[r] = S[r]+ G[r]
-    if(out$survey[r]==1){
-      out$R[r]= min(num_kf-out$Y[r], N[r])
-    }
+#     if(out$survey[r]==1){
+#       out$R[r]= min(num_kf-out$Y[r], N[r])
+#     }
   }
   det=rep(detection, num_years*num_surveys)
   det=ifelse(out$perfect_survey==1, 1, det)
   out$S=S
   out$G=G
   out$N=N
   out$n = rbinom(num_years*num_surveys, N-out$R, det)
-  return(out)
+
+  out_kf=NULL
+  for(i in 1:num_years){
+    Y = matrix(nrow=num_surveys, ncol=num_kf)
+    Y[1,] = 1
+    for(j in 2:num_surveys){
+      Y[j,] = rbinom(n=num_kf, size=Y[j-1,], prob=survival_kf)
+    }
+    CH=as.vector(Y)
+    id = rep(paste(i, 1:num_kf, sep="-"), each=num_surveys)
+    survey=rep(1:num_surveys, num_kf)
+    out_kf = rbind(out_kf, data.frame(id, year=i, survey, CH))
+  }
+
+  return(list(dnm_data=out, kf_data=out_kf))
 }
diff --git a/demo/sim_data_example.R b/demo/sim_data_example.R
@@ -3,40 +3,44 @@ require(kfdnm)
 ###
 ### Create data to mimic wolf pack example
 ###
-
-data=NULL
-num_groups=20
+num_groups=5
+dnm_data = NULL
+kf_data = NULL
 for(i in 1:num_groups){
-  data = rbind(data,
-              sim_group(
-                num_kf = 3, 
-                num_years = 5, 
-                num_surveys = 10, 
-                recruit_rate = 4, 
-                init_rate = 6, 
-                survival_dnm = 0.95, 
-                survival_kf = 0.95, 
-                perfect_survey_rate = 0.1, 
-                detection = 0.5
-              )
+  data_list = sim_data(
+    num_kf = 3, 
+    num_years = 5, 
+    num_surveys = 10, 
+    recruit_rate = 4, 
+    init_rate = 6, 
+    survival_dnm = 0.95, 
+    survival_kf = 0.95, 
+    perfect_survey_rate = 0.1, 
+    detection = 0.5
   )
+  data_list$kf_data$id = paste(i, data_list$kf_data$id, sep="-")
+  data_list$dnm_data$group=i
+  dnm_data = rbind(dnm_data, data_list$dnm_data) 
+  kf_data = rbind(kf_data, data_list$kf_data)
 }
-data$group = rep(1:num_groups, each=5*10)
+kf_data$id = factor(kf_data$id)
+dnm_data$group = factor(dnm_data$group)
+
 
 ###
 ### Create list of fixed parameters
 ###
 
-fixed_list = list(survival=rep(NA, nrow(data)), 
-                  recruit=ifelse(data$year!=1 & data$survey==1, NA, 0),
-                  det=ifelse(data$perfect==1, 1, NA)
-                  )
+fixed_list = list(
+  recruit=ifelse(dnm_data$year!=1 & dnm_data$survey==1, NA, 0),
+  detection=ifelse(dnm_data$perfect==1, 1, NA)
+)
 
 ###
 ### Make likelihood function
 ### 
-n2ll_kfdnm = make_ikfdnm_likelihood(dnm_survival=~1, dnm_recruit=~1, dnm_det=~1, kf_survival_effects=NULL, 
-                                  fixed_list=fixed_list, data=data, N_max=50)
+n2ll_kfdnm = make_ikfdnm_likelihood(survival=~1, recruit=~1, detection=~1, kf_survival_effects=NULL, 
+                                    fixed_list=fixed_list, kf_data=kf_data, dnm_data=dnm_data, N_max=50)
 
 ###
 ### Optimize and obtain estimates and variance-covariance matrix
@@ -56,8 +60,8 @@ se = sqrt(diag(2*solve(mle$hessian)))
 ###
 
 # omega
-cat("Annual survival: \n")
-cat(plogis(par[1])^10, "(", plogis(par[1]-2*se[1])^10, ",",plogis(par[1]+2*se[1])^10, ")\n")
+cat("Survival: \n")
+cat(plogis(par[1]), "(", plogis(par[1]-2*se[1]), ",",plogis(par[1]+2*se[1]), ")\n")
 
 # gamma
 cat("Recruitment rate:")

diff --git a/man/ikfdnm_mcmc.Rd b/man/ikfdnm_mcmc.Rd