Rで標準偏差と標準誤差の違い

# 必要なライブラリを読み込む
library(ggplot2)
library(dplyr)

# 乱数のシードを固定
seed <- 20250802
set.seed(seed)

# 1. 母集団の定義
# 平均50, 標準偏差10の正規分布に従う10万個のデータ
population_mean <- 50
population_sd <- 10
population_data <- rnorm(100000, mean = population_mean, sd = population_sd)

# 母集団の実際の平均と標準偏差を確認
cat("--- 母集団の情報 ---\n\n")
cat(paste0("母集団の真の平均: ", population_mean, "\n"))
cat(paste0("母集団の真の標準偏差: ", population_sd, "\n\n"))

# 2. シミュレーションの設定
n_sim <- 1000 # シミュレーション（サンプリング）の回数
n_small <- 10 # 小さいサンプルサイズ
n_large <- 100 # 大きいサンプルサイズ

# 3. シミュレーションの実行
# 小さいサンプルサイズ (n=10) でのシミュレーション
# replicate関数を使い、1000回サンプリングを繰り返す
results_small <- replicate(n_sim, {
  # 母集団からn_small個のデータをランダムにサンプリング
  sample_data <- sample(population_data, n_small)
  # サンプルの平均と標準偏差を計算（sd()は分母n-1で計算される）
  c(mean = mean(sample_data), sd = sd(sample_data))
})

# 大きいサンプルサイズ (n=100) でのシミュレーション
results_large <- replicate(n_sim, {
  # 母集団からn_large個のデータをランダムにサンプリング
  sample_data <- sample(population_data, n_large)
  # サンプルの平均と標準偏差を計算（sd()は分母n-1で計算される）
  c(mean = mean(sample_data), sd = sd(sample_data))
})

# 結果をデータフレームに変換
# t()で行と列を入れ替える
df_small <- as.data.frame(t(results_small))
df_large <- as.data.frame(t(results_large))


# 4. 結果の分析と出力
cat("--- シミュレーション結果の分析 ---\n\n")

# << 標準偏差の比較 >>
# 各シミュレーションで得られた1000個の「標本標準偏差」の平均値
mean_sd_small <- mean(df_small$sd)
mean_sd_large <- mean(df_large$sd)

cat("<< 標準偏差 (個々のサンプルのばらつき) の比較 >>\n\n")
cat(paste0("サンプルサイズ n=", n_small, " の場合、標本標準偏差の平均値は ", round(mean_sd_small, 3), " です。\n"))
cat(paste0("サンプルサイズ n=", n_large, " の場合、標本標準偏差の平均値は ", round(mean_sd_large, 3), " です。\n"))
cat(paste0("→ サンプルサイズが大きくなっても、標本標準偏差の値自体は母集団の標準偏差(", population_sd, ")から大きくは変わりません。\n\n"))

# << 標準誤差の比較 >>
# ここでは統計量として「標本平均」を扱っています
# (A) シミュレーションから直接計算した標準誤差（＝1000個の標本平均の標準偏差）
se_sim_small <- sd(df_small$mean)
se_sim_large <- sd(df_large$mean)

# (B) 理論式から計算した標準誤差 (s / sqrt(n))
# 1000個のサンプルから計算した標準誤差の平均値
se_formula_small <- mean(df_small$sd / sqrt(n_small))
se_formula_large <- mean(df_large$sd / sqrt(n_large))


cat("<< 標本平均の標準誤差 (標本平均のばらつき) の比較 >>\n\n")
cat(paste0("サンプルサイズ n=", n_small, " の場合:\n"))
cat(paste0("  (A) 標本平均の標準偏差（シミュレーションによるSE）: ", round(se_sim_small, 3), "\n"))
cat(paste0("  (B) 式(s/√n)から求めたSEの平均値　　　　　　: ", round(se_formula_small, 3), "\n\n"))

cat(paste0("サンプルサイズ n=", n_large, " の場合:\n"))
cat(paste0("  (A) 標本平均の標準偏差（シミュレーションによるSE）: ", round(se_sim_large, 3), "\n"))
cat(paste0("  (B) 式(s/√n)から求めたSEの平均値　　　　　　: ", round(se_formula_large, 3), "\n\n"))

cat("→ サンプルサイズが大きくなると、標本平均のばらつき（標準誤差）は明らかに小さくなります。\n")
cat("→ サンプルサイズの大小に関わらず、シミュレーションで得られた標本平均のばらつき(A)は、理論式(B)とほぼ一致します。\n\n")


# 5. 結果のプロット
# プロット用のデータフレームを作成
plot_data <- bind_rows(
  mutate(df_small, sample_size = paste0("n = ", n_small)),
  mutate(df_large, sample_size = paste0("n = ", n_large))
) %>%
  # 凡例の順序をnが小さい順にする
  mutate(sample_size = factor(sample_size, levels = c(paste0("n = ", n_small), paste0("n = ", n_large))))

# ggplotで標本平均の分布を可視化
g <- ggplot(plot_data, aes(x = mean, fill = sample_size)) +
  # 密度プロットを半透明で描画
  geom_density(alpha = 0.6, color = "black") +

  # 母集団の真の平均を黒の破線で示す
  geom_vline(
    xintercept = population_mean,
    color = "black",
    linetype = "dashed",
    linewidth = 1
  ) +

  # 凡例とラベルの設定
  labs(
    title = "標本平均の分布とサンプルサイズ",
    subtitle = "サンプルサイズが大きくなると、標本平均のばらつき（＝標準誤差）が小さくなる様子",
    x = "標本平均",
    y = "密度",
    fill = "サンプルサイズ"
  ) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "bottom")

# プロットを表示
print(g)