add link to r-universe

comp_labs.Rmd

@@ -19,9 +19,9 @@ Thursday from 11:20-12:20 FISH 213
 
 ### Resources
 
-* [FISH 550 lab book](https://atsa-es.github.io/atsa-labs/)
+* [FISH 550 lab book](https://atsa-es.github.io/atsa-labs/) &  [atsa-es package repository](https://atsa-es.r-universe.dev/packages)
 
-* [GitHub repo for assignments](https://github.com/atsa-es/fish550-2023)
+* [GitHub repo for assignments](https://github.com/atsa-es/fish550-2023) & [Team Lab Write-ups](https://atsa-es.github.io/fish550-2023/)
 
 * [Discusssion board](https://github.com/atsa-es/fish550-2023/discussions)
 

docs/Lectures/TMB/uni.hpp

@@ -0,0 +1,81 @@
+/// @file uni.hpp
+// univariate time series
+
+#ifndef uni_hpp
+#define uni_hpp
+
+#include "marssTMB/LOM.hpp"
+
+#undef TMB_OBJECTIVE_PTR
+#define TMB_OBJECTIVE_PTR obj
+
+template<class Type>
+Type uni(objective_function<Type>* obj) {
+
+  using namespace density;
+
+  DATA_VECTOR(Y);
+  DATA_INTEGER(n);
+  DATA_INTEGER(est_drift);
+  DATA_INTEGER(est_rho);
+  DATA_IVECTOR(keep);
+  DATA_STRUCT(par, LOM); // list of model matrices
+
+  PARAMETER(u); // drift
+  PARAMETER(logit_rho); // inv-logit(b)
+  PARAMETER(log_obs_sd); // obs error
+  PARAMETER(log_pro_sd); // pro error
+  PARAMETER_VECTOR(x); // n x 1
+  //PARAMETER(x0);
+
+  Type obs_sigma = exp(log_obs_sd);
+  Type pro_sigma = exp(log_pro_sd);
+  Type rho;
+  if(est_rho==1) rho = exp(logit_rho)/(1+exp(logit_rho));
+
+  // random effects / penalties
+  Type nll=0; // total function to maximize
+  vector<Type> pred(n); // predictions
+
+  // process deviations
+  for(int i = 0; i < n; i++) {
+    nll -= dnorm(x(i), Type(0.0), pro_sigma, true);
+  }
+  // initial conditions
+  pred(0) = x(0);
+  // process / random walk
+  for(int i = 1; i < n; i++) {
+    pred(i) = pred(i-1);
+    if(est_rho == 1) pred(i) = rho * pred(i);
+    if(est_drift == 1) pred(i) = pred(i) + u;
+    pred(i) = pred(i) + x(i); // process stochasticity
+  }
+
+  // observation likelihood
+  for(int i = 0; i < n; i++) {
+    if(keep(i) == 1) nll -= dnorm(Y(i), pred(i), obs_sigma, true); 
+  }
+
+  if(est_rho) {
+    ADREPORT(rho); 
+    REPORT(rho); 
+  }
+  if(est_drift) {
+    ADREPORT(u); 
+    REPORT(u); 
+  }
+  ADREPORT(obs_sigma); 
+  REPORT(obs_sigma);
+  ADREPORT(pro_sigma); 
+  REPORT(pro_sigma);
+  ADREPORT(pred); 
+  REPORT(pred);
+  REPORT(par(0));
+  // end
+  return (nll);
+}
+
+#undef TMB_OBJECTIVE_PTR
+#define TMB_OBJECTIVE_PTR this
+
+#endif

docs/Lectures/TMB/uniTMB.R

@@ -0,0 +1,68 @@
+#' Fit a univariate MARSS model with TMB.
+#' 
+#' @param y the data. Can have NAs.
+#' @param estimate_drift estimate the u parameter
+#' @param estimate_rho estimate b parameter
+#'
+#' @return A data frame with estimates and se's
+#' @example inst/examples/uni_example.R
+#' @author Eric Ward and edited by Eli Holmes.
+#' @export
+uniTMB <- function(y, estimate_drift = TRUE, estimate_rho = FALSE){
+
+parameters <- list(
+  log_obs_sd = 0,
+  log_pro_sd = 0,
+  x = rep(0, length(y)),
+  u = 0,
+  #x0 = 0,
+  logit_rho = 0
+)
+
+# Map off parameters not being estimated
+tmb_map <- list(x = as.factor(c(NA,1:(length(y)-1))))
+if(estimate_drift == FALSE) tmb_map <- c(tmb_map, list(u = factor(NA)))
+if(estimate_rho == FALSE) tmb_map <- c(tmb_map, list(logit_rho = factor(NA)))
+
+# Create TMB data
+data_list <- list(Y = y, n = length(y),
+                  est_drift = as.numeric(estimate_drift),
+                  est_rho = as.numeric(estimate_rho),
+                  keep = ifelse(!is.na(y),1,0),
+                  par = list(matrix(1,2,2), matrix(2,2,2)),
+                  model = "uni")
+
+# Create object for fitting
+obj <- TMB::MakeADFun(
+  data = data_list,
+  map = tmb_map,
+  random = "x",
+  parameters = parameters,
+  DLL = "marssTMB_TMBExports",
+  silent = TRUE
+)
+
+# Do the fitting with stats::nlminb, sometimes need to change default control args if not converging
+pars <- stats::nlminb(
+  start = obj$par, objective = obj$fn,
+  gradient = obj$gr
+)
+
+par_summary <- summary(TMB::sdreport(obj))
+
+indx <- grep("pred", rownames(par_summary))
+df <- data.frame(
+  pred = as.numeric(par_summary[indx,"Estimate"]),
+  se = as.numeric(par_summary[indx,"Std. Error"]),
+  y = y,
+  t = 1:length(y)
+)
+vals <- c("u", "obs_sigma", "pro_sigma")
+pardf <- data.frame(
+  name = vals,
+  estimate = as.numeric(par_summary[vals, "Estimate"]),
+  se = as.numeric(par_summary[vals,"Std. Error"])
+)
+
+return(list(df = df, coef = pardf, all = par_summary))
+}