Rの関数:rmultireg {bayesm}

# ライブラリの読み込み
library(bayesm)
library(MASS) # 多変量正規分布の乱数生成に使用

# 乱数シードの設定
seed <- 20251209
set.seed(seed)

# データの次元設定
n_obs <- 500 # 観測数 (n)
n_vars <- 3 # 説明変数の数 (k, 切片含む)
n_resp <- 2 # 目的変数の数 (m)

# 説明変数行列 X の生成 (n x k)
# 第1列は切片項として1、残りは標準正規分布から生成
X <- cbind(1, matrix(rnorm(n_obs * (n_vars - 1)), nrow = n_obs))
colnames(X) <- c("切片", "変数1", "変数2")

# 真の回帰係数行列 B_true (k x m)
# これが推定のターゲットとなる真の値です
B_true <- matrix(c(
  2.0,  0.5, # 切片: Y1は2.0, Y2は0.5
  1.5, -1.0, # 変数1: Y1に正の影響, Y2に負の影響
  -0.8, 0.8 # 変数2: Y1に負の影響, Y2に正の影響
), nrow = n_vars, byrow = TRUE)
colnames(B_true) <- c("Y1", "Y2")
rownames(B_true) <- c("切片", "変数1", "変数2")

# 真の誤差分散共分散行列 Sigma_true (m x m)
# Y1とY2の間には相関があると仮定します
Sigma_true <- matrix(c(
  1.0, 0.5,
  0.5, 1.0
), nrow = n_resp)

# 誤差項 U の生成
U <- mvrnorm(n = n_obs, mu = rep(0, n_resp), Sigma = Sigma_true)

# 目的変数行列 Y の生成 (モデル: Y = XB + U)
Y <- X %*% B_true + U
colnames(Y) <- c("Y1", "Y2")

cat("--- 真のパラメータ B ---\n")
print(B_true)
cat("\n--- 真のパラメータ Sigma ---\n")
print(Sigma_true)
cat("\n")

# 事前分布のハイパーパラメータ設定
# ここでは比較的「情報の少ない」事前分布を設定します

# Bの事前平均 Bbar (0で設定)
Bbar <- matrix(0, nrow = n_vars, ncol = n_resp)

# Bの事前精度行列 A (小さな値 = 大きな分散 = 弱い事前情報)
A <- 0.01 * diag(n_vars)

# Sigmaの事前分布の自由度 nu (最小限の値)
nu <- n_resp + 3

# Sigmaの事前尺度行列 V (事前情報を弱めるため、単位行列を使用)
V <- nu * diag(n_resp)

# rmultireg は1回の呼び出しで事後分布から1つのサンプルを返します。
# 事後分布全体を把握するため、繰り返し呼び出して多数のサンプルを蓄積します。

R_draws <- 10000 # シミュレーション回数
B_samples <- array(0, dim = c(n_vars, n_resp, R_draws))
Sigma_samples <- array(0, dim = c(n_resp, n_resp, R_draws))

for (r in 1:R_draws) {
  # rmultireg の実行
  out <- rmultireg(Y, X, Bbar, A, nu, V)

  # 結果の格納
  B_samples[, , r] <- out$B
  Sigma_samples[, , r] <- out$Sigma
}

# 事後平均の計算
B_posterior_mean <- apply(B_samples, c(1, 2), mean)
Sigma_posterior_mean <- apply(Sigma_samples, c(1, 2), mean)

colnames(B_posterior_mean) <- colnames(Y)
rownames(B_posterior_mean) <- colnames(X)

cat("--- 推定された B の事後平均 ---\n")
print(round(B_posterior_mean, 4))

cat("\n--- 推定された Sigma の事後平均 ---\n")
print(round(Sigma_posterior_mean, 4))

# 変数名と目的変数名の取得
var_names <- rownames(B_posterior_mean)
resp_names <- colnames(B_posterior_mean)

# 各目的変数ごとにループして結果を表示
for (j in 1:length(resp_names)) {
  cat(paste0("=== 目的変数 ", resp_names[j], " に対する回帰係数の要約 ===\n"))

  # 結果を格納するための行列を作成
  # 列構成: 事後平均, 下側2.5%, 上側97.5%
  summary_table <- matrix(NA, nrow = length(var_names), ncol = 3)
  rownames(summary_table) <- var_names
  colnames(summary_table) <- c("事後平均", "下側2.5%", "上側97.5%")

  # 各説明変数について計算
  for (i in 1:length(var_names)) {
    # 該当するパラメータの全サンプルを抽出
    draws <- B_samples[i, j, ]

    # 平均と分位点の計算
    post_mean <- mean(draws)
    ci <- quantile(draws, probs = c(0.025, 0.975))

    # 行列に格納
    summary_table[i, 1] <- post_mean
    summary_table[i, 2] <- ci[1]
    summary_table[i, 3] <- ci[2]
  }

  # 結果の出力
  print(round(summary_table, 4))
  cat("\n")
}

# 真の値の再確認（比較用）
cat("=== 参考: 真のパラメータ値 (B_true) ===\n")
print(B_true)