Add yearly totals

nwfsc-cb · Jul 10, 2023 · 7a78dc1 · 7a78dc1
1 parent a31adfe
commit 7a78dc1
Showing 1 changed file with 49 additions and 1 deletion.
diff --git a/src/phenomix.cpp b/src/phenomix.cpp
@@ -237,7 +237,7 @@ Type objective_function<Type>::operator() ()
   vector<Type> alpha1(nLevels), alpha2(nLevels);
   vector<Type> lower25(nLevels), upper75(nLevels);
   vector<Type> range(nLevels); // 75th - 25th percentile
-  int i;
+  int i, t;
   int n = y.size();
   vector<Type> log_dens(n);
   vector<Type> pred(n);
@@ -357,6 +357,49 @@ Type objective_function<Type>::operator() ()
     }
   }
 
+  // this is for the cumulative annual predictions
+  vector<Type> year_log_tot(nLevels);
+  vector<Type> year_tot(nLevels);
+  for(i = 0; i < nLevels; i++) {
+    year_log_tot(i) = 0;
+    year_tot(i) = 0;
+    for(t = 1; t < 366; t++) {
+
+        if(asymmetric == 1) {
+          // model is asymmetric, left side smaller / right side bigger
+          if(tail_model==0) {
+            log_dens(i) = ddnorm(t, mu(years(i)-1), sigma1(years(i)-1), sigma2(years(i)-1));
+          }
+          if(tail_model==1) {
+            log_dens(i) = ddt(t, mu(years(i)-1), sigma1(years(i)-1), sigma2(years(i)-1), Type(tdf_1), Type(tdf_2));
+          }
+          if(tail_model==2) {
+            log_dens(i) = ddgnorm(t, mu(years(i)-1), alpha1(years(i)-1), alpha2(years(i)-1), beta_1, beta_2, sigma1(years(i)-1), sigma2(years(i)-1));
+          }
+          year_log_tot(i) = year_log_tot(i) + log_dens(i) + theta(years(i)-1);
+          year_tot(i) = year_tot(i) + exp(log_dens(i) + theta(years(i)-1));
+
+        } else {
+          if(tail_model==0) {
+            // model is symmetric around mu, gaussian tails
+            log_dens(i) = dnorm(t, mu(years(i)-1), sigma1(years(i)-1), true);
+          } else {
+            if(tail_model==1) {
+              // model is symmetric around mu, student-t tails
+              log_dens(i) = dt((t - mu(years(i)-1)) / sigma1(years(i)-1), tdf_1, true) - log(sigma1(years(i)-1));
+            } else {
+              // gnorm, copied from maryclare/gnorm
+              // alpha = sqrt( var * gamma(1/beta) / gamma(3/beta) ), alpha = sigma(1)*beta_ratio(1)
+              log_dens(i) = dgnorm(t, mu(years(i)-1), alpha1(years(i)-1), beta_1);
+            }
+          }
+          year_log_tot(i) = year_log_tot(i) + log_dens(i) + theta(years(i)-1);
+          year_tot(i) = year_tot(i) + exp(log_dens(i) + theta(years(i)-1));
+        }
+    }
+  }
+
+
   // this is the likelihood
   Type s1 = 0;
   Type s2 = 0;
@@ -399,6 +442,11 @@ Type objective_function<Type>::operator() ()
   REPORT(b_mu);
   ADREPORT(b_sig1); // mean of curves by year
   REPORT(b_sig1);
+  ADREPORT(year_tot); // mean of curves by year
+  REPORT(year_tot);
+  ADREPORT(year_log_tot); // mean of curves by year
+  REPORT(year_log_tot);
+
   if(family != 2 && family != 4) {
     ADREPORT(obs_sigma); // obs sd (or phi, NB)
     REPORT(obs_sigma);