Rで確率分布:グンベル分布

# 必要なライブラリを読み込みます
library(evd)
library(ggplot2)
library(dplyr)
library(tidyr)

# 1. 描画範囲となるx軸の値を生成
x_vals <- seq(-5, 15, length.out = 1000)

# 2. 異なるパラメータを持つグンベル分布の確率密度を計算
# dgumbel(x, loc, scale) を使用
df <- tibble(
  x = x_vals
) %>%
  mutate(
    `μ=0, σ=1` = dgumbel(x, loc = 0, scale = 1),
    `μ=2, σ=1` = dgumbel(x, loc = 2, scale = 1),
    `μ=0, σ=2` = dgumbel(x, loc = 0, scale = 2)
  )

# 3. ggplotで描画しやすいように、データを「ロングフォーマット」に変換
df_long <- df %>%
  pivot_longer(
    cols = -x,
    names_to = "parameters",
    values_to = "density"
  ) %>%
  # 凡例の順序を調整
  mutate(parameters = factor(parameters, levels = c(
    "μ=0, σ=1",
    "μ=2, σ=1",
    "μ=0, σ=2"
  )))

# 4. 各分布に割り当てる色を定義
manual_colors <- c(
  `μ=0, σ=1` = "blue",
  `μ=2, σ=1` = "red",
  `μ=0, σ=2` = "darkgreen"
)

# 5. ggplotを使用してチャートを描画
p <- ggplot(df_long, aes(x = x, y = density, color = parameters)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = manual_colors) +
  labs(
    title = "パラメータ(μ, σ)によるグンベル分布の変化",
    subtitle = "最大値の分布をモデル化するのに適した、右に歪んだ形状を持つ",
    x = "xの値",
    y = "確率密度",
    color = "パラメータ (μ, σ)"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "gray40")
  )

# チャートの表示
print(p)

# シミュレーションの設定
seed <- 20250808
set.seed(seed)
num_simulations <- 10000
sample_sizes <- c(10, 100, 10000)

# シミュレーションの実行
simulation_results <- lapply(sample_sizes, function(n) {
  max_values <- replicate(num_simulations, max(rnorm(n)))
  tibble(
    n = paste0("n = ", n),
    max_value = max_values
  )
}) %>%
  bind_rows()

# 理論的なグンベル分布のパラメータを計算
gumbel_params <- tibble(
  n_val = sample_sizes
) %>%
  mutate(
    c_n = sqrt(2 * log(n_val)),
    d_n = c_n - (log(log(n_val)) + log(4 * pi)) / (2 * c_n),
    loc = d_n,
    scale = 1 / c_n,
    n = paste0("n = ", n_val)
  )

# 可視化
p_simulation <- ggplot(simulation_results, aes(x = max_value)) +
  # ヒストグラムを描画
  geom_histogram(aes(y = after_stat(density)), bins = 50, fill = "lightblue", color = "black") +

  # 理論的なグンベル分布の曲線を描画
  # 別のデータフレーム(gumbel_params)を使うため、geom_line() を使用
  geom_line(
    data = . %>% distinct(n) %>% left_join(gumbel_params, by = "n") %>%
      rowwise() %>%
      mutate(x_range = list(seq(min(simulation_results$max_value[simulation_results$n == n]),
        max(simulation_results$max_value[simulation_results$n == n]),
        length.out = 200
      ))) %>%
      unnest(x_range) %>%
      mutate(density = dgumbel(x_range, loc = loc, scale = scale)),
    aes(x = x_range, y = density, color = "理論値 (グンベル分布)"),
    linewidth = 1.2,
    # ggplot() の aes(x = max_value) を継承しない
    inherit.aes = FALSE
  ) +

  # ファセットで分割
  facet_wrap(~n, scales = "free") +
  scale_color_manual(values = c("理論値 (グンベル分布)" = "red")) +
  labs(
    title = "フィッシャー・ティペット・グネデンコの定理のシミュレーション",
    subtitle = "正規分布からの最大値の分布は、サンプルサイズ(n)が大きいほどグンベル分布に近づく",
    x = "各サンプルセットの最大値",
    y = "確率密度",
    color = ""
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "gray40", size = 12),
    strip.text = element_text(face = "bold", size = 12)
  )

# チャートの表示
print(p_simulation)