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

# ターゲットにしたい複雑な関数（ガンマ分布の形）
target_pdf <- function(x) {
  if (x < 0) return(0)
  return(x^2 * exp(-x)) # 数式で解くのが面倒な形
}

# 2. 散歩の記録用ベクトルを用意
x_samples <- numeric(N)
x_samples[1] <- 1.0 # スタート地点

# 3. MHアルゴリズムで散歩スタート！
set.seed(42)
for (i in 2:N) {
  current_x <- x_samples[i - 1]
  
  # いまの場所の近く（標準正規分布）から「次の候補地」を提案する
  proposed_x <- current_x + rnorm(1, mean = 0, sd = 0.5)
  
  # 高さを比較するための「受容確率 α (アルファ)」を計算
  # 次の高さ / いまの高さ
  alpha <- target_pdf(proposed_x) / target_pdf(current_x)
  
  # サイコロを振って移動するか判定
  if (runif(1) < alpha) {
    x_samples[i] <- proposed_x # 合格！新しい場所に移動
  } else {
    x_samples[i] <- current_x  # 不合格！元の場所にとどまる
  }
}

# 4. 散歩の軌跡が、ちゃんと関数の形になっているか確認（最初の数歩を捨てるバーンイン）
clean_samples <- x_samples[1001:N]

# グラフで答え合わせ
hist(clean_samples, breaks = 50, probability = TRUE, 
     main = "MHアルゴリズムでハッキングした乱数", xlab = "X", col = "skyblue")
curve(x^2 * exp(-x) / 2, add = TRUE, col = "red", lwd = 3) # 理論上の曲線