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

int n;                    // Integer
real x;                   // Real number
complex z;                // Complex number

vector[N] v;              // Column vector
row_vector[N] r;          // Row vector
matrix[M, N] A;           // Matrix

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

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

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);

y ~ bernoulli(theta);
y ~ binomial(n, theta);
y ~ poisson(lambda);
y ~ neg_binomial_2(mu, phi);
y ~ categorical(theta);

y ~ multi_normal(mu, Sigma);        // Sigma is COVARIANCE
y ~ multi_normal_cholesky(mu, L);
y ~ lkj_corr(eta);

// GOOD - Efficient
y ~ normal(mu, sigma);

// BAD - Slow
for (n in 1:N) y[n] ~ normal(mu[n], sigma);

parameters {
  vector[J] theta_raw;
}
transformed parameters {
  vector[J] theta = mu + tau * theta_raw;
}
model {
  theta_raw ~ std_normal();
}

// These are equivalent:
y ~ normal(mu, sigma);
target += normal_lpdf(y | mu, sigma);

// 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

library(cmdstanr)
mod <- cmdstan_model("model.stan")
fit <- mod$sample(data = stan_data, chains = 4)
fit$summary()
fit$cmdstan_diagnose()

# 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?

fit <- mod$sample(data = stan_data, chains = 4, adapt_delta = 0.95)

fit$summary()              # Rhat, ESS
fit$cmdstan_diagnose()     # Divergences, treedepth
library(bayesplot)
mcmc_rank_hist(fit$draws()) # Ranked traceplots (preferred)

y_rep <- fit$draws("y_rep", format = "matrix")
library(bayesplot)
ppc_dens_overlay(y, y_rep[1:100, ])

library(loo)
loo1 <- loo(fit1$draws("log_lik"))
loo2 <- loo(fit2$draws("log_lik"))
loo_compare(loo1, loo2)

# 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))

# 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 {
  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);
  }
}