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 = "#756bb1", size = 1.5) +
  geom_smooth(
    method = "lm", se = FALSE,
    color = "#d7191c", linewidth = 1, linetype = "dashed"
  ) +
  labs(
    title = "母集団：x と y の関係",
    subtitle = sprintf("相関係数 ρ = %.3f（負の相関）", rho),
    x = "補助変数 x（気温偏差）",
    y = "目標変数 y（害虫被害量）"
  ) +
  theme_bw() +
  theme(plot.title = element_text(size = 13, face = "bold"), plot.subtitle = element_text(size = 10))

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

p2 <- ggplot(df_long, aes(x = estimate, fill = method)) +
  geom_histogram(aes(y = after_stat(density)),
    bins = 60, alpha = 0.55, position = "identity"
  ) +
  geom_vline(
    xintercept = mu_y, color = "black",
    linewidth = 1.2, linetype = "solid"
  ) +
  scale_fill_manual(values = c(
    "単純推定" = "#fc8d59",
    "比推定" = "#2c7bb6",
    "積推定" = "#1a9641"
  )) +
  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 = 10),
    legend.position = "top"
  )

# --- 図3: 相対効率 vs 相関係数 ---
p3 <- ggplot(df_eff, aes(x = rho, y = efficiency, color = method)) +
  geom_line(linewidth = 1.2) +
  geom_hline(
    yintercept = 1, linetype = "dashed",
    color = "gray40", linewidth = 0.8
  ) +
  geom_vline(
    xintercept = threshold_ratio,
    color = "#2c7bb6", linetype = "dotted", linewidth = 1
  ) +
  geom_vline(
    xintercept = threshold_product,
    color = "#1a9641", linetype = "dotted", linewidth = 1
  ) +
  geom_vline(
    xintercept = rho,
    color = "gray20", linetype = "solid", linewidth = 0.9
  ) +
  annotate("text",
    x = threshold_ratio - 0.03, y = 4.5,
    label = sprintf("比推定\n閾値 ρ=%.2f", threshold_ratio),
    color = "#2c7bb6", size = 3.2, hjust = 1
  ) +
  annotate("text",
    x = threshold_product + 0.03, y = 4.5,
    label = sprintf("積推定\n閾値 ρ=%.2f", threshold_product),
    color = "#1a9641", size = 3.2, hjust = 0
  ) +
  annotate("text",
    x = rho + 0.03, y = 2.5,
    label = sprintf("今回\nρ=%.2f", rho),
    color = "gray20", size = 3.2, hjust = 0
  ) +
  scale_color_manual(values = c(
    "比推定" = "#2c7bb6",
    "積推定" = "#1a9641"
  )) +
  coord_cartesian(ylim = c(0, 6)) +
  labs(
    title = "相関係数と相対効率（理論曲線）",
    subtitle = "相対効率 = Var(単純推定) / Var(各推定法)　＞1 なら\n単純推定より有利",
    x = "相関係数 ρ(x, y)",
    y = "相対効率",
    color = "推定法"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(size = 13, face = "bold"), plot.subtitle = element_text(size = 10),
    legend.position = "top"
  )

# --- 図4: 標本サイズと RMSE ---
n_seq <- seq(10, 200, 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,
      product = 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)),
      sqrt(mean((sims["product", ] - 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",
    "積推定" = "#1a9641"
  )) +
  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 = 10),
    legend.position = "top"
  )

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

combined