Rでピアソンの経験則

# ============================================================
# ピアソンの経験則シミュレーション：5分布の可視化
# ggplot2 による密度プロット
# ============================================================

library(ggplot2)
library(dplyr)
library(tidyr)

# ------------------------------------------------------------
# データ生成 & 統計量計算
# ------------------------------------------------------------
samples_list <- list()
vlines_list <- list()

for (d in dists) {
  x <- d$gen()
  mn <- mean(x)
  md <- median(x)
  mo <- d$mode
  pred <- mn - (1 / 3) * (mn - mo)
  sk <- mean((x - mn)^3) / sd(x)^3
  rel <- (pred - md) / md * 100

  label <- sprintf(
    "歪度: %.3f\n相対誤差: %+.2f%%",
    sk, rel
  )

  samples_list[[d$name]] <- data.frame(x = x, dist = d$name)

  vlines_list[[d$name]] <- data.frame(
    dist = d$name,
    xintercept = c(mn, md, mo, pred),
    stat = factor(
      c("平均値", "実際の中央値", "最頻値", "ピアソン予測"),
      levels = c("平均値", "最頻値", "実際の中央値", "ピアソン予測")
    ),
    label = label
  )
}

df_samples <- bind_rows(samples_list)
df_vlines <- bind_rows(vlines_list)

# facet の表示順を固定
dist_order <- sapply(dists, `[[`, "name")
df_samples$dist <- factor(df_samples$dist, levels = dist_order)
df_vlines$dist <- factor(df_vlines$dist, levels = dist_order)

# ラベル用（各パネル右上に歪度・相対誤差を表示）
df_labels <- df_vlines |>
  group_by(dist) |>
  slice(1) |>
  ungroup() |>
  select(dist, label)

# ------------------------------------------------------------
# 配色・線種・線幅の定義
# ------------------------------------------------------------
line_colors <- c(
  "平均値" = "orange",
  "最頻値" = "green",
  "実際の中央値" = "blue",
  "ピアソン予測" = "red"
)
line_types <- c(
  "平均値" = "dashed",
  "最頻値" = "dashed",
  "実際の中央値" = "solid",
  "ピアソン予測" = "dotted"
)
line_widths <- c(
  "平均値" = 1.0,
  "最頻値" = 1.0,
  "実際の中央値" = 1.2,
  "ピアソン予測" = 1.2
)

# ------------------------------------------------------------
# line_types・line_widths から凡例用ベクトルを動的生成
# override.aes の順序は factor levels に合わせる
# ------------------------------------------------------------
legend_levels <- levels(df_vlines$stat)
legend_linetypes <- unname(line_types[legend_levels])
legend_linewidths <- unname(line_widths[legend_levels])

# geom_vline を linewidth の値ごとにグループ化して描画
unique_widths <- unique(line_widths)
geom_vline_layers <- lapply(unique_widths, function(w) {
  stats_w <- names(line_widths)[line_widths == w]
  geom_vline(
    data = df_vlines[df_vlines$stat %in% stats_w, ],
    aes(xintercept = xintercept, color = stat, linetype = stat),
    linewidth = w
  )
})

# ------------------------------------------------------------
# プロット
# ------------------------------------------------------------
p <- ggplot(df_samples, aes(x = x)) +

  # 密度曲線
  geom_density(
    fill      = "#2a78d6",
    alpha     = 0.12,
    color     = "#2a78d6",
    linewidth = 0.5
  ) +

  # 縦線（linewidth ごとに分割したレイヤーを追加）
  geom_vline_layers +

  # 歪度・相対誤差のラベル（右上）
  geom_text(
    data = df_labels,
    aes(label = label),
    x = Inf, y = Inf,
    hjust = 1.08, vjust = 1.5,
    size = 3, color = "#52514e",
    inherit.aes = FALSE,
    lineheight = 1.3
  ) +

  # スケール（override.aes を line_types・line_widths から動的参照）
  scale_color_manual(
    name = NULL,
    values = line_colors,
    guide = guide_legend(
      override.aes = list(
        linetype  = legend_linetypes,
        linewidth = legend_linewidths
      )
    )
  ) +
  scale_linetype_manual(name = NULL, values = line_types, guide = "none") +

  # ファセット（分布ごとに独立した x 軸）
  facet_wrap(~dist, scales = "free", ncol = 2) +

  # 軸ラベル
  labs(
    title    = "ピアソンの経験則：5分布における中央値の予測",
    subtitle = "密度曲線と各統計量の位置（n = 100,000）",
    x        = "値",
    y        = "密度",
    caption  = "ピアソン予測：中央値 ≈ 平均値 − (1/3)(平均値 − 最頻値)"
  ) +

  # テーマ
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(
      size = 14, face = "bold",
      margin = margin(b = 4)
    ),
    plot.subtitle = element_text(
      size = 11, color = "#52514e",
      margin = margin(b = 12)
    ),
    plot.caption = element_text(
      size = 9, color = "#52514e",
      hjust = 0
    ),
    strip.text = element_text(
      size = 10, face = "bold",
      lineheight = 1.2
    ),
    legend.position = "bottom",
    legend.text = element_text(size = 10),
    legend.key.width = unit(1.8, "cm"),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "#e8e7e1"),
    axis.text = element_text(size = 9, color = "#52514e"),
    axis.title = element_text(size = 10, color = "#52514e"),
    plot.margin = margin(12, 16, 12, 16)
  )

p