0Pricing
R Academy · レッスン

R で Stan モデルを書く

Stan 構文でデータブロック、パラメーター、モデルブロックを定義します。

「R で Stan モデルを書く」はCoddyKit上の無料R Academyレッスンです。 これはレッスン2/4です。 下記で完全なレッスンを無料で読むことができます。その後、ブラウザ内の組み込みコードエディタと24時間対応のAIチューターでハンズオン演習できます。 これはR Academy学習パスの一部であり、ウェブとCoddyKitアプリ全体で進捗が同期されます。 R Academyコースには全4レッスンが含まれています。

Stanとは

Stanは、ベイズ統計モデリングのための確率的プログラミング言語です。RStanはStanのR用インターフェースです。モデルをStanの言語（C++に似た言語）で記述すると、Stanが効率的なC++コードにコンパイルし、ハミルトニアンモンテカルロ（HMC）サンプリングを実行して事後分布を近似します。

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

Stanモデルの構造

Stanモデルには、最大6つの名前付きブロックがあります。data、transformed data、parameters、transformed parameters、model、generated quantitiesです。必須のブロックはdata、parameters、modelの3つです。

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{}：入力データの宣言

dataブロックでは、Rからモデルに渡されるすべての外部データを宣言します。型にはint、real、vector[N]、matrix[M,N]、arrayなどがあります。<lower=0>のような制約は、実行時にチェックされます。

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{}：パラメータの型

parametersブロックでは、Stanがサンプリングする未知量を宣言します。このブロックの制約によってパラメータ空間を定義します。正の量には<lower=0>、確率には<lower=0, upper=1>を指定します。Stanは、制約付きパラメータに対する対数確率のヤコビアン補正を自動的に適用します。

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{}：事前分布と尤度

modelブロックでは、~（チルダ）構文を使って対数事後分布を累積します。y ~ normal(mu, sigma)は、正規分布のPDFの対数をtargetに加える処理の省略記法です。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)

単純な正規モデルのフィッティング

完全で最小限のStanワークフローを示します。モデル文字列を定義し、データリストを準備して、stan()を呼び出し、print()で結果を確認します。Stanは初回にモデルをコンパイルし、以降の実行のためにキャッシュします。

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

Stanでの線形回帰

Stanを使うと、ベイズ線形回帰を簡単に実行できます。係数に正規分布の事前分布を指定し、sigmaには指数分布または半Cauchy分布の事前分布を指定します。事後分布からは、点推定値だけでなく、各係数の不確実性の完全な分布が得られます。

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

transformed parameters{}ブロック

transformed parametersブロックでは、サンプリングしたパラメータから派生量を計算します。主なパラメータを解釈しやすい状態に保ちながら、数値的に安定した形でモデルをパラメータ化する場合（例：コレスキー分解や対数変換）に便利です。

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

generated quantities{}ブロック

generated quantitiesブロックはサンプリング後に実行され、事後予測サンプル（y_rep）、LOO-CV用の対数尤度、要約用の変換パラメータなど、追加の量を計算します。ここでの値は事後予測分布からサンプリングされます。

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

RからStanへのデータの渡し方

stan()のdata引数には、名前付きのRリストを指定します。名前は、Stanのdata{}ブロックで宣言した変数名と完全に一致していなければなりません。ベクトルはStanのvector[N]に、整数はintに、Rの行列はStanの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')
}

コンパイル済みモデルの表示

フィッティング後にprint(fit)を使うと、すべてのパラメータについて事後分布の要約（mean、se_mean、sd、quantiles、Rhat、n_eff）を確認できます。stan_plot(fit)では区間プロットを、traceplot(fit)ではチェーンの混合状態を確認できます。

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

確認問題

Stanモデルで、パラメータを厳密に正（例：標準偏差）に制約したいとします。正しい宣言はどれですか。

まとめ：Stanモデルの記述

重要なポイント：

  • Stanモデルには、基本となる3つのブロックがあります：data{}、parameters{}、model{}
  • 制約：<lower=0>、<upper=1>、確率には<lower=0, upper=1>
  • 型：int、real、vector[N]、matrix[M,N]、simplex[K]
  • modelブロック：事前分布と尤度にはチルダ構文y ~ normal(mu, sigma)を使用します
  • 派生量にはtransformed parameters{}を、サンプリング後の計算にはgenerated quantities{}を使用します
  • Stanのブロック名と完全に一致する名前付きRリストとしてデータを渡します
  • rstan_options(auto_write = TRUE)によってコンパイル済みモデルがキャッシュされます
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')

よくある質問

「R で Stan モデルを書く」レッスンは無料ですか？

はい。「R で Stan モデルを書く」の完全なテキストはこのウェブで無料で読めます。インタラクティブに演習し（組み込みコードエディタと24時間対応のAIチューター）、R Academyコースの残りをアンロックするには、CoddyKit PROにアップグレードしてください。 R Academyコースには全4レッスンが含まれています。

「R で Stan モデルを書く」で何を学びますか？

Stan 構文でデータブロック、パラメーター、モデルブロックを定義します。 ブラウザで直接実行するハンズオンコードでR Academyを演習し、24時間対応のAIチューターがレッスンを進める中での質問に答えます。

R Academyを始めるのに経験は必要ですか？

事前経験は必要ありません。CoddyKitのR Academyは初級者から上級者向けに構成されているため、ここから始めるか最初から始めて、自分のペースで進むことができます。 これはレッスン2/4です。

「R で Stan モデルを書く」レッスンにはどのくらい時間がかかりますか？

ほとんどのCoddyKitレッスンは約5～10分かかります。各レッスンはコンパクトでインタラクティブなので、着実に進歩し、ウェブとアプリ全体で正確に前回の場所から再開できます。

このR Academyレッスンでコードを書いて実行できますか？

はい。すべてのR Academyレッスンに組み込みコードエディタが含まれているため、ブラウザでリアルコードを書いて実行し、即座のAIフィードバックを取得できます。ローカル設定は不要です。

このコースのすべてのレッスン

  1. ベイズ的思考の基礎
  2. R で Stan モデルを書く
  3. MCMC サンプリングと診断
  4. 事後予測チェック
← R Academyに戻る