Vb model add

.gitignore

@@ -60,9 +60,15 @@ inst/doc
 *.so
 src/rmot.dll
 src/*.tmp
+inst/stan/*.exe
 
 # Bugged directory from snapshot package
 tests/testthat/_snaps
 
+#Testing material
+testing.r
+testingOutput/
+vBtesting.r
+
 # Data prep files
 data-raw/

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"
 
@@ -105,8 +107,8 @@ rmot_power_single_ind <- function(){
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
-       y_bar = NULL,
        y_0_obs = NULL,
+       y_bar = NULL,
        model = "power_single_ind")
 }
 
@@ -122,7 +124,39 @@ rmot_power_multi_ind <- function(){
        obs_index = NULL,
        time = NULL,
        ind_id = NULL,
-       y_bar = NULL,
        y_0_obs = NULL,
+       y_bar = 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,
+       y_bar = 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,
+       y_bar = 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/power_multi_ind.stan

@@ -3,11 +3,11 @@
 functions{
   //DE function for use with Runge-Kutta method
   //pars = [ind_coeff, ind_power]
-  real DE(real y, vector pars){
+  real DE(real y, array[] real pars){
     return pars[1] * pow((y/pars[3]), -pars[2]);
   }
 
-  real rk4_step(real y, vector pars, real interval){
+  real rk4_step(real y, array[] real pars, real interval){
     real k1;
     real k2;
     real k3;
@@ -24,7 +24,7 @@ functions{
     return y_hat;
   }
 
-  real rk4(real y, vector pars, real interval, real step_size){
+  real rk4(real y, array[] real pars, real interval, real step_size){
     int steps;
     real duration;
     real y_hat;
@@ -57,8 +57,8 @@ data {
   int obs_index[n_obs];
   real time[n_obs];
   int ind_id[n_obs];
-  real y_bar;
   real y_0_obs[n_ind];
+  real y_bar;
 }
 
 // The parameters accepted by the model.
@@ -81,7 +81,7 @@ parameters {
 // The model to be estimated.
 model {
   real y_hat[n_obs];
-  vector[3] pars;
+  array[3] real pars;
 
   for(i in 1:n_obs){
     // Initialise the parameters for the observation
@@ -123,7 +123,7 @@ model {
 generated quantities{
   real y_hat[n_obs];
   real Delta_hat[n_obs];
-  vector[3] pars;
+  array[3] real pars;
 
   for(i in 1:n_obs){
     // Initialise the parameters for the observation

inst/stan/power_single_ind.stan

@@ -3,11 +3,11 @@
 functions{
   //DE function for use with Runge-Kutta method
   //pars = [ind_coeff, ind_power]
-  real DE(real y, vector pars){
+  real DE(real y, array[] real pars){
     return pars[1] * pow((y/pars[3]), -pars[2]);
   }
 
-  real rk4_step(real y, vector pars, real interval){
+  real rk4_step(real y, array[] real pars, real interval){
     real k1;
     real k2;
     real k3;
@@ -24,7 +24,7 @@ functions{
     return y_hat;
   }
 
-  real rk4(real y, vector pars, real interval, real step_size){
+  real rk4(real y, array[] real pars, real interval, real step_size){
     int steps;
     real duration;
     real y_hat;
@@ -55,16 +55,16 @@ data {
   real y_obs[n_obs];
   int obs_index[n_obs];
   real time[n_obs];
-  real y_bar;
   real y_0_obs;
+  real y_bar;
 }
 
 // The parameters accepted by the model.
 parameters {
   //Individual level
   real<lower=0> ind_y_0;
-  real<lower=0> ind_coeff;
-  real<lower=0> ind_power;
+  real<lower=0> ind_beta_0;
+  real<lower=0> ind_beta_1;
 
   //Global level
   real<lower=0> global_error_sigma;
@@ -73,10 +73,10 @@ parameters {
 // The model to be estimated.
 model {
   real y_hat[n_obs];
-  vector[3] pars;
+  array[3] real pars;
 
-  pars[1] = ind_coeff;
-  pars[2] = ind_power;
+  pars[1] = ind_beta_0;
+  pars[2] = ind_beta_1;
   pars[3] = y_bar;
 
   for(i in 1:n_obs){
@@ -97,8 +97,8 @@ model {
   //Priors
   //Individual level
   ind_y_0 ~ normal(y_0_obs, global_error_sigma);
-  ind_coeff ~lognormal(0, 1);
-  ind_power ~lognormal(0, 1);
+  ind_beta_0 ~lognormal(0, 1);
+  ind_beta_1 ~lognormal(0, 1);
 
   //Global level
   global_error_sigma ~cauchy(0,5);
@@ -107,10 +107,15 @@ model {
 generated quantities{
   real y_hat[n_obs];
   real Delta_hat[n_obs];
-  vector[3] pars;
+  array[3] real pars;
+  real ind_coeff;
+  real ind_power;
+
+  ind_coeff = ind_beta_0 / pow(y_bar, -ind_beta_1);
+  ind_power = ind_beta_1;
 
-  pars[1] = ind_coeff;
-  pars[2] = ind_power;
+  pars[1] = ind_beta_0;
+  pars[2] = ind_beta_1;
   pars[3] = y_bar;
 
   for(i in 1:n_obs){

inst/stan/vb_multi_ind.stan

@@ -0,0 +1,150 @@
+//Growth function
+functions{
+  //Growth function for use with Runge-Kutta method
+  //pars = (growth_par, max_size)
+  real DE(real y, array[] real pars){ //change number of pars
+    return pars[1] - (pars[2] * (y-pars[3])); //growth function
+  }
+
+  real rk4_step(real y, array[] real 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, array[] real 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];
+  real y_bar;
+}
+
+// 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];
+  array[3] real 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]];
+    pars[3] = y_bar;
+
+    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];
+  array[3] real 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]];
+    pars[3] = y_bar;
+
+    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]);
+    }
+  }
+}