Menulis Model Stan di R
Tentukan blok data, parameter, dan blok model dalam sintaks Stan.
Menulis Model Stan di R adalah pelajaran R Academy gratis di CoddyKit. Ini adalah pelajaran 2 dari 4. Kamu bisa membaca pelajaran lengkapnya di bawah secara gratis — lalu praktikkan langsung di browser dengan editor kode bawaan dan tutor AI 24/7. Ini adalah bagian dari jalur belajar R Academy, dan progresmu tersinkronisasi di web dan aplikasi CoddyKit. Kursus R Academy mencakup 4 pelajaran total.
Apa Itu Stan?
Stan adalah bahasa pemrograman probabilistik untuk pemodelan statistik Bayesian. RStan adalah antarmuka R ke Stan. Anda menulis model dalam bahasa Stan (mirip C++), lalu Stan mengompilasinya menjadi kode C++ yang efisien dan menjalankan pengambilan sampel Hamiltonian Monte Carlo (HMC) untuk mendekati posterior.
# Stan installation check
library(rstan)
# Check version
cat('RStan version:', as.character(packageVersion('rstan')), '\n')
# Enable parallel chains
options(mc.cores = parallel::detectCores())
# Reuse compiled models across sessions
rstan_options(auto_write = TRUE)
cat('Stan ready. Cores:', parallel::detectCores(), '\n')Struktur Model Stan
Model Stan memiliki hingga enam blok bernama: data, transformed data, parameters, transformed parameters, model, dan generated quantities. Tiga blok esensialnya adalah data, parameters, dan model.
library(rstan)
# Stan model as a character string in R
stan_code <- '
data {
int<lower=0> N; // number of observations
vector[N] x; // predictor
vector[N] y; // response
}
parameters {
real alpha; // intercept
real beta; // slope
real<lower=0> sigma; // noise (must be positive)
}
model {
// Priors
alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ exponential(1);
// Likelihood
y ~ normal(alpha + beta * x, sigma);
}
'
cat('Stan model defined as a string in R\n')
cat('Blocks: data, parameters, model\n')data{}: Mendeklarasikan Data Masukan
Blok data mendeklarasikan semua data eksternal yang diterima model dari R. Jenis yang tersedia mencakup int, real, vector[N], matrix[M,N], dan array. Batasan seperti <lower=0> diperiksa saat program berjalan.
library(rstan)
# Comprehensive data block examples
data_block_examples <- '
data {
// Scalars
int<lower=1> N; // at least 1 observation
int<lower=2> K; // at least 2 groups
// Constrained scalars
real<lower=0, upper=1> rate; // probability
// Vectors
vector[N] y; // continuous response
array[N] int<lower=0, upper=1> z; // binary outcomes
// Matrix
matrix[N, K] X; // design matrix
// Integer array
array[N] int group; // group membership 1..K
}
'
cat(data_block_examples)parameters{}: Jenis Parameter
Blok parameters mendeklarasikan hal-hal yang tidak diketahui yang akan diambil sampelnya oleh Stan. Batasan dalam blok ini menentukan ruang parameter: <lower=0> untuk kuantitas positif, <lower=0, upper=1> untuk probabilitas. Stan secara otomatis menerapkan koreksi Jacobian probabilitas log untuk parameter yang dibatasi.
library(rstan)
# Common parameter declarations
param_examples <- '
parameters {
// Unconstrained
real mu; // mean
vector[K] beta; // regression coefficients
// Positive (sigma, lambda, variance)
real<lower=0> sigma;
real<lower=0> lambda;
// Probability
real<lower=0, upper=1> theta;
// Simplex (sums to 1, for mixture weights)
simplex[K] pi;
// Correlation matrix
corr_matrix[K] Omega;
// Cholesky factor of covariance
cholesky_factor_cov[K] L_Sigma;
}
'
cat(param_examples)model{}: Prior dan Likelihood
Blok model mengakumulasikan log-posterior melalui sintaks ~ (tilde). y ~ normal(mu, sigma) adalah bentuk singkat untuk menambahkan log PDF Normal ke target. Anda dapat menuliskan kenaikan probabilitas log secara eksplisit dengan target += normal_lpdf(y | mu, sigma).
library(rstan)
model_block_example <- '
model {
// --- Priors ---
mu ~ normal(0, 10); // weakly informative
sigma ~ cauchy(0, 2.5); // half-Cauchy for scale
beta ~ normal(0, 1); // standardised coefficients
// --- Likelihood ---
// Tilde notation (most common)
y ~ normal(mu + X * beta, sigma);
// Equivalent explicit notation:
// target += normal_lpdf(y | mu + X * beta, sigma);
// For loop (less common but valid)
// for (i in 1:N)
// target += normal_lpdf(y[i] | mu, sigma);
}
'
cat(model_block_example)Menyesuaikan Model Normal Sederhana
Berikut alur kerja Stan minimal yang lengkap: definisikan string model, siapkan daftar data, panggil stan(), lalu periksa hasilnya dengan print(). Stan mengompilasi model pada penggunaan pertama dan menyimpannya di tembolok untuk penggunaan berikutnya.
library(rstan)
# Stan model: estimate mean and SD of a normal distribution
normal_model <- '
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
mu ~ normal(0, 10);
sigma ~ exponential(0.1);
y ~ normal(mu, sigma);
}
'
# Simulate data
set.seed(42)
y_data <- rnorm(50, mean = 5, sd = 2)
# Fit the model
fit <- stan(
model_code = normal_model,
data = list(N = length(y_data), y = y_data),
chains = 2,
iter = 1000,
warmup = 500,
refresh = 0 # suppress iteration output
)
print(fit, pars = c('mu', 'sigma'))Regresi Linear dalam Stan
Stan membuat regresi linear Bayesian menjadi mudah: tentukan prior normal untuk koefisien dan prior eksponensial atau half-Cauchy untuk sigma. Posterior memberikan distribusi ketidakpastian lengkap untuk setiap koefisien, bukan hanya estimasi titik.
library(rstan)
lin_reg_model <- '
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ exponential(1);
y ~ normal(alpha + beta * x, sigma);
}
generated quantities {
vector[N] y_rep; // posterior predictive
for (i in 1:N)
y_rep[i] = normal_rng(alpha + beta * x[i], sigma);
}
'
set.seed(7)
n <- 80
x <- rnorm(n); y <- 2 + 3 * x + rnorm(n, 0, 1)
fit <- stan(model_code = lin_reg_model,
data = list(N = n, x = x, y = y),
chains = 2, iter = 1000, refresh = 0)
print(fit, pars = c('alpha', 'beta', 'sigma'))Blok transformed parameters{}
Blok transformed parameters menghitung kuantitas turunan dari parameter yang diambil sampelnya. Kuantitas ini berguna untuk memparametrisasi model dengan cara yang stabil secara numerik (misalnya, dekomposisi Cholesky dan transformasi log), sekaligus menjaga agar parameter utama tetap mudah diinterpretasikan.
library(rstan)
# Model with transformed parameters
trans_param_example <- '
data {
int<lower=0> N;
array[N] int<lower=0> y; // counts
}
parameters {
real log_lambda; // work in log space for stability
}
transformed parameters {
real<lower=0> lambda;
lambda = exp(log_lambda); // transform back to original scale
}
model {
log_lambda ~ normal(1, 2); // prior on log scale
y ~ poisson(lambda);
}
'
set.seed(42)
y_counts <- rpois(30, lambda = 5)
fit_poisson <- stan(
model_code = trans_param_example,
data = list(N = length(y_counts), y = y_counts),
chains = 2, iter = 1000, refresh = 0
)
print(fit_poisson, pars = c('log_lambda', 'lambda'))Blok generated quantities{}
Blok generated quantities dijalankan setelah pengambilan sampel untuk menghitung kuantitas tambahan: sampel prediktif posterior (y_rep), log-likelihood untuk LOO-CV, atau parameter yang ditransformasi untuk ringkasan. Nilai di sini diambil sampelnya dari distribusi prediktif posterior.
library(rstan)
# Use generated quantities for posterior predictive checks
model_with_gq <- '
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
mu ~ normal(0, 10);
sigma ~ exponential(0.5);
y ~ normal(mu, sigma);
}
generated quantities {
vector[N] y_rep; // replicated datasets
real mean_y_rep; // mean of replicated data
for (i in 1:N)
y_rep[i] = normal_rng(mu, sigma);
mean_y_rep = mean(y_rep);
}
'
set.seed(1)
y_obs <- rnorm(40, 3, 1.5)
fit <- stan(model_code = model_with_gq,
data = list(N = length(y_obs), y = y_obs),
chains = 2, iter = 1000, refresh = 0)
print(fit, pars = c('mu', 'sigma', 'mean_y_rep'))Meneruskan Data dari R ke Stan
Argumen data untuk stan() adalah daftar R bernama. Nama-namanya harus sama persis dengan nama variabel yang dideklarasikan dalam blok data{} Stan. Vektor menjadi vector[N] Stan; bilangan bulat menjadi int; matriks R menjadi matrix[M,N] Stan.
library(rstan)
# Data preparation: names must match Stan data block exactly
set.seed(42)
n <- 60
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 1.5 + 2 * x1 - 0.8 * x2 + rnorm(n, 0, 0.5)
# Build the design matrix
X <- cbind(x1, x2) # 60 x 2 matrix
# Named list passed to stan(data = ...)
stan_data <- list(
N = n, # int
K = ncol(X), # int
X = X, # matrix[N, K]
y = y # vector[N]
)
cat('Stan data list elements:\n')
for (nm in names(stan_data)) {
cat(' ', nm, ': class =', class(stan_data[[nm]]),
'dim =', paste(dim(stan_data[[nm]]), collapse = 'x'),
'\n')
}Melihat Model yang Dikompilasi
Setelah penyesuaian, gunakan print(fit) untuk melihat ringkasan posterior (mean, se_mean, sd, kuantil, Rhat, n_eff) untuk semua parameter. Gunakan stan_plot(fit) untuk plot interval visual dan traceplot(fit) untuk menilai pencampuran rantai.
library(rstan)
# Reuse the simple normal model from scene 6
normal_model <- '
data { int<lower=0> N; vector[N] y; }
parameters { real mu; real<lower=0> sigma; }
model {
mu ~ normal(0, 10);
sigma ~ exponential(0.1);
y ~ normal(mu, sigma);
}
'
set.seed(5)
fit <- stan(model_code = normal_model,
data = list(N = 50, y = rnorm(50, 7, 3)),
chains = 2, iter = 1000, refresh = 0)
# Detailed summary table
print(fit)
# Extract as data frame
posterior_df <- as.data.frame(fit)
cat('\nPosterior samples shape:', nrow(posterior_df),
'rows x', ncol(posterior_df), 'cols\n')Pemeriksaan Cepat
Dalam model Stan, Anda ingin membatasi suatu parameter agar bernilai positif secara ketat (misalnya, simpangan baku). Deklarasi mana yang benar?
Rangkuman: Menulis Model Stan
Hal-hal penting:
- Model Stan memiliki tiga blok penting:
data{},parameters{},model{} - Batasan:
<lower=0>,<upper=1>,<lower=0, upper=1>untuk probabilitas - Jenis:
int,real,vector[N],matrix[M,N],simplex[K] - Blok model: gunakan sintaks tilde
y ~ normal(mu, sigma)untuk prior dan likelihood transformed parameters{}untuk kuantitas turunan;generated quantities{}untuk proses setelah pengambilan sampel- Teruskan data sebagai daftar R bernama yang namanya sama persis dengan nama blok Stan
rstan_options(auto_write = TRUE)menyimpan model yang telah dikompilasi di tembolok
library(rstan)
# Stan model skeleton
model_skeleton <- '
data { int N; vector[N] y; }
parameters { real mu; real<lower=0> sigma; }
model { mu ~ normal(0,10); sigma ~ exponential(1); y ~ normal(mu, sigma); }
'
cat(model_skeleton)
cat('\nFit with: stan(model_code = model_skeleton, data = list(N=..., y=...), chains=4)\n')Pertanyaan yang Sering Diajukan
Apakah pelajaran “Menulis Model Stan di R” gratis?
Ya — teks lengkap “Menulis Model Stan di R” gratis dibaca di sini di web. Untuk praktiknya secara interaktif (editor kode bawaan dan tutor AI 24/7) dan buka sisa kursus R Academy, upgrade ke CoddyKit PRO. Kursus R Academy mencakup 4 pelajaran total.
Apa yang akan aku pelajari di “Menulis Model Stan di R”?
Tentukan blok data, parameter, dan blok model dalam sintaks Stan. Kamu berlatih R Academy dengan kode praktik yang langsung kamu jalankan di browser, dan tutor AI 24/7 menjawab pertanyaanmu saat kamu mengerjakan pelajaran ini.
Apakah aku perlu pengalaman untuk memulai R Academy?
Tidak diperlukan pengalaman sebelumnya. R Academy di CoddyKit dirancang untuk pemula hingga pelajar tingkat lanjut, jadi kamu bisa memulai di sini atau dari awal dan belajar sesuai kecepatan kamu sendiri. Ini adalah pelajaran 2 dari 4.
Berapa lama pelajaran “Menulis Model Stan di R” memakan waktu?
Sebagian besar pelajaran CoddyKit memakan waktu sekitar 5–10 menit. Setiap pelajaran ringkas dan interaktif, jadi kamu membuat kemajuan stabil dan melanjutkan dari tempat kamu tinggalkan di web dan aplikasi.
Bisakah aku menulis dan menjalankan kode dalam pelajaran R Academy ini?
Ya. Setiap pelajaran R Academy menyertakan editor kode bawaan, jadi kamu menulis dan menjalankan kode nyata langsung di browser dan mendapatkan umpan balik AI instan — tidak diperlukan penyiapan lokal.
Semua pelajaran dalam kursus ini
- Pengantar Pemikiran Bayesian
- Menulis Model Stan di R
- Pengambilan Sampel MCMC dan Diagnostik
- Pemeriksaan Prediktif Posterior