数理工房

100本ノック / シミュレーション / シミュレーション100本ノック

製造業の確率過程入門｜マルコフ連鎖から粒子フィルタまでPython実践10選

設備の「次」を確率で捉える：製造業の確率過程シミュレーション10本ノック（No.071〜No.080）

タイトル・概要

設備は正常・劣化・故障を行き来し、故障や受注は不規則に発生し、センサ値にはノイズが混ざります。本記事では、こうした時間とともに変わる不確実性を扱う確率過程を、架空の精密部品工場を題材に実装します。

No.071〜No.080を通じて、状態遷移、長期稼働率、異常兆候、故障件数、部品交換、外部価格、設備劣化、逐次状態推定を、保全・在庫・投資判断につなげます。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

はじめに：この記事で扱う製造業の実務課題

対象は、主要加工設備の突発停止を減らしながら、保全部品と点検工数を過剰にしない意思決定です。平均故障間隔だけでは、「いま劣化状態にいる設備が次にどうなるか」「来月の故障件数が上振れする確率はどれほどか」「ノイズを含むセンサから真の状態をどう推定するか」を答えられません。

本記事では確率過程を、予測値を一つ当てる技術ではなく、複数の将来シナリオとその確率を比較するための意思決定基盤として扱います。

現場でよくある状況

  • 設備台帳には故障履歴があるが、正常・劣化・故障の遷移構造は整理されていない
  • 月平均の故障件数は分かるが、保全部品が不足する確率は分からない
  • 振動や温度の閾値監視では、誤報と見逃しが多い
  • 市況価格や需要の変動を、設備投資評価の感度分析に入れられていない
  • 高度なモデルを導入しても、パラメータの根拠と更新手順が運用に落ちない

なぜこの問題は判断が難しいのか

平均値が同じでも、発生のばらつき、状態への依存、観測ノイズが異なれば必要な対策は変わります。確率過程では時点 tt の状態を確率変数 XtX_t とし、系列 (Xt)t≥0(X_t)_{t\geq0} 全体を扱います。

実務上は、(1)何を状態と定義するか、(2)時間を離散・連続のどちらで扱うか、(3)観測できない状態をどう推定するか、(4)推定誤差を意思決定にどう織り込むか、の設計が重要です。モデル精度だけでなく、停止損失、点検費、欠品費などの非対称な損失で評価する必要があります。

今回扱うノックの全体像

No.テーマ製造業での問い主な出力
071マルコフ連鎖設備状態は次期にどう遷移するか状態確率の推移
072定常分布長期的な故障・劣化比率は何％か長期状態比率
073隠れマルコフモデルセンサから見えない劣化を推定できるかフィルタ確率
074ポアソン過程故障・呼出件数は何件になるか件数分布
075更新過程交換を繰り返す部品の更新数はどうなるか更新数・間隔
076ブラウン運動累積誤差の揺らぎをどう表すか経路分布
077幾何ブラウン運動外部価格の比率変動をどう評価するか価格シナリオ
078確率微分方程式平均回帰する設備指標をどう表すか劣化指標経路
079カルマンフィルタノイズを除いて連続状態を推定できるか推定値・区間
080粒子フィルタ非線形な観測でも状態推定できるか粒子分布・推定値

Python 環境の準備

外部データには依存せず、NumPy・pandas・matplotlibだけを使います。乱数生成器は固定seedで初期化し、再実行しても同じ結果になるようにします。グラフ内のラベルはMarkdown変換時の環境差を避けるため英語表記とします。

%matplotlib inline
import sys
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from IPython.display import display

SEED = 20260712
rng = np.random.default_rng(SEED)
plt.rcParams.update({"figure.figsize": (8, 4), "axes.grid": True})
print(f"Python     : {sys.version.split()[0]}")
print(f"NumPy      : {np.__version__}")
print(f"pandas     : {pd.__version__}")
print(f"matplotlib : {matplotlib.__version__}")
print(f"random seed: {SEED}")
Python     : 3.13.1
NumPy      : 2.5.1
pandas     : 3.0.3
matplotlib : 3.11.0
random seed: 20260712

架空データの作成

精密部品工場の主力加工設備を想定し、90日分の「真の劣化度」と振動センサ値を生成します。真の劣化度は検証用にのみ保持し、実運用では観測できない潜在状態です。後半のフィルタでは、観測値だけからこの状態を推定します。

days = 90
t = np.arange(days)
true_health = np.zeros(days)
true_health[0] = 1.0
for i in range(1, days):
    true_health[i] = max(0, true_health[i-1] + 0.025 + rng.normal(0, 0.045))
vibration = 1.2 + 0.85 * true_health + 0.18 * true_health**2 + rng.normal(0, 0.22, days)
temperature = 42 + 3.2 * true_health + rng.normal(0, 0.8, days)
sensor_df = pd.DataFrame({"day": t, "true_degradation": true_health,
                          "vibration_mm_s": vibration, "temperature_C": temperature})
display(sensor_df.head().round(3))
print(f"rows={len(sensor_df)}, missing={sensor_df.isna().sum().sum()}")

fig, ax = plt.subplots()
ax.plot(t, vibration, label="Observed vibration", alpha=0.75)
ax.plot(t, 1.2 + 0.85*true_health + 0.18*true_health**2, label="Noise-free signal")
ax.set_title("Synthetic equipment sensor data")
ax.set_xlabel("Day"); ax.set_ylabel("Vibration [mm/s]")
ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
day true_degradation vibration_mm_s temperature_C
0 0 1.000 2.179 44.119
1 1 1.058 2.443 46.809
2 2 1.114 2.385 45.303
3 3 1.158 2.244 44.773
4 4 1.201 2.654 45.412 
rows=90, missing=0

png

No.071：マルコフ連鎖で設備状態の推移を予測する

実務での意味

点検時の設備状態を Normal（正常）、Degraded（劣化）、Failed（故障）の3段階に分け、翌週の状態確率を予測します。保全計画を「故障する／しない」の一点予測ではなく、停止リスクの確率として共有できます。

分析・モデル化の考え方

マルコフ連鎖は「次の状態が現在の状態だけに依存する」という近似です。遷移行列 PP の要素 pij=P(Xt+1=j∣Xt=i)p_{ij}=P(X_{t+1}=j\mid X_t=i) を定義すると、状態確率ベクトルは πt+1=πtP\boldsymbol{\pi}_{t+1}=\boldsymbol{\pi}_tP で更新されます。故障後は修理により正常または劣化へ戻る設定です。

Pythonで確認する

states = ["Normal", "Degraded", "Failed"]
P = np.array([[0.88, 0.11, 0.01], [0.18, 0.68, 0.14], [0.72, 0.28, 0.00]])
display(pd.DataFrame(P, index=states, columns=states).style.format("{:.0%}"))
pi = np.array([1.0, 0.0, 0.0])
history = [pi.copy()]
for _ in range(12):
    pi = pi @ P
    history.append(pi.copy())
markov_df = pd.DataFrame(history, columns=states).assign(week=np.arange(13))
display(markov_df.iloc[[0, 1, 4, 8, 12]].set_index("week").round(4))
fig, ax = plt.subplots()
for s in states: ax.plot(markov_df.week, markov_df[s], marker="o", label=s)
ax.set_title("Equipment-state probabilities by Markov chain")
ax.set_xlabel("Week"); ax.set_ylabel("Probability")
ax.set_ylim(0, 1); ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
  Normal Degraded Failed
Normal 88% 11% 1%
Degraded 18% 68% 14%
Failed 72% 28% 0%
Normal Degraded Failed
week
0 1.0000 0.0000 0.0000
1 0.8800 0.1100 0.0100
4 0.7253 0.2372 0.0374
8 0.6872 0.2687 0.0441
12 0.6820 0.2730 0.0450

png

結果の読み取り

正常開始でも、時間が進むと劣化・故障の確率が蓄積し、やがて一定の構成比に近づきます。翌週だけでなく4〜12週先の故障確率を示すことで、予備品手配や計画停止枠の確保を前倒しできます。ただし、設備年齢や負荷の影響が強い場合は状態を細分化するか、時間依存の遷移確率が必要です。

No.072：定常分布で長期の設備状態比率を見積もる

実務での意味

同じ運用を長く続けたとき、設備群の何％が正常・劣化・故障になるかを見積もります。長期の保全人員、代替能力、停止損失の基準値になります。

分析・モデル化の考え方

定常分布 π∗\boldsymbol{\pi}^* は π∗=π∗P\boldsymbol{\pi}^*=\boldsymbol{\pi}^*P、∑iπi∗=1\sum_i\pi_i^*=1 を満たします。初期状態にかかわらず同じ分布へ収束する条件の下では、短期シミュレーションの揺れに左右されない長期KPIを得られます。

Pythonで確認する

A = np.vstack([P.T - np.eye(len(states)), np.ones(len(states))])
b = np.r_[np.zeros(len(states)), 1.0]
stationary, *_ = np.linalg.lstsq(A, b, rcond=None)
stationary_df = pd.DataFrame({"state": states, "stationary_probability": stationary})
display(stationary_df.style.format({"stationary_probability": "{:.2%}"}))
weekly_loss = np.array([0, 350_000, 2_400_000])
expected_loss = stationary @ weekly_loss
print(f"Long-run expected weekly loss: JPY {expected_loss:,.0f}")
fig, ax = plt.subplots()
ax.bar(states, stationary, color=["#4c9f70", "#e0a458", "#d1495b"])
ax.set_title("Long-run equipment-state mix")
ax.set_xlabel("State"); ax.set_ylabel("Stationary probability")
ax.set_ylim(0, 1); ax.grid(True, axis="y", alpha=.3); plt.tight_layout(); plt.show()
  state stationary_probability
0 Normal 68.12%
1 Degraded 27.37%
2 Failed 4.51% 
Long-run expected weekly loss: JPY 204,076

png

結果の読み取り

定常状態の劣化・故障比率を、1週あたりの状態別損失と掛けると長期期待損失を算出できます。保全策変更後の遷移行列でも同じ計算を行い、期待損失の低下額と施策費を比較すると投資判断になります。定常分布は「今週の予測」ではなく、運用ルールを変えない場合の長期基準です。

No.073：隠れマルコフモデルで潜在的な劣化を推定する

実務での意味

設備の真の状態は直接見えず、振動判定（Low / Medium / High）だけが観測される状況を扱います。単発の異常値に反応せず、過去の観測系列を含めて劣化確率を更新できます。

分析・モデル化の考え方

隠れマルコフモデル（HMM）は潜在状態の遷移確率と、状態ごとの観測確率を組み合わせます。前向きフィルタでは、予測 α^t=αt−1P\hat{\alpha}_t=\alpha_{t-1}P と観測による補正 αt∝α^t⊙b(yt)\alpha_t\propto\hat{\alpha}_t\odot b(y_t) を反復します。

Pythonで確認する

emissions = np.array([[0.82, 0.16, 0.02], [0.18, 0.62, 0.20], [0.03, 0.22, 0.75]])
obs_names = ["Low", "Medium", "High"]
observations = np.array([0,0,1,0,1,1,2,1,2,2,1,2])
alpha = np.array([0.85, 0.14, 0.01])
filtered = []
for y in observations:
    pred = alpha @ P
    alpha = pred * emissions[:, y]
    alpha /= alpha.sum()
    filtered.append(alpha.copy())
hmm_df = pd.DataFrame(filtered, columns=states)
hmm_df["observation"] = [obs_names[i] for i in observations]
hmm_df.index = np.arange(1, len(hmm_df)+1)
display(hmm_df.tail().round(3))
fig, ax = plt.subplots()
ax.plot(hmm_df.index, hmm_df.Degraded + hmm_df.Failed, marker="o", label="P(Degraded or Failed)")
ax.axhline(.7, color="red", linestyle="--", label="Inspection threshold")
ax.set_title("Hidden-state risk inferred from vibration categories")
ax.set_xlabel("Inspection sequence"); ax.set_ylabel("Filtered probability")
ax.set_ylim(0,1); ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
Normal Degraded Failed observation
8 0.175 0.782 0.044 Medium
9 0.032 0.556 0.412 High
10 0.051 0.597 0.352 High
11 0.162 0.791 0.046 Medium
12 0.031 0.556 0.412 High

png

結果の読み取り

High観測が続くと劣化または故障の事後確率が上がり、途中のMediumだけでは急に正常へ戻りません。例えばリスク確率70％を精密点検の基準にできます。閾値は誤報率ではなく、見逃し時の停止損失と点検費の比較で決めるのが実務的です。

No.074：ポアソン過程で故障呼出件数を見積もる

実務での意味

独立した設備群から一定率で発生する故障呼出件数をモデル化し、保全要員や予備品の不足確率を求めます。

分析・モデル化の考え方

単位時間あたり発生率を λ\lambda とするポアソン過程では、期間 tt の件数 N(t)N(t) は P(N(t)=k)=e−λt(λt)k/k!P(N(t)=k)=e^{-\lambda t}(\lambda t)^k/k! に従い、到着間隔は平均 1/λ1/\lambda の指数分布です。発生率が季節・シフトで変わる場合は非一様ポアソン過程などへ拡張します。

Pythonで確認する

lambda_day = 0.42
sim_counts = rng.poisson(lambda_day * 30, size=20_000)
spares = 17
summary = pd.Series(sim_counts).describe(percentiles=[.5,.8,.9,.95,.99]).to_frame("30-day failures")
display(summary.round(2))
print(f"P(failures > {spares} spares) = {(sim_counts > spares).mean():.2%}")
fig, ax = plt.subplots()
bins = np.arange(sim_counts.min(), sim_counts.max()+2)-.5
ax.hist(sim_counts, bins=bins, density=True, color="#4c78a8", alpha=.8)
ax.axvline(spares, color="red", linestyle="--", label=f"Spare stock={spares}")
ax.set_title("Simulated monthly failure-call count")
ax.set_xlabel("Failures in 30 days"); ax.set_ylabel("Probability density")
ax.grid(True, axis="y", alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
30-day failures
count 20000.00
mean 12.60
std 3.55
min 1.00
50% 12.00
80% 16.00
90% 17.00
95% 19.00
99% 21.00
max 28.00 
P(failures > 17 spares) = 8.95%

png

結果の読み取り

平均件数だけを在庫数にすると、約半数の月で不足し得ます。分位点と欠品確率を併記すれば、「95％の月をカバーする在庫」などサービス水準で議論できます。故障が連鎖する、設備間で共通原因がある、発生率が上昇している場合はポアソン仮定が過小評価になるため注意が必要です。

No.075：更新過程で交換部品のライフサイクルを再現する

実務での意味

部品は故障・交換後に寿命がリセットされ、再び使用されます。この繰り返しをモデル化し、年間交換回数と予算を見積もります。

分析・モデル化の考え方

独立同分布の寿命 T1,T2,…T_1,T_2,\ldots に対し、更新時点は Sn=∑i=1nTiS_n=\sum_{i=1}^nT_i、時刻 tt までの更新数は N(t)=max⁡{n:Sn≤t}N(t)=\max\{n:S_n\le t\} です。ここでは寿命のばらつきを表しやすいワイブル分布を使います。

Pythonで確認する

shape, scale = 2.3, 95.0
def renewal_count(horizon, generator):
    elapsed = 0.0; count = 0; times = []
    while True:
        elapsed += scale * generator.weibull(shape)
        if elapsed > horizon: break
        count += 1; times.append(elapsed)
    return count, times
renewals = np.array([renewal_count(365, rng)[0] for _ in range(10_000)])
one_count, one_times = renewal_count(365, np.random.default_rng(SEED+1))
display(pd.Series(renewals).value_counts(normalize=True).sort_index().rename("probability").to_frame().style.format("{:.2%}"))
print(f"Expected annual replacements: {renewals.mean():.2f}")
fig, ax = plt.subplots()
ax.step([0]+one_times+[365], np.arange(one_count+2), where="post")
ax.scatter(one_times, np.arange(1, one_count+1), color="red", zorder=3)
ax.set_title("One sample path of component renewals")
ax.set_xlabel("Day"); ax.set_ylabel("Cumulative replacements")
ax.grid(True, alpha=.3); plt.tight_layout(); plt.show()
  probability
1 0.07%
2 5.21%
3 28.66%
4 38.92%
5 19.98%
6 6.07%
7 0.98%
8 0.11% 
Expected annual replacements: 3.96

png

結果の読み取り

交換回数は平均だけでなく分布で把握できるため、年間予算の上振れや保全部品の必要量を評価できます。新品同等に戻らない修理、予防交換、複数部品の依存故障がある場合は、単純な更新過程から仮想年齢モデルや競合リスクへ拡張します。

No.076：ブラウン運動で累積する微小な変動を表す

実務での意味

寸法補正誤差や熱変位など、多数の小さな影響が積み重なる量の不確実性を経路として評価します。

分析・モデル化の考え方

標準ブラウン運動 WtW_t は W0=0W_0=0、独立増分を持ち、Wt+Δt−Wt∼N(0,Δt)W_{t+\Delta t}-W_t\sim N(0,\Delta t) です。ドリフト付き過程 Xt=μt+σWtX_t=\mu t+\sigma W_t とすれば、系統的なずれとランダムな揺らぎを分離できます。ただし物理的境界や平均回帰は含まれません。

Pythonで確認する

n_paths, n_steps, horizon = 5000, 120, 8.0
dt = horizon / n_steps
mu, sigma = 0.018, 0.085
increments = mu*dt + sigma*np.sqrt(dt)*rng.normal(size=(n_paths, n_steps))
paths = np.c_[np.zeros(n_paths), np.cumsum(increments, axis=1)]
time = np.linspace(0, horizon, n_steps+1)
limit = 0.30
print(f"P(|error at horizon| > {limit:.2f}) = {(np.abs(paths[:,-1]) > limit).mean():.2%}")
fig, ax = plt.subplots()
for p in paths[:30]: ax.plot(time, p, color="#4c78a8", alpha=.18)
q05, q50, q95 = np.quantile(paths, [.05,.5,.95], axis=0)
ax.plot(time, q50, color="black", label="Median")
ax.fill_between(time, q05, q95, alpha=.25, label="90% range")
ax.axhline(limit, color="red", linestyle="--"); ax.axhline(-limit, color="red", linestyle="--")
ax.set_title("Accumulated process error under Brownian motion")
ax.set_xlabel("Operating time [h]"); ax.set_ylabel("Accumulated error [mm]")
ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
P(|error at horizon| > 0.30) = 29.82%

png

結果の読み取り

終業時点の規格外確率だけでなく、途中の経路と90％範囲を確認できます。校正間隔は、許容リスクを超える時刻を基準に設計できます。一方、誤差が中心へ戻る制御系ならブラウン運動は長期分散を過大評価するため、No.078の平均回帰モデルが適します。

No.077：幾何ブラウン運動で外部価格のシナリオを作る

実務での意味

電力・原材料など、絶対額より変化率で動き、負値を取らない外部価格をシナリオ化します。調達契約や設備投資の採算感度に利用できます。

分析・モデル化の考え方

幾何ブラウン運動は dSt=μStdt+σStdWtdS_t=\mu S_tdt+\sigma S_tdW_t、解は St=S0exp⁡{(μ−σ2/2)t+σWt}S_t=S_0\exp\{(\mu-\sigma^2/2)t+\sigma W_t\} です。価格を保証する予測モデルではなく、仮定した変動率の下で採算の分布を作るモデルです。

Pythonで確認する

n_paths, months = 10000, 24
dt = 1/12
s0, annual_mu, annual_sigma = 100.0, 0.035, 0.22
z = rng.normal(size=(n_paths, months))
log_inc = (annual_mu-.5*annual_sigma**2)*dt + annual_sigma*np.sqrt(dt)*z
price_paths = s0*np.exp(np.c_[np.zeros(n_paths), np.cumsum(log_inc, axis=1)])
terminal = price_paths[:,-1]
display(pd.Series(terminal).quantile([.05,.25,.5,.75,.95]).rename("price_index_at_24m").to_frame().round(1))
print(f"P(price index > 130) = {(terminal > 130).mean():.2%}")
fig, ax = plt.subplots()
for p in price_paths[:35]: ax.plot(np.arange(months+1), p, alpha=.15)
ax.plot(np.arange(months+1), np.median(price_paths, axis=0), color="black", label="Median")
ax.set_title("Input-price scenarios by geometric Brownian motion")
ax.set_xlabel("Month"); ax.set_ylabel("Price index")
ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
price_index_at_24m
0.05 61.4
0.25 83.4
0.50 102.4
0.75 125.4
0.95 171.2 
P(price index > 130) = 21.32%

png

結果の読み取り

24か月後の中央値と上側分位点には差があり、単一の価格前提だけでは上振れリスクを見落とします。投資案件ごとに各経路の利益・NPVを計算すれば、赤字確率で比較できます。実価格には平均回帰、ジャンプ、季節性があるため、過去変動率の機械的な外挿は避けます。

No.078：確率微分方程式で平均回帰する劣化指標を表す

実務での意味

温度や制御偏差のように、外乱で揺れても制御により基準へ戻る指標を扱います。単なるランダムウォークより現実的な逸脱時間を評価できます。

分析・モデル化の考え方

Ornstein–Uhlenbeck過程 dXt=θ(μ−Xt)dt+σdWtdX_t=\theta(\mu-X_t)dt+\sigma dW_t をEuler–Maruyama法で離散化します。θ\theta は復元速度、μ\mu は長期平均、σ\sigma は外乱強度です。Δt\Delta t が粗すぎると数値誤差が増えるため刻み幅の検証が必要です。

Pythonで確認する

n_paths, n_steps, dt = 3000, 240, 1/24
theta, long_mean, sigma_ou, x0 = 1.4, 2.0, 0.55, 3.2
ou = np.zeros((n_paths, n_steps+1)); ou[:,0] = x0
for k in range(n_steps):
    ou[:,k+1] = ou[:,k] + theta*(long_mean-ou[:,k])*dt + sigma_ou*np.sqrt(dt)*rng.normal(size=n_paths)
threshold = 2.8
exceed_hours = (ou[:,1:] > threshold).sum(axis=1)
print(f"Mean hours above threshold: {exceed_hours.mean():.2f}")
print(f"P(above threshold for >= 8 hours): {(exceed_hours >= 8).mean():.2%}")
fig, ax = plt.subplots()
time_ou = np.arange(n_steps+1)*dt
for p in ou[:25]: ax.plot(time_ou, p, alpha=.15)
ax.plot(time_ou, np.median(ou, axis=0), color="black", label="Median")
ax.axhline(long_mean, color="green", linestyle="--", label="Long-run mean")
ax.axhline(threshold, color="red", linestyle="--", label="Alarm threshold")
ax.set_title("Mean-reverting equipment indicator")
ax.set_xlabel("Day"); ax.set_ylabel("Condition indicator")
ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
Mean hours above threshold: 10.48
P(above threshold for >= 8 hours): 58.40%

png

結果の読み取り

初期値が高くても平均へ戻る一方、外乱により閾値超過が繰り返されます。単発超過より「8時間以上の累積超過」のような運用KPIが安定します。パラメータは設備・運転条件ごとに推定し、残差の独立性と分布を確認する必要があります。

No.079：カルマンフィルタでセンサノイズを除去する

実務での意味

真の劣化度を直接測れないとき、過去の推定と新しいセンサ値を不確実性に応じて統合します。生値の移動平均より、推定誤差の幅を明示できる点が利点です。

分析・モデル化の考え方

線形状態空間モデル xt=Axt−1+wtx_t=Ax_{t-1}+w_t, yt=Hxt+vty_t=Hx_t+v_t を仮定します。予測誤差分散と観測誤差分散からカルマンゲイン KtK_t を計算し、信頼できる側へ重みを置きます。ここでは振動から線形変換した劣化観測を使います。

Pythonで確認する

y_linear = (vibration - 1.2) / 0.85
A, H, Q, R = 1.0, 1.0, 0.045**2 + 0.01**2, (0.22/0.85)**2
x_est, cov = 0.8, 0.4
estimates, variances = [], []
for y in y_linear:
    x_pred = A*x_est + 0.025
    p_pred = A*cov*A + Q
    gain = p_pred*H/(H*p_pred*H + R)
    x_est = x_pred + gain*(y-H*x_pred)
    cov = (1-gain*H)*p_pred
    estimates.append(x_est); variances.append(cov)
estimates, variances = np.array(estimates), np.array(variances)
rmse_raw = np.sqrt(np.mean((y_linear-true_health)**2))
rmse_kf = np.sqrt(np.mean((estimates-true_health)**2))
print(f"Raw-observation RMSE : {rmse_raw:.3f}")
print(f"Kalman-filter RMSE   : {rmse_kf:.3f}")
fig, ax = plt.subplots()
ax.plot(t, y_linear, color="gray", alpha=.35, label="Converted observation")
ax.plot(t, true_health, color="black", label="True degradation (validation only)")
ax.plot(t, estimates, color="#4c78a8", label="Kalman estimate")
sd = np.sqrt(variances)
ax.fill_between(t, estimates-1.96*sd, estimates+1.96*sd, alpha=.2, label="Approx. 95% interval")
ax.set_title("Linear state estimation by Kalman filter")
ax.set_xlabel("Day"); ax.set_ylabel("Degradation index")
ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
Raw-observation RMSE : 1.372
Kalman-filter RMSE   : 1.249

png

結果の読み取り

フィルタ推定は観測ノイズを抑え、真値に対するRMSEも比較できます。95％区間の上限が保全基準を超えたら点検するなど、不確実性を含むルールにできます。線形性・正規性が大きく崩れるとバイアスが生じるため、非線形観測には次の粒子フィルタを検討します。

No.080：粒子フィルタで非線形な設備状態を推定する

実務での意味

振動値が劣化度に対して二次的に増えるなど、観測関係が非線形な場合に潜在状態を逐次推定します。複数の可能性を粒子群として保持できるため、分布が歪む状況にも対応できます。

分析・モデル化の考え方

粒子フィルタは、状態分布を多数の粒子 xt(i)x_t^{(i)} と重み wt(i)w_t^{(i)} で近似します。予測、尤度による重み更新、正規化、リサンプリングを反復します。粒子数を増やすほど近似は安定しますが計算量も増えます。

Pythonで確認する

n_particles = 4000
pf_rng = np.random.default_rng(SEED+80)
particles = np.clip(pf_rng.normal(.8, .35, n_particles), 0, None)
pf_est, pf_low, pf_high, ess_history = [], [], [], []
for y in vibration:
    particles = np.clip(particles + 0.025 + pf_rng.normal(0, .045, n_particles), 0, None)
    predicted_y = 1.2 + .85*particles + .18*particles**2
    weights = np.exp(-.5*((y-predicted_y)/.22)**2)
    weights += 1e-300; weights /= weights.sum()
    pf_est.append(np.sum(weights*particles))
    order = np.argsort(particles); sorted_p = particles[order]; cum_w = np.cumsum(weights[order])
    pf_low.append(np.interp(.025, cum_w, sorted_p)); pf_high.append(np.interp(.975, cum_w, sorted_p))
    ess_history.append(1/np.sum(weights**2))
    positions = (pf_rng.random() + np.arange(n_particles))/n_particles
    indexes = np.searchsorted(np.cumsum(weights), positions)
    particles = particles[indexes]
pf_est, pf_low, pf_high = map(np.array, (pf_est, pf_low, pf_high))
print(f"Particle-filter RMSE: {np.sqrt(np.mean((pf_est-true_health)**2)):.3f}")
print(f"Mean effective sample size before resampling: {np.mean(ess_history):.0f}/{n_particles}")
fig, ax = plt.subplots()
ax.plot(t, true_health, color="black", label="True degradation (validation only)")
ax.plot(t, pf_est, color="#f58518", label="Particle-filter estimate")
ax.fill_between(t, pf_low, pf_high, color="#f58518", alpha=.2, label="95% particle interval")
ax.set_title("Nonlinear state estimation by particle filter")
ax.set_xlabel("Day"); ax.set_ylabel("Degradation index")
ax.grid(True, alpha=.3); ax.legend(); plt.tight_layout(); plt.show()
Particle-filter RMSE: 0.067
Mean effective sample size before resampling: 3199/4000

png

結果の読み取り

非線形な振動モデルを直接使いながら、真の劣化度に近い推定と区間が得られます。有効サンプルサイズは粒子の偏りを監視する診断指標です。実装時は、粒子数、尤度、リサンプリング方法を検証し、既知の故障履歴で「何日前に警報できたか」を評価します。

対象ノックを通して見える実務上の示唆

  1. 状態・件数・連続量を分ける：設備状態ならマルコフ連鎖、独立な発生件数ならポアソン過程、連続的な揺らぎなら確率微分方程式が出発点です。
  2. 観測値と真の状態を分ける：センサ値を状態そのものとみなさず、HMM、カルマンフィルタ、粒子フィルタで観測誤差を明示します。
  3. 平均ではなく分布で決める：期待故障数に加え、在庫不足確率、上側分位点、区間推定をKPIにします。
  4. 複雑さは段階的に上げる：説明可能な単純モデルを基準にし、残差や意思決定価値が改善するときだけ高度化します。
  5. モデルを費用へ接続する：確率を停止損失、点検費、在庫費、品質損失へ変換して施策を比較します。

実務導入する場合に必要なこと

  • 状態定義、故障コード、時刻粒度、交換・修理の記録基準を統一する
  • 欠測、センサ交換、稼働停止時間、運転条件の変更を履歴として残す
  • 学習期間と将来検証期間を分け、設備単位・期間単位のデータ漏洩を防ぐ
  • 精度だけでなく、警報リードタイム、見逃し損失、誤報対応工数で評価する
  • 現場知識でパラメータ範囲と物理的制約を確認する
  • ドリフト監視、定期再推定、モデル承認、アラート後の標準作業まで設計する

小さく始めるなら、重要設備1系列、1つの意思決定（例：7日以内の精密点検）に絞り、既存ルールとの並行評価から進めます。

まとめ

No.071〜No.080では、設備状態の遷移、長期比率、潜在状態、イベント件数、交換サイクル、連続的な揺らぎ、逐次推定を一通り実装しました。確率過程の価値は、未来を断定することではなく、起こり得る幅と確率を示し、保全・在庫・投資の判断基準を定量化することにあります。

最初の実務ステップは、モデル選定よりも「どの判断を、何日前に、どの損失基準で変えたいか」を明確にすることです。その問いに合わせて必要最小限の状態・観測・時間粒度を設計すると、分析が運用につながりやすくなります。

法人向けのご相談

数理工房では、製造業における設備保全、故障・需要シミュレーション、センサデータによる状態推定、意思決定KPIの設計まで、データ整備段階を含めてご支援します。PoCだけで終わらせず、現場運用・評価・更新まで見据えた設計をご相談いただけます。

📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。

Related Posts

関連する記事