R Academy · Lekcja

Pisanie modeli Stan w R

Definiuj bloki danych, parametry i blok modelu w składni Stan.

Pisanie modeli Stan w R to bezpłatna lekcja R Academy na CoddyKit. To lekcja 2 z 4. Możesz przeczytać całą lekcję poniżej za darmo — a potem ćwiczyć ją interaktywnie w przeglądarce z wbudowanym edytorem kodu i tutorem AI dostępnym 24/7. To część ścieżki edukacyjnej R Academy, a Twój postęp synchronizuje się między webem a aplikacją CoddyKit. Kurs R Academy zawiera 4 lekcji w sumie.

Czym jest Stan

Stan to probabilistyczny język programowania służący do modelowania statystycznego metodami bayesowskimi. RStan jest interfejsem R do systemu Stan. Model zapisujesz w języku Stan (podobnym do C++), a Stan kompiluje go do wydajnego kodu C++, który wykonuje próbkowanie metodą Monte Carlo z hamiltonowskim łańcuchem Markowa (HMC), aby przybliżyć rozkład a posteriori.

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

Struktura modelu Stan

Model Stan może zawierać maksymalnie sześć nazwanych bloków: data, transformed data, parameters, transformed parameters, model oraz generated quantities. Trzy niezbędne bloki to data, parameters i 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{}: deklarowanie danych wejściowych

Blok data deklaruje wszystkie dane zewnętrzne przekazywane do modelu z R. Dostępne typy obejmują int, real, vector[N], matrix[M,N] oraz array. Ograniczenia takie jak <lower=0> są sprawdzane w czasie działania programu.

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{}: Typy parametrów

Blok parameters deklaruje niewiadome, które Stan będzie próbkować. Ograniczenia w tym bloku definiują przestrzeń parametrów: <lower=0> dla wartości dodatnich, <lower=0, upper=1> dla prawdopodobieństw. Stan automatycznie stosuje poprawki jakobianu logarytmu prawdopodobieństwa dla parametrów z ograniczeniami.

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{}: Rozkłady a priori i wiarygodność

Blok model kumuluje logarytm rozkładu a posteriori za pomocą składni ~ (tyldy). y ~ normal(mu, sigma) to skrót oznaczający dodanie logarytmu funkcji gęstości rozkładu normalnego do target. Jawne przyrosty logarytmu prawdopodobieństwa można zapisywać za pomocą 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)

Dopasowywanie prostego modelu normalnego

Oto kompletny minimalny przebieg pracy ze Stanem: należy zdefiniować ciąg modelu, przygotować listę danych, wywołać stan() i przejrzeć wyniki za pomocą print(). Stan kompiluje model przy pierwszym użyciu i zapisuje go w pamięci podręcznej na potrzeby kolejnych uruchomień.

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

Regresja liniowa w Stan

Stan upraszcza bayesowską regresję liniową: należy określić normalne rozkłady a priori dla współczynników oraz wykładniczy lub półrozkład Cauchy’ego a priori dla sigma. Rozkład a posteriori daje pełny rozkład niepewności dla każdego współczynnika, a nie tylko estymaty punktowe.

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 oblicza wielkości pochodne na podstawie próbkowanych parametrów. Są one przydatne do parametryzowania modeli w sposób stabilny numerycznie (np. za pomocą rozkładów Cholesky’ego i transformacji logarytmicznych), przy jednoczesnym zachowaniu interpretowalności parametrów podstawowych.

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 jest wykonywany po próbkowaniu w celu obliczenia dodatkowych wielkości: próbek predykcyjnych a posteriori (y_rep), logarytmu wiarygodności na potrzeby LOO-CV lub przekształconych parametrów do podsumowań. Wartości w tym bloku są próbkowane z predykcyjnego rozkładu a posteriori.

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

Przekazywanie danych z R do Stan

Argument data funkcji stan() jest nazwaną listą R. Nazwy muszą dokładnie odpowiadać nazwom zmiennych zadeklarowanych w bloku Stan data{}. Wektory stają się obiektami Stan typu vector[N], liczby całkowite — typu int, a macierze R — obiektami Stan typu matrix[M,N].

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')
}

Wyświetlanie skompilowanego modelu

Po dopasowaniu użyj print(fit), aby wyświetlić podsumowania rozkładu a posteriori (mean, se_mean, sd, quantiles, Rhat, n_eff) dla wszystkich parametrów. Użyj stan_plot(fit), aby utworzyć wykres przedziałów, oraz traceplot(fit), aby ocenić mieszanie się łańcuchów.

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

Szybkie sprawdzenie

W modelu Stan chcesz wymusić, aby parametr był ściśle dodatni (np. odchylenie standardowe). Która deklaracja jest poprawna?

Podsumowanie: pisanie modeli Stan

Najważniejsze informacje:

Model Stan ma trzy podstawowe bloki: data{}, parameters{}, model{}
Ograniczenia: <lower=0>, <upper=1>, <lower=0, upper=1> dla prawdopodobieństw
Typy: int, real, vector[N], matrix[M,N], simplex[K]
Blok modelu: używaj składni tyldy y ~ normal(mu, sigma) dla rozkładów a priori i wiarygodności
transformed parameters{} służy do wielkości pochodnych, a generated quantities{} do obliczeń po próbkowaniu
Przekazuj dane jako nazwaną listę R, której nazwy dokładnie odpowiadają nazwom bloków Stan
rstan_options(auto_write = TRUE) zapisuje skompilowane modele w pamięci podręcznej

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

Często zadawane pytania

Czy lekcja „Pisanie modeli Stan w R” jest bezpłatna?

Tak — pełny tekst „Pisanie modeli Stan w R” jest dostępny za darmo tutaj w sieci. Aby ćwiczyć ją interaktywnie (wbudowany edytor kodu i tutor AI dostępny 24/7) i odblokować resztę kursu R Academy, przejdź na CoddyKit PRO. Kurs R Academy zawiera 4 lekcji w sumie.

Co nauczysz się w „Pisanie modeli Stan w R”?

Definiuj bloki danych, parametry i blok modelu w składni Stan. Ćwiczysz R Academy z praktycznym kodem, który uruchamiasz bezpośrednio w przeglądarce, a tutor AI dostępny 24/7 odpowiada na Twoje pytania podczas pracy nad lekcją.

Czy potrzebuję doświadczenia, aby zacząć R Academy?

Nie wymagamy żadnego doświadczenia. R Academy w CoddyKit jest strukturyzowany dla początkujących i zaawansowanych użytkowników, więc możesz zacząć tutaj lub od początku i uczyć się w swoim tempie. To lekcja 2 z 4.

Ile czasu zajmuje lekcja „Pisanie modeli Stan w R”?

Większość lekcji CoddyKit trwa około 5–10 minut. Każda lekcja to mały, interaktywny krok, dzięki czemu robisz systematyczne postępy i zawsze wracasz dokładnie do tego samego miejsca — na webie i w aplikacji.

Czy mogę pisać i uruchamiać kod w tej lekcji R Academy?

Tak. Każda lekcja R Academy zawiera wbudowany edytor kodu, więc piszesz i uruchamiasz prawdziwy kod bezpośrednio w przeglądarce i od razu otrzymujesz sprzężenie zwrotne od AI — bez konfiguracji na komputerze.

Wszystkie lekcje w tym kursie

← Powrót do R Academy