Rで確率分布:ラプラス分布

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

# ラプラス分布の確率密度関数(PDF)を定義
d_laplace <- function(x, mu, b) {
  1 / (2 * b) * exp(-abs(x - mu) / b)
}

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

# 2. 異なるパラメータを持つ分布の確率密度を計算
# ラプラス分布(μ=0, b=1)の分散は 2*b^2 = 2。
# よって比較対象の正規分布の標準偏差は sqrt(2)
b1_var <- 2 * 1^2
norm_sd <- sqrt(b1_var)

df <- tibble(
  x = x_vals
) %>%
  mutate(
    `ラプラス (μ=0, b=1)` = d_laplace(x, mu = 0, b = 1),
    `ラプラス (μ=5, b=1)` = d_laplace(x, mu = 5, b = 1),
    `ラプラス (μ=0, b=2)` = d_laplace(x, mu = 0, b = 2),
    `正規分布 (μ=0, σ=√2)` = dnorm(x, mean = 0, sd = norm_sd)
  )

# 3. ggplotで描画しやすいように、データを「ロングフォーマット」に変換
df_long <- df %>%
  pivot_longer(
    cols = -x,
    names_to = "distribution",
    values_to = "density"
  ) %>%
  # 凡例の順序を調整
  mutate(distribution = factor(distribution, levels = c(
    "ラプラス (μ=0, b=1)",
    "ラプラス (μ=5, b=1)",
    "ラプラス (μ=0, b=2)",
    "正規分布 (μ=0, σ=√2)"
  )))

# 4. 各分布に割り当てる色と線種を定義
manual_colors <- c(
  `ラプラス (μ=0, b=1)` = "blue",
  `ラプラス (μ=5, b=1)` = "red",
  `ラプラス (μ=0, b=2)` = "darkgreen",
  `正規分布 (μ=0, σ=√2)` = "black"
)

manual_linetypes <- c(
  `ラプラス (μ=0, b=1)` = "solid",
  `ラプラス (μ=5, b=1)` = "solid",
  `ラプラス (μ=0, b=2)` = "solid",
  `正規分布 (μ=0, σ=√2)` = "dashed"
)

# 5. ggplotを使用してチャートを描画
p <- ggplot(df_long, aes(x = x, y = density, color = distribution, linetype = distribution)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = manual_colors) +
  scale_linetype_manual(values = manual_linetypes) +
  labs(
    title = "ラプラス分布の変化と正規分布との比較",
    subtitle = "正規分布に比べピークが鋭く、裾が重いのが特徴",
    x = "xの値",
    y = "確率密度",
    color = "分布の種類",
    linetype = "分布の種類"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "gray40")
  ) +
  guides(color = guide_legend(nrow = 2))

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

p_log <- ggplot(df_long, aes(x = x, y = density, color = distribution, linetype = distribution)) +
  geom_line(linewidth = 1.2) +
  # 対数スケールに変更するレイヤーを追加。ただしdensity=0はプロットできないので、非常に小さい値をクリップ
  scale_y_log10(
    breaks = scales::trans_breaks("log10", function(x) 10^x),
    labels = scales::trans_format("log10", scales::math_format(10^.x))
  ) +
  scale_color_manual(values = manual_colors) +
  scale_linetype_manual(values = manual_linetypes) +
  labs(
    title = "ラプラス分布と正規分布の比較（y軸対数スケール）",
    subtitle = "対数スケールにしてラプラス分布の裾の重さ（減少の遅さ）を明確化",
    x = "xの値",
    y = "確率密度 (対数スケール)",
    color = "分布の種類",
    linetype = "分布の種類"
  ) +
  # y軸の表示範囲を指定。下限を10^-10に設定
  coord_cartesian(xlim = c(-8, 8), ylim = c(1e-10, NA)) + # NAは上限を自動設定
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "gray40")
  ) +
  guides(color = guide_legend(nrow = 2))

# 対数スケールのチャートを表示
print(p_log)