stan-fundamentals

# Stan Fundamentals ## When to Use This Skill - Writing new Stan models from scratch - Understanding Stan program structure - Learning Stan syntax and conventions - Translating models from other languages to Stan - Optimizing existing Stan code ## Program Structure Stan models have up to 7 blocks in this exact order: ```stan functions { } // User-defined functions data { } // Input data declarations transformed data { } // Data preprocessing parameters { } // Model parameters transformed parameters { } // Derived parameters model { } // Log probability generated quantities { } // Posterior predictions ``` All blocks are optional. Empty string is valid (but useless) Stan program. ## Type System Quick Reference ### Scalars ```stan int n; // Integer real x; // Real number complex z; // Complex number ``` ### Vectors and Matrices ```stan vector[N] v; // Column vector row_vector[N] r; // Row vector matrix[M, N] A; // Matrix ``` ### Arrays (Modern Syntax) ```stan array[N] real x; // 1D array of reals array[M, N] int y; // 2D array of integers array[J] vector[K] theta; // Array of vectors ``` ### Constrained Types ```stan real<lower=0> sigma; // Non-negative real<lower=0, upper=1> p; // Probability simplex[K] theta; // Sums to 1 ordered[K] c; // Ascending corr_matrix[K] Omega; // Correlation cov_matrix[K] Sigma; // Covariance cholesky_factor_corr[K] L_Omega; // Cholesky correlation ``` ## Key Distributions ### Continuous (SD parameterization!) ```stan y ~ normal(mu, sigma); // sigma is SD y ~ student_t(nu, mu, sigma); y ~ cauchy(mu, sigma); y ~ exponential(lambda); y ~ gamma(alpha, beta); y ~ beta(a, b); y ~ lognormal(mu, sigma); ``` ### Discrete ```stan y ~ bernoulli(theta); y ~ binomial(n, theta); y ~ poisson(lambda); y ~ neg_binomial_2(mu, phi); y ~ categorical(theta); ``` ### Multivariate ```stan y ~ multi_normal(mu, Sigma); // Sigma is COVARIANCE y ~ multi_normal_cholesky(mu, L); y ~ lkj_corr(eta); ``` ## Essential Patterns ### Vectorization ```stan // GOOD - Efficient y ~ normal(mu, sigma); // BAD - Slow for (n in 1:N) y[n] ~ normal(mu[n], sigma); ``` ### Non-Centered Parameterization ```stan parameters { vector[J] theta_raw; } transformed parameters { vector[J] theta = mu + tau * theta_raw; } model { theta_raw ~ std_normal(); } ``` ### Target Syntax ```stan // These are equivalent: y ~ normal(mu, sigma); target += normal_lpdf(y | mu, sigma); ``` ## Common Priors ```stan // Location parameters mu ~ normal(0, 10); // Scale parameters sigma ~ exponential(1); sigma ~ cauchy(0, 2.5); // half-Cauchy when sigma > 0 // Probabilities theta ~ beta(1, 1); // Uniform on (0,1) // Regression coefficients beta ~ normal(0, 2.5); // Correlation matrices Omega ~ lkj_corr(2); // eta=2 favors identity ``` ## R Integration (cmdstanr) ```r library(cmdstanr) mod <- cmdstan_model("model.stan") fit <- mod$sample(data = stan_data, chains = 4) fit$summary() fit$cmdstan_diagnose() ``` ## Bayesian Workflow (Statistical Rethinking) ### 1. Prior Predictive Check ```r # Simulate from priors before fitting n_sim <- 1000 prior_alpha <- rnorm(n_sim, 0, 10) prior_sigma <- rexp(n_sim, 1) # Plot: do these produce sensible y values? ``` ### 2. Fit Model ```r fit <- mod$sample(data = stan_data, chains = 4, adapt_delta = 0.95) ``` ### 3. Diagnostics ```r fit$summary() # Rhat, ESS fit$cmdstan_diagnose() # Divergences, treedepth library(bayesplot) mcmc_rank_hist(fit$draws()) # Ranked traceplots (preferred) ``` ### 4. Posterior Predictive Check ```r y_rep <- fit$draws("y_rep", format = "matrix") library(bayesplot) ppc_dens_overlay(y, y_rep[1:100, ]) ``` ### 5. Model Comparison ```r library(loo) loo1 <- loo(fit1$draws("log_lik")) loo2 <- loo(fit2$draws("log_lik")) loo_compare(loo1, loo2) ``` ## link vs sim Pattern ### link(): Uncertainty in mu (epistemic) ```r # Posterior of expected value post <- fit$draws(format = "df") mu <- post$alpha + post$beta * x_new # Matrix of mu samples mu_PI <- apply(mu, 2, quantile, c(0.055, 0.945)) ``` ### sim(): Prediction interval (epistemic + aleatoric) ```r # Includes observation noise y_sim <- rnorm(n_samples, mu, post$sigma) y_PI <- apply(y_sim, 2, quantile, c(0.055, 0.945)) ``` ## Generated Quantities Template Always include for diagnostics and model comparison: ```stan generated quantities { vector[N] log_lik; // For LOO/WAIC array[N] real y_rep; // For posterior predictive checks for (n in 1:N) { log_lik[n] = normal_lpdf(y[n] | mu[n], sigma); y_rep[n] = normal_rng(mu[n], sigma); } } ``` ## Diagnostic Checklist - [ ] Rhat < 1.01 for all parameters - [ ] ESS_bulk > 400 - [ ] ESS_tail > 400 - [ ] Zero divergences - [ ] Not hitting max_treedepth - [ ] Prior predictive produces sensible values - [ ] Posterior predictive matches data pattern ## Key Differences from BUGS | Feature | Stan | BUGS/JAGS | |---------|------|-----------| | Normal | `normal(mu, sigma)` SD | `dnorm(mu, tau)` precision | | MVN | `multi_normal(mu, Sigma)` cov | `dmnorm(mu, Omega)` precision | | Execution | Sequential (order matters) | Declarative (order doesn't matter) | | Sampling | HMC/NUTS | Gibbs/Metropolis |