Rで確率分布:フレシェ分布

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

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

# 2. 異なるパラメータを持つフレシェ分布の確率密度を計算
# dfrechet(x, loc, scale, shape) を使用
df <- tibble(
  x = x_vals
) %>%
  mutate(
    `α=1 (平均なし)` = dfrechet(x, loc = 0, scale = 1, shape = 1),
    `α=2 (分散なし)` = dfrechet(x, loc = 0, scale = 1, shape = 2),
    `α=5 (分散あり)` = dfrechet(x, loc = 0, scale = 1, shape = 5)
  )

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

# 4. 各分布に割り当てる色を定義
manual_colors <- c(
  `α=1 (平均なし)` = "blue",
  `α=2 (分散なし)` = "red",
  `α=5 (分散あり)` = "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 = "形状パラメータ (α)"
  ) +
  coord_cartesian(xlim = c(0, 8)) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "gray40")
  )

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

# 必要なパッケージ
library(patchwork)
library(actuar)

# シミュレーションの設定
seed <- 20250808
set.seed(seed)
num_simulations <- 10000
sample_sizes <- c(10, 100, 10000)
pareto_alpha <- 3.5 # 元のパレート分布の形状パラメータ

# シミュレーションの実行
simulation_results_fr <- lapply(sample_sizes, function(n) {
  # パレート分布(shape=3.5, scale=1)からn個サンプリングし、最大値を取得
  max_values <- replicate(num_simulations, max(rpareto(n, shape = pareto_alpha, scale = 1)))
  tibble(
    n_val = n,
    max_value = max_values
  )
}) %>%
  bind_rows()

# 理論的なフレシェ分布のパラメータを計算
# 元の分布がパレート(xm, α_p)の場合、最大値の分布は
# Frechet(shape = α_p, scale = xm*(n)^(1/α_p), loc = 0) で近似できる
frechet_params <- tibble(
  n_val = sample_sizes
) %>%
  mutate(
    shape = pareto_alpha,
    scale = (n_val)^(1 / pareto_alpha),
    loc = 0 # loc=0ではなく、loc=scale となる場合もあるが、ここでは単純化
  )

# 可視化

# 各サンプルサイズ（ファセット）に対応するプロットを個別に生成する関数
create_plot_for_n <- function(sample_n, sim_data, param_data) {
  # 対象となるnのデータとパラメータを抽出
  data_n <- sim_data %>% filter(n_val == sample_n)
  params_n <- param_data %>% filter(n_val == sample_n)

  # 描画範囲を決定 (99.9%タイル)して巨大な値を描画から除外。あくまでも描画のため。
  plot_range <- quantile(data_n$max_value, c(0, 0.999))

  # 個別のggplotオブジェクトを作成
  p <- ggplot(data_n, aes(x = max_value)) +
    geom_histogram(aes(y = after_stat(density)), bins = 50, fill = "lightgreen", color = "black") +
    stat_function(
      fun = dfrechet,
      args = list(loc = params_n$loc, scale = params_n$scale, shape = params_n$shape),
      aes(color = "理論値 (フレシェ分布)"),
      linewidth = 1.2,
      n = 5001
    ) +
    # 描画範囲を個別に設定
    coord_cartesian(xlim = plot_range) +
    scale_color_manual(values = c("理論値 (フレシェ分布)" = "purple")) +
    labs(title = paste("n =", sample_n)) +
    theme_minimal(base_size = 12) +
    theme(legend.position = "none") # 凡例は結合後に一つだけ表示

  return(p)
}

# 各サンプルサイズでプロットを生成
plot_n10 <- create_plot_for_n(sample_sizes[1], simulation_results_fr, frechet_params)
plot_n100 <- create_plot_for_n(sample_sizes[2], simulation_results_fr, frechet_params)
plot_n10000 <- create_plot_for_n(sample_sizes[3], simulation_results_fr, frechet_params)

# patchworkを使ってプロットを結合
p_simulation_fr <- (plot_n10 / plot_n100 / plot_n10000) +
  plot_layout(guides = "collect") + # 凡例をまとめる
  plot_annotation(
    title = "フィッシャー・ティペット・グネデンコの定理のシミュレーション (フレシェ)",
    subtitle = "パレート分布(重い裾)からの最大値の分布は、サンプルサイズ(n)が大きいほどフレシェ分布に近づく",
    theme = theme(
      plot.title = element_text(face = "bold", size = 16),
      plot.subtitle = element_text(color = "gray40", size = 12),
      legend.position = "top"
    )
  )

# 最終的なチャートの表示
print(p_simulation_fr)