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)

# 1. データ数（営業日数）
N <- 30

# 2. 真のパラメータを設定
true_intercept <- 50 # 基本売上（千円）
true_beta <- c(
  1.0, # 気温が1度上がるごとの売上増（千円）
  15.0, # 広告費1万円あたりの売上増（千円）
  20.0 # イベント開催時の売上増（千円）
)
true_sigma <- 5.0 # 誤差の標準偏差（千円）

# 3. 説明変数を生成
# 気温 (15〜35℃)
temperature <- runif(N, 15, 35)
# SNS広告費 (0〜5万円)
ad_cost <- runif(N, 0, 5)
# イベントの有無 (開催確率20%)
is_event <- rbinom(N, 1, 0.2)

# 4. 目的変数（売上）を生成
# 行列形式で説明変数をまとめる
X <- cbind(temperature, ad_cost, is_event)

# 売上 = 切片 + X * β + 誤差
sales <- true_intercept + X %*% true_beta + rnorm(N, mean = 0, sd = true_sigma)

# データフレームにまとめる
sim_data <- data.frame(
  sales = as.numeric(sales),
  temperature,
  ad_cost,
  is_event
)

cat("--- 生成したデータの一部を確認 ---\n")
head(sim_data)

model_code <- "
data {
  int<lower=0> N;          // データ数
  int<lower=0> K;          // 説明変数の数
  matrix[N, K] X;          // 説明変数の行列
  vector[N] y;             // 目的変数
}

parameters {
  real intercept;          // 切片
  vector[K] beta;          // 回帰係数
  real<lower=0> sigma;     // 誤差の標準偏差
}

model {
  // --- 事前分布 ---
  // 正規分布を仮定（弱情報事前分布）
  intercept ~ normal(0, 100);
  // スケールパラメータにはコーシー分布を仮定
  sigma ~ cauchy(0, 5);
  // 回帰係数(beta)に個別の事前分布を設定
  // オーナーの「肌感覚」をモデルに組み込む
  // beta[1]: 気温の効果。平均1.0(千円), 標準偏差0.5(千円)の正規分布を仮定
  beta[1] ~ normal(1.0, 0.5);
  // beta[2]: 広告費の効果。特に事前知識がないため、弱情報事前分布を適用
  beta[2] ~ normal(0, 100);
  // beta[3]: イベント開催の効果。平均20.0(千円), 標準偏差2.5(千円)の正規分布を仮定
  beta[3] ~ normal(20.0, 2.5);
  // --- 尤度（データが生成されるモデル） ---
  // yは、平均が (切片 + X*beta) で、標準偏差がsigmaの正規分布に従う
  y ~ normal(intercept + X * beta, sigma);
}

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

# Stanに渡すためのデータをリスト形式で準備
stan_data <- list(
  N = N,
  K = ncol(X),
  X = X,
  y =  as.vector(sales)
)

# StanモデルのコンパイルとMCMCの実行
fit <- stan(
  model_code = model_code,
  data = stan_data, # 入力データ
  seed = seed, # 乱数シード
  chains = 4, # MCMCのチェーン数
  iter = 2000, # 1チェーンあたりの繰り返し数
  warmup = 1000, # バーンイン期間
  thin = 1 # 間引き数
)

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