Rで比推定

# --- 図1: 母集団の散布図 ---
df_pop <- data.frame(x = x_pop, y = y_pop)
p1 <- ggplot(df_pop[sample(N, 300), ], aes(x = x, y = y)) +
  geom_point(alpha = 0.4, color = "#2c7bb6", size = 1.5) +
  geom_abline(
    slope = R_pop, intercept = 0,
    color = "#d7191c", linewidth = 1, linetype = "dashed"
  ) +
  labs(
    title = "母集団：x と y の関係",
    subtitle = sprintf("真の比 R = %.3f（赤点線：y = R·x）", R_pop),
    x = "補助変数 x（昨年収穫量 kg）",
    y = "目標変数 y（今年収穫量 kg）"
  ) +
  theme_bw() +
  theme(plot.title = element_text(size = 13, face = "bold"), plot.subtitle = element_text(size = 9))

# --- 図2: 推定量の分布比較（ヒストグラム）---
df_long <- df_sim |>
  select(simple, ratio) |>
  pivot_longer(everything(),
    names_to = "method",
    values_to = "estimate"
  ) |>
  mutate(method = recode(method,
    simple = "単純推定",
    ratio  = "比推定"
  ))

p2 <- ggplot(df_long, aes(x = estimate, fill = method)) +
  geom_histogram(aes(y = after_stat(density)),
    bins = 60, alpha = 0.6, position = "identity"
  ) +
  geom_vline(
    xintercept = mu_y, color = "black",
    linewidth = 1.2, linetype = "solid"
  ) +
  scale_fill_manual(values = c(
    "単純推定" = "#fc8d59",
    "比推定" = "#2c7bb6"
  )) +
  labs(
    title = "推定量の標本分布",
    subtitle = sprintf("黒線：真の母平均 μ_y = %.2f", mu_y),
    x = "推定値", y = "密度", fill = "推定法"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(size = 13, face = "bold"), plot.subtitle = element_text(size = 9),
    legend.position = "top"
  )

# --- 図3: 相対効率 vs 相関係数 ---
p3 <- ggplot(df_eff, aes(x = rho, y = efficiency)) +
  geom_line(color = "#1a9641", linewidth = 1.2) +
  geom_hline(
    yintercept = 1, linetype = "dashed",
    color = "gray40", linewidth = 0.8
  ) +
  geom_vline(
    xintercept = threshold_rho,
    color = "#d7191c", linetype = "dotted", linewidth = 1
  ) +
  geom_vline(
    xintercept = cor(x_pop, y_pop),
    color = "#2c7bb6", linetype = "solid", linewidth = 1
  ) +
  annotate("text",
    x = threshold_rho - 0.35, y = 2,
    label = sprintf("閾値\nρ=%.2f", threshold_rho),
    color = "#d7191c", size = 3.0, hjust = 0
  ) +
  annotate("text",
    x = cor(x_pop, y_pop) - 0.35,
    y = 0.4,
    label = sprintf("今回\nρ=%.2f", cor(x_pop, y_pop)),
    color = "#2c7bb6", size = 3.0, hjust = 0
  ) +
  coord_cartesian(ylim = c(0, 5)) +
  labs(
    title = "相関係数と比推定の相対効率",
    subtitle = "相対効率 = Var(単純推定) / Var(比推定)  ＞1 なら比推定が有利",
    x = "相関係数 ρ(x, y)",
    y = "相対効率"
  ) +
  theme_bw() +
  theme(plot.title = element_text(size = 13, face = "bold"), plot.subtitle = element_text(size = 9))

# --- 図4: 標本サイズと推定精度（RMSE）の比較 ---
n_seq <- seq(10, 600, by = 10)
rmse_df <- lapply(n_seq, function(ni) {
  sims <- replicate(2000, {
    idx <- sample(N, ni, replace = FALSE)
    xs <- x_pop[idx]
    ys <- y_pop[idx]
    c(
      simple = mean(ys),
      ratio = mean(ys) / mean(xs) * mu_x
    )
  })
  data.frame(
    n = ni,
    method = c("単純推定", "比推定"),
    rmse = c(
      sqrt(mean((sims["simple", ] - mu_y)^2)),
      sqrt(mean((sims["ratio", ] - mu_y)^2))
    )
  )
}) |> bind_rows()

p4 <- ggplot(rmse_df, aes(x = n, y = rmse, color = method)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 1.5) +
  scale_color_manual(values = c(
    "単純推定" = "#fc8d59",
    "比推定" = "#2c7bb6"
  )) +
  labs(
    title = "標本サイズと推定精度（RMSE）",
    subtitle = "標本サイズが小さいほど比推定の優位性が顕著",
    x = "標本サイズ n",
    y = "RMSE（二乗平均誤差の平方根）",
    color = "推定法"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(size = 13, face = "bold"), plot.subtitle = element_text(size = 9),
    legend.position = "top"
  )

# --- 統合出力 ---
combined <- (p1 + p2) / (p3 + p4) +
  plot_annotation(
    title = "比推定（Ratio Estimation）シミュレーション",
    subtitle = sprintf(
      "母集団サイズ N=%d　標本サイズ n=%d　シミュレーション回数 B=%d",
      N, n, B
    ),
    theme = theme(
      plot.title    = element_text(size = 16, face = "bold"),
      plot.subtitle = element_text(size = 11, color = "gray40")
    )
  )

combined