Vb model add (#18)

Added von Bertalanffy model, testing suite.

* Added code for von Bertalanffy model. Currently has problems with estimating a max size smaller than the starting size.

* Updated von Butterfly test, everything runs locally so here goes.

R/rmot_models.R

@@ -19,7 +19,9 @@ rmot_model <- function(model=NULL){
                    canham_single_ind = rmot_canham_single_ind(),
                    canham_multi_ind = rmot_canham_multi_ind(),
                    power_single_ind = rmot_power_single_ind(),
-                   power_multi_ind = rmot_power_multi_ind())
+                   power_multi_ind = rmot_power_multi_ind(),
+                   vb_single_ind = rmot_vb_single_ind(),
+                   vb_multi_ind = rmot_vb_multi_ind())
 
   class(output) <- "rmot_object"
 
@@ -126,3 +128,33 @@ rmot_power_multi_ind <- function(){
        y_0_obs = NULL,
        model = "power_multi_ind")
 }
+
+#' Data configuration template for von Bertalanffy growth single individual model
+#' @keywords internal
+#' @noRd
+
+rmot_vb_single_ind <- function(){
+  list(step_size = NULL,
+       n_obs = NULL,
+       y_obs = NULL,
+       obs_index = NULL,
+       time = NULL,
+       y_0_obs = NULL,
+       model = "vb_single_ind")
+}
+
+#' Data configuration template for von Bertalanffy growth single species model
+#' @keywords internal
+#' @noRd
+
+rmot_vb_multi_ind <- function(){
+  list(step_size = NULL,
+       n_obs = NULL,
+       n_ind = NULL,
+       y_obs = NULL,
+       obs_index = NULL,
+       time = NULL,
+       ind_id = NULL,
+       y_0_obs = NULL,
+       model = "vb_multi_ind")
+}

R/rmot_run.R

@@ -22,7 +22,9 @@ rmot_run <- function(model_template, ...) {
          canham_single_ind = rstan::sampling(stanmodels$canham_single_ind, data = model_template, ...),
          canham_multi_ind = rstan::sampling(stanmodels$canham_multi_ind, data = model_template, ...),
          power_single_ind = rstan::sampling(stanmodels$power_single_ind, data = model_template, ...),
-         power_multi_ind = rstan::sampling(stanmodels$power_multi_ind, data = model_template, ...))
+         power_multi_ind = rstan::sampling(stanmodels$power_multi_ind, data = model_template, ...),
+         vb_single_ind = rstan::sampling(stanmodels$vb_single_ind, data = model_template, ...),
+         vb_multi_ind = rstan::sampling(stanmodels$vb_multi_ind, data = model_template, ...))
 
   return(out)
 }

R/stanmodels.R

@@ -1,7 +1,7 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "linear", "power_multi_ind", "power_single_ind")
+stanmodels <- c("canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "linear", "power_multi_ind", "power_single_ind", "vb_multi_ind", "vb_single_ind")
 
 # load each stan module
 Rcpp::loadModule("stan_fit4canham_multi_ind_mod", what = TRUE)
@@ -11,6 +11,8 @@ Rcpp::loadModule("stan_fit4constant_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4linear_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4power_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4power_single_ind_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4vb_multi_ind_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4vb_single_ind_mod", what = TRUE)
 
 # instantiate each stanmodel object
 stanmodels <- sapply(stanmodels, function(model_name) {

inst/stan/vb_multi_ind.stan

@@ -0,0 +1,147 @@
+//Growth function
+functions{
+  //Growth function for use with Runge-Kutta method
+  //pars = (growth_par, max_size)
+  real DE(real y, vector pars){ //change number of pars
+    return pars[1] * (pars[2] - y); //growth function
+  }
+
+  real rk4_step(real y, vector pars, real interval){
+    real k1;
+    real k2;
+    real k3;
+    real k4;
+    real y_hat;
+
+    k1 = DE(y, pars);
+    k2 = DE(y+interval*k1/2.0, pars);
+    k3 = DE(y+interval*k2/2.0, pars);
+    k4 = DE(y+interval*k3, pars);
+
+    y_hat = y + (1.0/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4) * interval;
+
+    return y_hat;
+  }
+
+  real rk4(real y, vector pars, real interval, real step_size){
+    int steps;
+    real duration;
+    real y_hat;
+    real step_size_temp;
+
+    duration = 0;
+    y_hat = y;
+
+    while(duration < interval){
+      //Determine the relevant step size
+      step_size_temp = min([step_size, interval-duration]);
+
+      //Get next size estimate
+      y_hat = rk4_step(y_hat, pars, step_size_temp);
+
+      //Increment observed duration
+      duration = duration + step_size_temp;
+    }
+
+    return y_hat;
+  }
+}
+
+// Data structure
+data {
+  real step_size;
+  int n_obs;
+  int n_ind;
+  real y_obs[n_obs];
+  int obs_index[n_obs];
+  real time[n_obs];
+  int ind_id[n_obs];
+  real y_0_obs[n_ind];
+}
+
+// The parameters accepted by the model.
+parameters {
+  //Individual level
+  real<lower=0> ind_y_0[n_ind];
+  real<lower=0> ind_growth_par[n_ind];
+  real<lower=0> ind_max_size[n_ind];
+
+  //Species level
+  real species_growth_par_mean;
+  real<lower=0> species_growth_par_sd;
+  real species_max_size_mean;
+  real<lower=0> species_max_size_sd;
+
+  //Global level
+  real<lower=0> global_error_sigma;
+}
+
+// The model to be estimated.
+model {
+  real y_hat[n_obs];
+  vector[2] pars;
+
+  for(i in 1:n_obs){
+    // Initialise the parameters for the observation
+    pars[1] = ind_growth_par[ind_id[i]];
+    pars[2] = ind_max_size[ind_id[i]];
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0[ind_id[i]];
+    }
+
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        //Estimate next size
+        y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      }
+    }
+  }
+
+  //Likelihood
+  y_obs ~normal(y_hat, global_error_sigma);
+
+  //Priors
+  //Individual level
+  ind_y_0 ~normal(y_0_obs, global_error_sigma);
+  ind_growth_par ~lognormal(species_growth_par_mean, species_growth_par_sd);
+  ind_max_size ~lognormal(species_max_size_mean, species_max_size_sd);
+
+  //Species level
+  species_growth_par_mean ~normal(0, 2);
+  species_growth_par_sd ~cauchy(0, 2);
+  species_max_size_mean ~normal(max(log(y_obs)), 2);
+  species_max_size_sd ~cauchy(0, 2);
+
+  //Global level
+  global_error_sigma ~cauchy(0, 5);
+}
+
+generated quantities{
+  real y_hat[n_obs];
+  real Delta_hat[n_obs];
+  vector[2] pars;
+
+  for(i in 1:n_obs){
+    // Initialise the parameters for the observation
+    pars[1] = ind_growth_par[ind_id[i]];
+    pars[2] = ind_max_size[ind_id[i]];
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0[ind_id[i]];
+    }
+
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        //Estimate next size
+        y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+        Delta_hat[i] = y_hat[i+1] - y_hat[i];
+
+      } else {
+        Delta_hat[i] = DE(y_hat[i], pars)*(time[i] - time[i-1]);
+      }
+    } else {
+      Delta_hat[i] = DE(y_hat[i], pars)*(time[i] - time[i-1]);
+    }
+  }
+}

inst/stan/vb_single_ind.stan

@@ -0,0 +1,127 @@
+//Growth function
+functions{
+  //Growth function for use with Runge-Kutta method
+  //pars = (growth_par, max_size)
+  real DE(real y, vector pars){ //change number of pars
+    return pars[1] * (pars[2] - y); //growth function
+  }
+
+  real rk4_step(real y, vector pars, real interval){
+    real k1;
+    real k2;
+    real k3;
+    real k4;
+    real y_hat;
+
+    k1 = DE(y, pars);
+    k2 = DE(y+interval*k1/2.0, pars);
+    k3 = DE(y+interval*k2/2.0, pars);
+    k4 = DE(y+interval*k3, pars);
+
+    y_hat = y + (1.0/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4) * interval;
+
+    return y_hat;
+  }
+
+  real rk4(real y, vector pars, real interval, real step_size){
+    int steps;
+    real duration;
+    real y_hat;
+    real step_size_temp;
+
+    duration = 0;
+    y_hat = y;
+
+    while(duration < interval){
+      //Determine the relevant step size
+      step_size_temp = min([step_size, interval-duration]);
+
+      //Get next size estimate
+      y_hat = rk4_step(y_hat, pars, step_size_temp);
+
+      //Increment observed duration
+      duration = duration + step_size_temp;
+    }
+
+    return y_hat;
+  }
+}
+
+// Data structure
+data {
+  real step_size;
+  int n_obs;
+  real y_obs[n_obs];
+  int obs_index[n_obs];
+  real time[n_obs];
+  real y_0_obs;
+}
+
+// The parameters accepted by the model.
+parameters {
+  //Individual level
+  real<lower=0> ind_y_0;
+  real<lower=0> ind_growth_par;
+  real<lower=0> ind_max_size;
+
+  //Global level
+  real<lower=0> global_error_sigma;
+}
+
+// The model to be estimated.
+model {
+  real y_hat[n_obs];
+  vector[2] pars;
+
+  pars[1] = ind_growth_par;
+  pars[2] = ind_max_size;
+
+  for(i in 1:n_obs){
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0;
+    }
+
+    if(i < n_obs){
+      //Estimate growth rate
+      y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+    }
+  }
+
+  //Likelihood
+  y_obs ~ normal(y_hat, global_error_sigma);
+
+  //Priors
+  //Individual level
+  ind_y_0 ~ normal(y_0_obs, global_error_sigma);
+  ind_growth_par ~lognormal(0, 1);
+  ind_max_size ~lognormal(max(log(y_obs)), 1); //Take max obs. size as average value
+
+  //Global level
+  global_error_sigma ~cauchy(1,5);
+}
+
+generated quantities{
+  real y_hat[n_obs];
+  real Delta_hat[n_obs];
+  vector[2] pars;
+
+  pars[1] = ind_growth_par;
+  pars[2] = ind_max_size;
+
+  for(i in 1:n_obs){
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0;
+    }
+
+    if(i < n_obs){
+      //Estimate next size
+      y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      Delta_hat[i] = y_hat[i+1] - y_hat[i];
+
+    } else {
+      Delta_hat[i] = DE(y_hat[i], pars)*(time[i] - time[i-1]);
+    }
+  }
+}

src/RcppExports.cpp

@@ -19,6 +19,8 @@ RcppExport SEXP _rcpp_module_boot_stan_fit4constant_single_ind_mod();
 RcppExport SEXP _rcpp_module_boot_stan_fit4linear_mod();
 RcppExport SEXP _rcpp_module_boot_stan_fit4power_multi_ind_mod();
 RcppExport SEXP _rcpp_module_boot_stan_fit4power_single_ind_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4vb_multi_ind_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4vb_single_ind_mod();
 
 static const R_CallMethodDef CallEntries[] = {
     {"_rcpp_module_boot_stan_fit4canham_multi_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4canham_multi_ind_mod, 0},
@@ -28,6 +30,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_rcpp_module_boot_stan_fit4linear_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4linear_mod, 0},
     {"_rcpp_module_boot_stan_fit4power_multi_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4power_multi_ind_mod, 0},
     {"_rcpp_module_boot_stan_fit4power_single_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4power_single_ind_mod, 0},
+    {"_rcpp_module_boot_stan_fit4vb_multi_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4vb_multi_ind_mod, 0},
+    {"_rcpp_module_boot_stan_fit4vb_single_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4vb_single_ind_mod, 0},
     {NULL, NULL, 0}
 };