Rで確率分布:ワイブル分布

Rで 確率分布:ワイブル分布 を試みます。

本ポストはこちらの続きです。

Rで確率分布:指数分布
const typesetMath = (el) => { if (window.MathJax) { // MathJax Typeset window.MathJax.typeset(); } else if (window.katex...
www.saecanet.com
2025.07.31

1. ワイブル分布とは

ワイブル分布（Weibull Distribution）は、特に製品の寿命や故障時間の分析において広く用いられる連続確率分布です。指数分布が「無記憶性」を持つ、つまり故障率が常に一定であるという強い仮定を置くのに対し、ワイブル分布は時間の経過とともに故障率が変化する（増加する、減少する、または一定）ような、より柔軟な状況をモデル化することができます。

確率密度関数 (Probability Density Function, PDF)

ワイブル分布に従う確率変数 \(X\) の確率密度関数 \(f(x)\) は、2つのパラメータを用いて以下のように定義されます。

\[f(x | k, \lambda) = \dfrac{k}{\lambda} \left(\dfrac{x}{\lambda}\right)^{k-1} e^{-(x/\lambda)^k} \quad (x \ge 0)\]

この分布は、2つのパラメータによってその形状が決定されます。

  • \(k\) (または \(\beta\)): 形状パラメータ (shape parameter)

    • 分布の形状を決定するパラメータです。故障率が時間とともにどう変化するかを表現します。
      • \(k < 1\): 初期故障型 (Decreasing Failure Rate, DFR)。使用開始直後の故障率が最も高く、時間とともに減少していきます。初期不良が多い製品などが該当します。
      • \(k = 1\): 偶発故障型 (Constant Failure Rate, CFR)。故障率は時間によらず一定です。このとき、ワイブル分布は指数分布と一致します。
      • \(k > 1\): 摩耗故障型 (Increasing Failure Rate, IFR)。時間とともに故障率が増加していきます。摩耗や劣化によって故障しやすくなる機械製品などが該当します。
  • \(\lambda\) (または \(\eta\)): 尺度パラメータ (scale parameter)

    • 分布の横方向の伸縮を決定します。分布の特性寿命とも呼ばれ、累積故障確率が \(1 - 1/e \approx 63.2\%\) に達する時間を示します。
    • この値が大きくなると、分布は右に引き伸ばされ、寿命が全体的に長くなります。

主な特徴

  • 定義域: 確率変数は常に非負の値 (\(x \ge 0\)) をとります。
  • 形状の多様性: 形状パラメータ \(k\) の値を変えるだけで、指数分布や、正規分布に似た左右対称な形状、あるいは左に歪んだ形状など、非常に多様な形の分布を表現できます。特に \(k\) が3.6程度のとき、正規分布とよく似た形状になります。
  • 他の分布との関係:

    • 指数分布: 形状パラメータ \(k=1\) のとき、ワイブル分布はレートパラメータが \(1/\lambda\) の指数分布になります。
    • レイリー分布: 形状パラメータ \(k=2\) のとき、レイリー分布と呼ばれる分布になります。

2. ワイブル分布の応用例

  • 信頼性工学・製造業

    • ベアリング、コンデンサ、機械部品などの製品寿命を分析し、保証期間を設定したり、メンテナンス計画を立てたりするのに使われます。バスタブ曲線（初期故障期、偶発故障期、摩耗故障期の3段階からなる故障率曲線）の各期間をモデル化できます。
  • 生存時間分析

    • 医学や生物学の分野で、患者の生存時間や、ある処置後の再発までの時間を分析するために用いられます。
  • 気象学

    • 風速の分布をモデル化する標準的な手法として広く利用されています。これにより、風力発電のポテンシャル評価などが行われます。
  • 保険数理

    • 大規模な保険金請求が発生するまでの時間をモデル化するのに使われることがあります。

3. R言語によるシミュレーション

ここでは、尺度パラメータ \(\lambda\) を固定し、形状パラメータ \(k\) を変更した4つのワイブル分布を1枚のチャートに描画します。これにより、形状パラメータが分布の形状と故障率のタイプに与える影響を視覚的に理解します。

  • ケース1: k=0.7, λ=1 （初期故障型: DFR）
  • ケース2: k=1.0, λ=1 （偶発故障型: CFR, 指数分布と同一）
  • ケース3: k=2.0, λ=1 （摩耗故障型: IFR, レイリー分布）
  • ケース4: k=5.0, λ=1 （摩耗故障型: IFR, より対称に近い形状）

Rコード

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

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

# 2. 異なるパラメータを持つワイブル分布の確率密度を計算
# dweibull(x, shape, scale) を使用。shapeが形状k, scaleが尺度λ
df <- tibble(
  x = x_vals
) %>%
  mutate(
    `k=0.7 (初期故障型)` = dweibull(x, shape = 0.7, scale = 1),
    `k=1.0 (偶発故障型)` = dweibull(x, shape = 1.0, scale = 1),
    `k=2.0 (摩耗故障型)` = dweibull(x, shape = 2.0, scale = 1),
    `k=5.0 (摩耗故障型)` = dweibull(x, shape = 5.0, scale = 1)
  )

# 3. ggplotで描画しやすいように、データを「ロングフォーマット」に変換
df_long <- df %>%
  pivot_longer(
    cols = -x,
    names_to = "parameters",
    values_to = "density"
  ) %>%
  # 凡例の順序を調整
  mutate(parameters = factor(parameters, levels = c(
    "k=0.7 (初期故障型)",
    "k=1.0 (偶発故障型)",
    "k=2.0 (摩耗故障型)",
    "k=5.0 (摩耗故障型)"
  )))

# 4. 各分布に割り当てる色を定義
manual_colors <- c(
  `k=0.7 (初期故障型)` = "blue",
  `k=1.0 (偶発故障型)` = "red",
  `k=2.0 (摩耗故障型)` = "darkgreen",
  `k=5.0 (摩耗故障型)` = "purple"
)

# 5. ggplotを使用してチャートを描画
p <- ggplot(df_long, aes(x = x, y = density, color = parameters)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = manual_colors) +
  labs(
    title = "形状パラメータ(k)によるワイブル分布の変化 (λ=1で固定)",
    subtitle = "kの値によって、初期故障型(k<1)、偶発故障型(k=1)、摩耗故障型(k>1)を表現できる",
    x = "時間・寿命 (x)",
    y = "確率密度",
    color = "形状パラメータ (k)"
  ) +
  coord_cartesian(ylim = c(0, 2.0)) + # y軸の範囲を調整
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "gray40", size = 12)
  ) +
  guides(color = guide_legend(nrow = 2))

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

Figure 1 の解説

上記のRコードを実行すると、4つのワイブル分布が描画されたチャート Figure 1 が生成されます。

  • k=0.7 (青線): 初期故障型。\(x=0\) で確率密度が無限大に発散し、その後急速に減少します。これは、使い始めの非常に早い段階で故障する確率が最も高いことを示しています。
  • k=1.0 (赤線): 偶発故障型。\(x=0\) で最大値をとり、その後単調に減少します。これは、レートパラメータが1の指数分布と完全に同じ形状です。
  • k=2.0 (緑線): 摩耗故障型。分布は右に歪んだ山型をしています。故障はしばらく起こりにくいですが、ある時点を過ぎると故障する確率が上がっていく様子を表しています。
  • k=5.0 (紫線): 摩耗故障型。\(k\)がさらに大きくなったことで、分布のばらつきが小さくなり、より左右対称に近い（正規分布に似た）尖った山型になっています。これは、ほとんどの製品が特定の寿命（この場合 \(x=1\) の周辺）に集中して故障することを示唆しています。

以上です。