fix random effects

NOAA-FIMS · Jul 26, 2024 · c93e9f7 · c93e9f7
1 parent 4d5cd14
commit c93e9f7
Show file tree

Hide file tree

Showing 4 changed files with 14 additions and 9 deletions.
diff --git a/inst/include/distributions/ar1_lpdf.hpp b/inst/include/distributions/ar1_lpdf.hpp
@@ -32,10 +32,13 @@ struct AR1LPDF : public DensityComponentBase<Type> {
             rho[0] = 1 / (1 + exp(-logit_rho[0])) * 2 - 1;
         }
 
-        sd[0] = exp(log_sd[0]);
-        Rcout << "AR1 sd is: " << sd[0] << std::endl;
-        Rcout << "AR1 rho is:" << rho[0] << std::endl;
+        sd.resize(log_sd.size());
+        for(size_t i=0; i<log_sd.size(); i++){
+            sd[i] = exp(log_sd[i]);
+        }
+
         log_likelihood = -density::SCALE(density::AR1(rho[0]), sd[0])(this->observed_value);
+        this->log_likelihood_vec.resize(1);
         this->log_likelihood_vec[0] = log_likelihood;
 
         return(log_likelihood);

diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp
@@ -93,7 +93,7 @@ class PopulationInterface  : public PopulationInterfaceBase {
         } else {
           pop->u.resize(this->u.size());
           for(int i =0; i < this->u.size(); i++){
-            pop->u[i] = this->u[i];
+            pop->u[i] = this->u[i].value;
             if(this->u[i].is_random_effect){
               model->random_effects.push_back(&(pop)->u[i]);
             }

diff --git a/tests/testthat/test-random-effects.R b/tests/testthat/test-random-effects.R
@@ -70,7 +70,8 @@ Data <- list(
 
 #create a parameter list
 Parameters <- list(
-  p = get_parameter_vector()
+  p = get_parameter_vector(),
+  re = get_random_effects()
 )
 
 obj <- TMB::MakeADFun(Data, Parameters, DLL="ModularTMBExample", trace = TRUE)
@@ -159,7 +160,7 @@ Pop$ages<-ages
 Pop$set_growth(vonB$get_id())
 init.re <- rnorm(nobs)
 Pop$u <- new(VariableVector, init.re, length(length.data))
-Pop$u$set_all_random_effects(TRUE) 
+#Pop$u$set_all_random_effects(TRUE) 
 
 DataLL <- new(NormalLPDF)
 
@@ -173,13 +174,13 @@ DataLL$simulate_flag <- TRUE
 DataLL$set_distribution_links("data", Pop$get_id(), Pop$get_module_name(), "length")
 
 Ar1LL <- new(AR1LPDF)
-Ar1LL$rho <- new(VariableVector, 1.0, 1)
+Ar1LL$rho <- new(VariableVector, .99, 1)
 Ar1LL$log_sd <- new(VariableVector, 1)
 Ar1LL$log_sd[1]$value <- 0 
 Ar1LL$log_sd[1]$estimable <- TRUE
 Ar1LL$input_type <- "re"
 Ar1LL$observed_value <- new(VariableVector, rep(0, nobs), nobs)
-
+Ar1LL$observed_value$set_all_random_effects(TRUE)
 Ar1LL$set_distribution_links("re", Pop$get_id(), Pop$get_module_name(), "u")
 
 

diff --git a/tests/testthat/test-single-nll.R b/tests/testthat/test-single-nll.R
@@ -60,7 +60,8 @@ CreateModel()
 
 #create a data list (data set above)
 Data <- list(
-  y = get_data_vector()
+  y = get_data_vector(),
+  re = get_random_effects_vector()
 )
 
 #create a parameter list