RStan:単回帰分析と外挿

# RStanの読み込み
library(rstan)

stan_output <- "D:/stan_output"

# Stanの並列計算設定
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# rstan のバージョン確認
cat("--- パッケージ rstan のバージョン確認 ---\n")
sapply(X = c("rstan"), packageVersion)

seed <- 20250629
set.seed(seed)

# ---- シミュレーションデータの作成 ----
# ストーリーに基づいたパラメータを設定
N <- 30 # 観測日数
alpha_true <- 20 # 気温0度の時の基礎売上（切片）
beta_true <- 5 # 気温が1度上がるごとの売上増加数（傾き）
sigma_true <- 8 # 日々のばらつき（標準偏差）

# 説明変数：気温データ (15度〜25度の範囲で観測)
temperature <- runif(N, 15, 25)

# 目的変数：クロワッサンの売上個数データ
# y = a + b*x + ε (εは正規分布に従う誤差)
sales <- round(alpha_true + beta_true * temperature + rnorm(N, 0, sigma_true))
sales[sales < 0] <- 0 # 売上はマイナスにならない

# データフレームにまとめる
bakery_data <- data.frame(temperature, sales)

# 生成したデータの一部を確認
head(bakery_data)

model_code <- "
data {
  int<lower=0> N;          // データ数
  vector[N] x;             // 説明変数 (temperature)
  vector[N] y;             // 目的変数 (sales)
}

parameters {
  real alpha;              // 切片 (intercept)
  real beta;               // 傾き (slope)
  real<lower=0> sigma;     // 誤差の標準偏差
}

model {
  // 事前分布 (比較的情報量の少ない弱情報事前分布)
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  sigma ~ cauchy(0, 5);

  // 尤度 (単回帰モデルの本体)
  // y は 平均 alpha + beta * x, 標準偏差 sigma の正規分布に従う
  y ~ normal(alpha + beta * x, sigma);
}

generated quantities {
    // 予測分布を生成するための部分
    vector[N] y_pred;
    for (n in 1:N) {
        y_pred[n] = normal_rng(alpha + beta * x[n], sigma);
    }
}
"

# Stanに渡すためのデータをリスト形式で準備
data_list <- list(
  N = nrow(bakery_data),
  x = bakery_data$temperature,
  y = bakery_data$sales
)

# StanモデルのコンパイルとMCMCの実行
fit <- stan(
  model_code = model_code,
  data = data_list, # データ
  seed = seed, # 乱数シード
  iter = 2000, # サンプリング回数
  warmup = 1000, # バーンイン期間
  chains = 4 # MCMCのチェーン数
)

# stanfit オブジェクトの保存
setwd(stan_output)
saveRDS(object = fit, file = "stan_fit.rds")

# stanfit オブジェクトの読み込み
setwd(stan_output)
fit <- readRDS("stan_fit.rds")

cat("--- 結果の要約 ---\n")
summary(fit, pars = c("alpha", "beta", "sigma"), probs = c(0.025, 0.5, 0.975))$summary
sampler_params <- get_sampler_params(fit, inc_warmup = FALSE)
divergences <- sum(sapply(sampler_params, function(x) sum(x[, "divergent__"])))
cat("\nDivergent Transitionsの数:", divergences, "\n")

# 予測したい気温のベクトル
target_temps <- c(28, 30, 35)

# 各気温での予測分布を生成し、要約する
# lapply を使って、各気温に対して同じ処理を繰り返す
prediction_summary_list <- lapply(target_temps, function(temp) {
  # 1. 予測の平均(mu)を各MCMCサンプルで計算
  #    mcmc_samplesにはalphaとbetaのサンプルが4000個入っている
  pred_mu <- mcmc_samples$alpha + mcmc_samples$beta * temp

  # 2. 予測分布(y_pred)を生成
  #    各MCMCサンプルのパラメータ(pred_mu, sigma)を使って、
  #    正規分布から新しいデータをシミュレーションする
  y_pred <- rnorm(
    n = length(pred_mu), # サンプル数 (4000)
    mean = pred_mu,
    sd = mcmc_samples$sigma
  )

  # 3. 予測分布の要約統計量を計算し、データフレームとして返す
  #    売上個数は整数なのでround()で丸める
  data.frame(
    "気温(度)" = temp,
    "予測売上_平均" = round(mean(y_pred)),
    "予測売上_中央値" = round(median(y_pred)),
    "95%予測区間_下限" = round(quantile(y_pred, 0.025)),
    "95%予測区間_上限" = round(quantile(y_pred, 0.975)),
    check.names = FALSE, # 列名に特殊文字を許す
    row.names = NULL
  )
})

# リストを一つのデータフレームにまとめる
prediction_summary_df <- do.call(rbind, prediction_summary_list)

cat("--- 予測結果 ---\n")
print(prediction_summary_df)