update

DESCRIPTION

@@ -38,7 +38,7 @@ Encoding: UTF-8
 LazyData: true
 URL: https://atsa-es.github.io/atsar/
 BugReports: https://github.com/atsa-es/atsar/issues
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 SystemRequirements: GNU make
 Biarch: true
 VignetteBuilder: knitr

R/fit_stan.R

@@ -186,6 +186,15 @@ fit_stan <- function(y, x = NA, model_name = NA,
     n_pos <- length(row_indx_pos)
     y <- y[which(!is.na(y))]
 
+    est_A <- rep(1, M)
+    for(i in 1:max(marss$states)) {
+      indx <- which(marss$states==i)
+      est_A[indx[1]] <- 0
+    }
+    est_A <- which(est_A > 0)
+    est_A <- c(est_A, 0)
+    n_A <- length(est_A) - 1
+
     data = list("N"=N,"M"=M, "y"=y,
                      "states"=marss$states, "S" = max(marss$states), "obsVariances"=marss$obsVariances,
                      "n_obsvar" = max(marss$obsVariances), "proVariances" = proVariances,
@@ -194,6 +203,8 @@ fit_stan <- function(y, x = NA, model_name = NA,
                      "n_pos" = n_pos,
                      "col_indx_pos" = col_indx_pos,
                      "row_indx_pos" = row_indx_pos,
+                     "est_A" = est_A,
+                     "n_A" = n_A,
                      "y_int"=round(y),
                      "family"=1,
                      "est_trend" = as.numeric(est_trend),
@@ -202,6 +213,7 @@ fit_stan <- function(y, x = NA, model_name = NA,
     pars = c("pred", "log_lik","sigma_process","sigma_obs","x0")
     if(marss$est_B) pars = c(pars, "B")
     if(marss$est_trend) pars = c(pars, "U")
+    if(n_A > 0) pars = c(pars,"A")
   }
   if (map_estimation == FALSE) {
     out <- rstan::sampling(

inst/stan/marss.stan

@@ -14,6 +14,8 @@ data {
   int<lower=0> row_indx_pos[n_pos];
   int<lower=0> est_trend; 
   int<lower=0> est_B;
+  int<lower=0> n_A;
+  int<lower=0> est_A[n_A+1];
   vector[n_pos] y; // data
   int y_int[n_pos];
   int family; // 1 = normal, 2 = binomial, 3 = poisson, 4 = gamma, 5 = lognormal
@@ -23,6 +25,7 @@ parameters {
   vector[S] pro_dev[N-1];
   vector[n_trends * est_trend] U;
   matrix[S*est_B,S*est_B] B;
+  vector[n_A] A; // offsets
   real<lower=0> sigma_process[S];
   real<lower=0> sigma_obs[n_obsvar];
 }
@@ -32,6 +35,10 @@ transformed parameters {
   //matrix[N, S] x;
   matrix[S,S] Bmat;
   vector[S] Uvec;
+  vector[M] Avec;
+
+  for(i in 1:M) Avec[i] = 0;
+  for(i in 1:n_A) Avec[est_A[i]] = A[i];
 
   for(i in 1:S) {
     if(est_trend) {
@@ -58,32 +65,27 @@ transformed parameters {
       x[t,] = x[t,] + Uvec;
     }
   }
-  // random walk in states
-  // for(s in 1:S) {
-  //  x[1,s] = x0[s]; // initial state, vague prior below
-  //  for(t in 2:N) {
-  //   x[t,s] = x[t-1,s] + pro_dev[t-1,s];
-  //   if(est_trend == 1) x[t,s] = x[t,s] + U[trends[s]];
-  //  }
-  // }
+
   // map predicted states to time series
   for(m in 1:M) {
     for(t in 1:N) {
-      pred[t,m] = x[t,states[m]];
+      pred[t,m] = x[t,states[m]] + Avec[m];
     }
   }
 }
 model {
-  //x0 ~ normal(0,10);
+  x0 ~ normal(0, 3); // initial states
+  A ~ normal(0, 3); // A offsets
+
   for(i in 1:n_obsvar) {
-    sigma_obs[i] ~ student_t(3,0,2);
+    sigma_obs[i] ~ student_t(3,0,1); // observation var sigma
   }
   for(s in 1:n_provar) {
-    sigma_process[s] ~ student_t(3,0,2); // process var
+    sigma_process[s] ~ student_t(3,0,1); // process var sigma
   }
   if(est_trend==1){
     for(i in 1:n_trends) {
-      U[i] ~ normal(0,1);
+      U[i] ~ normal(0,1); // optional trends
     }
   }
   for(s in 1:S) {