Rで確率分布:レヴィ分布

# 必要なライブラリを読み込みます
library(ggplot2)
library(dplyr)
library(tidyr)
library(rmutil) # 関数dlevyを利用

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

# 2. 異なるパラメータを持つレヴィ分布の確率密度を計算
df <- tibble(
  x = x_vals
) %>%
  mutate(
    # μ=0, c=1 の場合 (m=0)
    `μ=0, c=1` = {
      # 1. 結果を格納するベクトルを0で初期化
      res <- numeric(length(x))
      # 2. 条件を満たすデータの位置を特定
      idx <- x > 0
      # 3. 条件を満たすデータだけで計算し、結果を代入
      res[idx] <- dlevy(y = x[idx], m = 0, s = 1)
      # 最終的な結果ベクトルを返す
      res
    },

    # μ=2, c=1 の場合 (m=2)
    `μ=2, c=1` = {
      res <- numeric(length(x))
      idx <- x > 2
      # x[idx]が空でない場合のみ計算する
      if (any(idx)) {
        res[idx] <- dlevy(y = x[idx], m = 2, s = 1)
      }
      res
    },

    # μ=0, c=2 の場合 (m=0)
    `μ=0, c=2` = {
      res <- numeric(length(x))
      idx <- x > 0
      res[idx] <- dlevy(y = x[idx], m = 0, s = 2)
      res
    }
  )

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

# 4. 各分布に割り当てる色を定義
manual_colors <- c(
  `μ=0, c=1` = "blue",
  `μ=2, c=1` = "red",
  `μ=0, c=2` = "darkgreen"
)

# 5. ggplotを使用してチャートを描画
p <- ggplot(df_long, aes(x = x, y = density, color = parameters)) +
  # 密度がNAになる部分（定義域外）はプロットしない
  geom_line(data = . %>% filter(!is.na(density) & density > 1e-10), linewidth = 1.2) + # y軸0対応
  scale_color_manual(values = manual_colors) +
  labs(
    title = "レヴィ分布の変化",
    subtitle = "位置(μ)は分布の下限、尺度(c)は広がりを決定する",
    x = "xの値",
    y = "確率密度",
    color = "パラメータ (μ, c)"
  ) +
  coord_cartesian(xlim = c(0, 10), ylim = c(0, 0.5)) + # 見やすい範囲に調整
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "gray40")
  ) +
  scale_x_continuous(breaks = scales::pretty_breaks(10))

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