# 1. データをリセットして歩数を設定
rm(list = ls())
N <- 50000 

# 2. 【本物の重回帰用データ】X1 と X2 の2つを用意！
set.seed(42)
X1 <- runif(20, min = 1, max = 10)
X2 <- runif(20, min = 5, max = 15)

X1<-scale(X1)
X2<-scale(X2)

# 【正解の設定】 切片=5, X1の傾き=3, X2の傾き=-2, ノイズ=2
Y <- 5 + 3 * X1 - 2 * X2 + rnorm(20, mean = 0, sd = 2)

# 3. 記録用シート（4列に増やします）
# 1:切片, 2:X1の傾き, 3:X2の傾き, 4:ノイズ
samples <- matrix(0, nrow = N, ncol = 4)
colnames(samples) <- c("beta0", "beta1", "beta2", "sigma")

# スタート地点
samples[1, ] <- c(0, 0, 0, 10) 

# ==========================================
# 📐 「重回帰」のフィット感（対数尤度）を計算する関数
# ==========================================
calc_log_likelihood <- function(b0, b1, b2, sigma) {
  if (sigma <= 0) return(-Inf) 
  
  # ★修正ポイント：X1 と X2 の両方を使って予測する！
  Y_pred <- b0 + b1 * X1 + b2 * X2
  
  log_prob <- sum(dnorm(Y, mean = Y_pred, sd = sigma, log = TRUE))
  return(log_prob)
}

# ==========================================
# 🚀 スクラッチMHアルゴリズム（4次元の散歩）
# ==========================================
for (i in 2:N) {
  # いまの場所
  curr_b0    <- samples[i - 1, "beta0"]
  curr_b1    <- samples[i - 1, "beta1"]
  curr_b2    <- samples[i - 1, "beta2"] # ★追加
  curr_sigma <- samples[i - 1, "sigma"]
  
  # 1歩先の「次の候補地」を提案（b2もランダムに動かす）
  prop_b0    <- curr_b0    + rnorm(1, mean = 0, sd = 0.1)
  prop_b1    <- curr_b1    + rnorm(1, mean = 0, sd = 0.5)
  prop_b2    <- curr_b2    + rnorm(1, mean = 0, sd = 0.5) # ★追加
  prop_sigma <- curr_sigma + rnorm(1, mean = 0, sd = 0.05)
  
  # フィット感をそれぞれ計算（引数にb2を追加）
  curr_fit <- calc_log_likelihood(curr_b0, curr_b1, curr_b2, curr_sigma)
  prop_fit <- calc_log_likelihood(prop_b0, prop_b1, prop_b2, prop_sigma)
  
  # 移動確率 α の計算（対数なので引き算）
  log_alpha <- prop_fit - curr_fit
  alpha <- exp(log_alpha)
  
  # サイコロを振って移動の合否判定
  if (runif(1) < alpha) {
    samples[i, ] <- c(prop_b0, prop_b1, prop_b2, prop_sigma)
  } else {
    samples[i, ] <- c(curr_b0, curr_b1, curr_b2, curr_sigma)
  }
}

# 4. 前半の迷走期間を捨てる
clean_samples <- samples[10001:N, ]

# ==========================================
# 📊 答え合わせ
# ==========================================
cat("======================================\n")
cat(" 【正写真：スクラッチMCMC重回帰の結果】\n")
cat(" 切片  (正解は 5)  :", mean(clean_samples[, "beta0"]), "\n")
cat(" 傾き1 (正解は 3)  :", mean(clean_samples[, "beta1"]), "\n")
cat(" 傾き2 (正解は -2) :", mean(clean_samples[, "beta2"]), "\n") # ★追加
cat(" ノイズ(正解は 2)  :", mean(clean_samples[, "sigma"]), "\n")
cat("======================================\n")

#distribution of coefficients
#5. 【追加】4つのパラメーターの分布（ヒストグラム）を一画面で描画
par(mfrow = c(2, 2)) # 画面を2×2の4分割にする魔法のコマンド
for (name in colnames(clean_samples)) {
  hist(clean_samples[, name], main = paste(name, "の分布"), xlab = "値", col = "lightgreen", breaks = 35)
}
par(mfrow = c(1, 1))

# ==========================================
# 📊 追加：各サンプルデータにおける Y の予測分布（事後予測分布）
# ==========================================

# 1. 2万歩分の散歩データから、データ20個それぞれの「Yの予測値」をすべて計算する
# 行数が20000（MCMCの歩数）、列数が20（データ数）の行列を作ります
num_data <- length(Y)
num_mcmc <- nrow(clean_samples)
Y_pred_matrix <- matrix(0, nrow = num_mcmc, ncol = num_data)

for (m in 1:num_mcmc) {
  # そのステップでのパラメーターを取り出す
  b0 <- clean_samples[m, "beta0"]
  b1 <- clean_samples[m, "beta1"]
  b2 <- clean_samples[m, "beta2"]
  
  # 20個のデータ全体の予測値を計算して、行列のm行目に保存
  Y_pred_matrix[m, ] <- b0 + b1 * X1 + b2 * X2
}

# 2. 描画の設定
# 実際のYの値の範囲に合わせて、グラフの縦軸の限界（ylim）を自動調整
y_min <- min(c(Y, Y_pred_matrix)) - 2
y_max <- max(c(Y, Y_pred_matrix)) + 2

# 3. 20個のサンプルの予測分布を「箱ひげ図」で一列に並べる
boxplot(Y_pred_matrix, 
        main = "各サンプルにおける Y の事後予測分布", 
        xlab = "サンプル番号 (1 〜 20)", 
        ylab = "Y の値（予測の分布）", 
        col = "lightblue", 
        ylim = c(y_min, y_max),
        names = 1:num_data)

# ★超重要：MCMCの予測の箱の上に、「実際の生のYの値」を赤い丸でプロットする！
points(1:num_data, Y, col = "red", pch = 19, cex = 1.3)

# 凡例をつけておく
legend("topright", legend = c("実際のYの値", "MCMCの予測分布"), 
       col = c("red", "lightblue"), pch = c(19, 15), pt.cex = c(1.3, 2))# 📊 答え合わせ
# ==========================================
cat("======================================\n")
cat(" 【正写真：スクラッチMCMC重回帰の結果】\n")
cat(" 切片  (正解は 5)  :", mean(clean_samples[, "beta0"]), "\n")
cat(" 傾き1 (正解は 3)  :", mean(clean_samples[, "beta1"]), "\n")
cat(" 傾き2 (正解は -2) :", mean(clean_samples[, "beta2"]), "\n") # ★追加
cat(" ノイズ(正解は 2)  :", mean(clean_samples[, "sigma"]), "\n")
cat("======================================\n")


par(mfrow = c(2, 2)) # 描画エリアを2×2の4画面に分割

# 各パラメータのトレースプロットを描画[1.Caterpillar Plot]
plot(samples[, "beta0"], type = "l", main = "beta0 (Intercept)", xlab = "Iteration", ylab = "Value")
plot(samples[, "beta1"], type = "l", main = "beta1 (X1 slope)", xlab = "Iteration", ylab = "Value")
plot(samples[, "beta2"], type = "l", main = "beta2 (X2 slope)", xlab = "Iteration", ylab = "Value")
plot(samples[, "sigma"], type = "l", main = "sigma (Noise SD)", xlab = "Iteration", ylab = "Value")

# 例として、もう一つ別のスタート地点（初期値）から走らせてみる
samples2 <- matrix(0, nrow = N, ncol = 4)
colnames(samples2) <- c("beta0", "beta1", "beta2", "sigma")
samples2[1, ] <- c(10, 5, -5, 2) # 違う初期値からスタート


# バーンイン期間（例：最初の10000歩）を捨てる
if (!require(coda)) {install.packages("coda");library(coda)}
burn_in <- 10001:N
chain1 <- mcmc(samples[burn_in, ])
chain2 <- mcmc(samples2[burn_in, ])

# 2つのチェーンをまとめる
mcmc_chains <- mcmc.list(chain1, chain2)

# Gelman-Rubin の収束判定（R-hat）を計算[2,Gelman-Rubin]
gelman.diag(mcmc_chains)

#####
# ==========================================
# 🚀 2つ目のチェーン（初期値を変えてもう一度散歩する）
# ==========================================
# 違う初期値からスタート
samples2 <- matrix(0, nrow = N, ncol = 4)
colnames(samples2) <- c("beta0", "beta1", "beta2", "sigma")
samples2[1, ] <- c(10, 5, -5, 2) 

# ★抜けていた samples2 用の MCMC ループを実行！
for (i in 2:N) {
  # いまの場所
  curr_b0    <- samples2[i - 1, "beta0"]
  curr_b1    <- samples2[i - 1, "beta1"]
  curr_b2    <- samples2[i - 1, "beta2"]
  curr_sigma <- samples2[i - 1, "sigma"]
  
  # 1歩先の「次の候補地」を提案（歩幅はsamples1と同じ）
  prop_b0    <- curr_b0    + rnorm(1, mean = 0, sd = 0.1)
  prop_b1    <- curr_b1    + rnorm(1, mean = 0, sd = 0.3)
  prop_b2    <- curr_b2    + rnorm(1, mean = 0, sd = 0.3)
  prop_sigma <- curr_sigma + rnorm(1, mean = 0, sd = 0.05)
  
  # フィット感をそれぞれ計算
  curr_fit <- calc_log_likelihood(curr_b0, curr_b1, curr_b2, curr_sigma)
  prop_fit <- calc_log_likelihood(prop_b0, prop_b1, prop_b2, prop_sigma)
  
  # 移動確率 α の計算
  log_alpha <- prop_fit - curr_fit
  alpha <- exp(log_alpha)
  
  # サイコロを振って移動の合否判定
  if (runif(1) < alpha) {
    samples2[i, ] <- c(prop_b0, prop_b1, prop_b2, prop_sigma)
  } else {
    samples2[i, ] <- c(curr_b0, curr_b1, curr_b2, curr_sigma)
  }
}
#=====
# ==========================================
# 🎯 追加：現在の「採択率（Acceptance Rate）」を計算して表示
# ==========================================
# unique()関数で「値が変わった（移動した）回数」を数えて全体の歩数(N)で割ります
rate1 <- length(unique(samples[, "beta1"])) / N
rate2 <- length(unique(samples2[, "beta1"])) / N

cat("\n--- MCMC Acceptance Rate (目標: 23%〜50%) ---\n")
cat(" チェーン1の採択率 :", round(rate1 * 100, 1), "%\n")
cat(" チェーン2の採択率 :", round(rate2 * 100, 1), "%\n")
cat("--------------------------------------------\n")
#====

# ==========================================
# 📊 Gelman-Rubin の収束判定（R-hat）の計算
# ==========================================
if (!require(coda)) { install.packages("coda"); library(coda) }

# バーンイン期間（最初の10000歩）を捨てる
burn_in <- 25001:N #10001:N
chain1 <- mcmc(samples[burn_in, ])
chain2 <- mcmc(samples2[burn_in, ])

# 2つのチェーンをまとめる
mcmc_chains <- mcmc.list(chain1, chain2)

# R-hat を計算（これで正しい Point est. が出ます！）
gelman.diag(mcmc_chains)


# ==========================================
# 📈 視覚的チェック：2つのチェーンが合流したか確認
# ==========================================
par(mfrow = c(2, 2)) # 画面を4分割

# 2つの線を重ねてプロット（黒＝チェーン1, 赤＝チェーン2）
plot(samples[, "beta0"], type = "l", main = "beta0 (Intercept)", xlab = "Iteration", ylab = "Value")
lines(samples2[, "beta0"], col = "red")

plot(samples[, "beta1"], type = "l", main = "beta1 (X1 slope)", xlab = "Iteration", ylab = "Value")
lines(samples2[, "beta1"], col = "red")

plot(samples[, "beta2"], type = "l", main = "beta2 (X2 slope)", xlab = "Iteration", ylab = "Value")
lines(samples2[, "beta2"], col = "red")

plot(samples[, "sigma"], type = "l", main = "sigma (Noise SD)", xlab = "Iteration", ylab = "Value")
lines(samples2[, "sigma"], col = "red")

par(mfrow = c(1, 1)) # 画面分割を戻す