100本ノック / 確率統計 / 確率・統計Python100本ノック

観測データから工程の真の姿を推定する：製造業の統計的推定10本ノック

2020-05-13

観測データから工程の真の姿を推定する：製造業の統計的推定10本ノック

製造現場で観測できる寸法、不適合数、設備停止件数は、母集団の一部にすぎません。本記事では、架空の精密部品工場を題材に、最尤推定、尤度の可視化、数値最適化、勾配、ニュートン法、フィッシャー情報量、ブートストラップ、EMアルゴリズムを実装します。

目標は数式を計算することではなく、「工程中心はどこか」「不適合率は何%か」「追加測定は必要か」「混在した工程を分けて考えるべきか」という意思決定へつなげることです。対象は確率・統計Python実装100本ノックの No.041〜No.050 です。

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

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

精密シャフトを生産する架空工場で、品質保証部が翌月の工程条件と検査計画を決める場面を考えます。手元にあるのは、抜取測定した外径、ロットごとの不適合数、日ごとの突発停止件数、設備IDが欠けた締付トルクです。

知りたい母平均や不適合率は直接見えません。観測値から未知のパラメータを推定し、その不確かさを示し、推定モデルが現場の発生機構と整合するかを確認する必要があります。本記事では、点推定だけで終わらず、推定精度とモデル利用条件まで扱います。

現場でよくある状況

平均寸法は規格中心に近いが、偶然のずれか工程偏りか判断しづらい
不適合率を「不適合数÷検査数」で報告しているが、標本数の違いが無視される
停止件数の平均を置いているものの、発生回数モデルとして妥当か確認していない
最適化ライブラリの答えを採用しているが、目的関数や収束を説明できない
設備IDの欠損により、異なる工程状態が一つの分布として集計される

統計的推定は、こうした状況で「観測データを最もよく説明する工程パラメータ」と「その推定をどこまで信頼できるか」を整理する共通言語になります。

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

同じ推定値でも、標本数が20と2,000では精度が異なります。また、正規分布・二項分布・ポアソン分布など、データの生成過程に合わない確率モデルを置けば、精密な計算でも判断を誤ります。

観測値を $x_1,\ldots,x_n$ 、未知パラメータを $\theta$ とすると、尤度は

L(\theta)=\prod_{i=1}^{n} f(x_i\mid\theta)

です。最尤推定値は $\hat{\theta}=\arg\max_\theta L(\theta)$ ですが、実装では積の桁落ちを避けるため対数尤度 $\ell(\theta)=\log L(\theta)$ を使います。さらに、独立性、分布形、データ欠損、測定系の妥当性は別途検証が必要です。

今回扱うノックの全体像

No.	テーマ	製造業での問い
041	最尤推定（正規分布）	外径の工程中心とばらつきはどこか
042	最尤推定（二項分布）	検査数が違うロットを含む全体不適合率は何%か
043	最尤推定（ポアソン分布）	1日あたりの突発停止発生率は何件か
044	scipy.optimizeによる最適化	閉形式がない場合も同じ推定問題を解けるか
045	尤度関数の可視化	どのパラメータ範囲がデータと整合するか
046	勾配計算	パラメータをどちらへ動かせば適合が改善するか
047	ニュートン法	曲率を使って停止発生率へ高速に収束できるか
048	フィッシャー情報量	測定数を増やすと推定精度はどう変わるか
049	ブートストラップ	理論式に頼りすぎず不確かさを評価できるか
050	EMアルゴリズム	設備ラベルなしで混在工程を推定できるか

前半で基本モデルを推定し、中盤で計算原理と精度を確認し、最後に潜在変数を持つモデルへ拡張します。

Python 環境の準備

NumPyで計算、pandasで表、SciPyで確率分布と最適化、matplotlibで可視化します。japanize_matplotlib は日本語ラベルの表示にだけ使い、seabornや外部データは使いません。乱数生成器は一度だけ作り、seedを固定します。

import platform
import numpy as np
import pandas as pd
import scipy
from scipy import stats
from scipy.optimize import minimize
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib

SEED = 20260711
rng = np.random.default_rng(SEED)
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.unicode_minus"] = False

print(f"Python     : {platform.python_version()}")
print(f"NumPy      : {np.__version__}")
print(f"pandas     : {pd.__version__}")
print(f"SciPy      : {scipy.__version__}")
print(f"Matplotlib : {matplotlib.__version__}")
print(f"random seed: {SEED}")

Python     : 3.13.1
NumPy      : 2.5.1
pandas     : 3.0.3
SciPy      : 1.18.0
Matplotlib : 3.11.0
random seed: 20260711

架空データの作成

4種類のデータをPython内で生成します。外径は安定工程の正規分布、ロット不適合数は二項分布、日次停止件数はポアソン分布、締付トルクは2設備の正規混合分布を仮定します。実務ではこの仮定を、管理図、時系列、層別、測定システム解析などで点検してから用います。

# 1) 外径測定、2) ロット別不適合、3) 日次停止件数
diameter = rng.normal(loc=20.018, scale=0.042, size=180)
lot_sizes = rng.integers(80, 151, size=36)
defect_counts = rng.binomial(lot_sizes, p=0.028)
stoppages = rng.poisson(lam=1.7, size=60)

# 4) 設備IDが欠けた締付トルク（検証用の真のラベルは分析には使わない）
n_torque = 320
true_machine = rng.choice([0, 1], size=n_torque, p=[0.62, 0.38])
torque = rng.normal(
    loc=np.where(true_machine == 0, 48.8, 52.6),
    scale=np.where(true_machine == 0, 0.75, 1.05),
)

data_summary = pd.DataFrame({
    "データ": ["外径", "ロット検査", "日次停止", "締付トルク"],
    "観測単位": ["個", "ロット", "日", "個"],
    "標本数": [len(diameter), len(lot_sizes), len(stoppages), len(torque)],
    "観測値の要約": [
        f"平均 {diameter.mean():.3f} mm",
        f"不適合 {defect_counts.sum()} / {lot_sizes.sum()}",
        f"平均 {stoppages.mean():.2f} 件/日",
        f"平均 {torque.mean():.2f} N·m",
    ],
})
data_summary

	データ	観測単位	標本数	観測値の要約
0	外径	個	180	平均 20.018 mm
1	ロット検査	ロット	36	不適合 114 / 4079
2	日次停止	日	60	平均 1.87 件/日
3	締付トルク	個	320	平均 50.03 N·m

No.041：最尤推定（正規分布）— 外径の工程中心とばらつきを推定する

実務での意味

連続量の品質特性が概ね正規分布に従うとき、平均 $\mu$ は工程中心、標準偏差 $\sigma$ は短期ばらつきの代表値になります。規格内に入った割合だけでなく、工程がどこに位置し、どれだけ散らばるかを把握する基礎です。

分析・モデル化の考え方

独立な観測 $x_i\sim N(\mu,\sigma^2)$ の対数尤度を最大化すると、

\hat\mu=\frac{1}{n}\sum_i x_i,\qquad \hat\sigma^2=\frac{1}{n}\sum_i(x_i-\hat\mu)^2

となります。最尤分散の分母は $n$ であり、不偏分散の $n-1$ とは目的が違います。

Pythonで確認する

mu_hat = diameter.mean()
sigma_hat = diameter.std(ddof=0)
normal_mle = pd.DataFrame({
    "推定対象": ["平均 mu", "標準偏差 sigma"],
    "最尤推定値": [mu_hat, sigma_hat],
    "単位": ["mm", "mm"],
})
display(normal_mle.round(4))

x_grid = np.linspace(diameter.min() - 0.03, diameter.max() + 0.03, 300)
plt.hist(diameter, bins=16, density=True, alpha=0.55, edgecolor="white", label="観測外径")
plt.plot(x_grid, stats.norm.pdf(x_grid, mu_hat, sigma_hat), lw=2.2, label="推定正規分布")
plt.axvline(mu_hat, color="tab:red", ls="--", label=f"推定平均={mu_hat:.3f}")
plt.title("外径データと最尤推定した正規分布")
plt.xlabel("外径 [mm]")
plt.ylabel("確率密度")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	推定対象	最尤推定値	単位
0	平均 mu	20.0179	mm
1	標準偏差 sigma	0.0430	mm

png

結果の読み取り

推定平均は規格値そのものではなく、今回の標本が最も支持する工程中心です。ヒストグラムと推定曲線が大きく乖離していないことを目視できますが、正規性を保証するものではありません。工程調整は平均だけで決めず、規格限界、測定誤差、時系列ドリフト、標準偏差の変化を併せて確認します。

No.042：最尤推定（二項分布）— ロット横断の不適合率を推定する

実務での意味

ロット $j$ の検査数 $n_j$ と不適合数 $k_j$ から、共通の不適合確率 $p$ を推定します。ロット率の単純平均では検査数の小さいロットを過大に重み付けするため、全検査個数を基準に集約します。

分析・モデル化の考え方

$K=\sum_j k_j$ 、 $N=\sum_j n_j$ とすれば、二項尤度は $L(p)\propto p^K(1-p)^{N-K}$ 、最尤推定値は

$\hat p=\frac{K}{N}$

です。これは各個体の不適合確率が共通で独立という仮定に基づきます。工程・製品・検査員で率が違う場合は層別や回帰モデルが必要です。

Pythonで確認する

total_defects = defect_counts.sum()
total_inspected = lot_sizes.sum()
p_hat = total_defects / total_inspected
lot_rates = defect_counts / lot_sizes

binomial_result = pd.DataFrame({
    "指標": ["総検査数", "総不適合数", "二項MLE", "ロット率の単純平均"],
    "値": [total_inspected, total_defects, p_hat, lot_rates.mean()],
})
display(binomial_result.round(5))

plt.scatter(lot_sizes, lot_rates * 100, alpha=0.75, label="ロット別")
plt.axhline(p_hat * 100, color="tab:red", lw=2, label=f"全体MLE={p_hat:.2%}")
plt.title("ロット別検査数と不適合率")
plt.xlabel("ロット検査数 [個]")
plt.ylabel("不適合率 [%]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	指標	値
0	総検査数	4079.00000
1	総不適合数	114.00000
2	二項MLE	0.02795
3	ロット率の単純平均	0.02672

png

結果の読み取り

二項MLEは「総不適合数÷総検査数」で、検査数に応じた加重平均です。散布図では、検査数が少ないロットほど率が離散的に揺れやすいことも分かります。この全体率だけでロットを順位付けせず、ロットごとの工程条件と信頼区間を確認するのが安全です。

No.043：最尤推定（ポアソン分布）— 日次の突発停止発生率を推定する

実務での意味

一定期間に発生する比較的まれな事象には、停止件数、故障件数、欠点数などがあります。1日あたりの平均発生率 $\lambda$ を推定すれば、保全要員や予備品の必要量を考える基準になります。

分析・モデル化の考え方

$Y_i\sim\mathrm{Poisson}(\lambda)$ なら、

$P(Y_i=y_i)=\frac{e^{-\lambda}\lambda^{y_i}}{y_i!},\qquad \hat\lambda=\bar y$

です。ポアソン分布は平均と分散が等しいという強い仮定を持ちます。過分散、ゼロ過多、曜日効果、連続故障がある場合は別モデルを検討します。

Pythonで確認する

lambda_hat = stoppages.mean()
poisson_check = pd.DataFrame({
    "指標": ["平均（lambda MLE）", "標本分散", "分散 / 平均", "ゼロ停止日の割合"],
    "値": [lambda_hat, stoppages.var(ddof=1), stoppages.var(ddof=1) / lambda_hat, np.mean(stoppages == 0)],
})
display(poisson_check.round(3))

k = np.arange(0, stoppages.max() + 1)
observed = np.bincount(stoppages, minlength=len(k))[:len(k)] / len(stoppages)
plt.bar(k - 0.18, observed, width=0.36, alpha=0.7, label="観測割合")
plt.bar(k + 0.18, stats.poisson.pmf(k, lambda_hat), width=0.36, alpha=0.7, label="推定Poisson")
plt.title("日次停止件数の観測分布と推定分布")
plt.xlabel("1日あたり停止件数 [件]")
plt.ylabel("割合・確率")
plt.grid(True, axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	指標	値
0	平均（lambda MLE）	1.867
1	標本分散	2.219
2	分散 / 平均	1.189
3	ゼロ停止日の割合	0.183

png

結果の読み取り

平均停止件数が最尤推定値になります。標本分散÷平均が1から大きく外れる場合は、単一の一定発生率という仮定を疑います。ここで得た発生率は要員計画の入力候補ですが、「1件あたり復旧時間」は含まれないため、工数計画には停止時間モデルも必要です。

No.044：scipy.optimizeによる最適化 — 制約付きの数値推定を実装する

実務での意味

標準的な分布では推定値を式で書けますが、打切り、複数要因、複雑な制約が入ると閉形式が得られないことがあります。数値最適化を使えば、負の対数尤度を目的関数として同じ枠組みで解けます。

分析・モデル化の考え方

最大化は $-\ell(\theta)$ の最小化に置き換えます。 $\sigma>0$ を必ず満たすよう $\eta=\log\sigma$ を最適化し、最後に $\sigma=e^\eta$ へ戻します。初期値、スケーリング、収束フラグ、境界解の確認が実務上重要です。

Pythonで確認する

def normal_nll(params, x):
    mu, log_sigma = params
    sigma = np.exp(log_sigma)
    return -np.sum(stats.norm.logpdf(x, loc=mu, scale=sigma))

opt_result = minimize(
    normal_nll,
    x0=np.array([20.0, np.log(0.05)]),
    args=(diameter,),
    method="Nelder-Mead",
    options={"xatol": 1e-11, "fatol": 1e-11, "maxiter": 2000},
)
mu_opt, sigma_opt = opt_result.x[0], np.exp(opt_result.x[1])
optimization_check = pd.DataFrame({
    "項目": ["収束", "反復回数", "最適化 mu", "解析解 mu", "最適化 sigma", "解析解 sigma"],
    "値": [opt_result.success, opt_result.nit, mu_opt, mu_hat, sigma_opt, sigma_hat],
})
optimization_check

	項目	値
0	収束	True
1	反復回数	77
2	最適化 mu	20.017902
3	解析解 mu	20.017902
4	最適化 sigma	0.043047
5	解析解 sigma	0.043047

結果の読み取り

数値解がNo.041の解析解とほぼ一致すれば、目的関数と変数変換の実装を検算できます。successだけでなく、初期値を変えて同じ解へ到達するか、勾配が十分小さいか、業務上あり得る範囲かも確認します。最適化器はモデルの妥当性までは保証しません。

No.045：尤度関数の可視化 — 推定値の周辺にある不確かさを見る

実務での意味

最尤推定値は一点ですが、近い尤度を持つ候補もデータと整合します。尤度の幅を見ることで、工程平均を小数点以下何桁まで議論できるか、追加測定が必要かを直感的に把握できます。

分析・モデル化の考え方

$\sigma$ を推定値に固定し、候補 $\mu$ ごとの相対対数尤度

$\ell(\mu)-\ell(\hat\mu)$

を描きます。1パラメータの尤度比近似では、95%の目安は $-\chi^2_{1,0.95}/2\approx-1.92$ です。これは有限標本で厳密な区間とは限りません。

Pythonで確認する

mu_candidates = np.linspace(mu_hat - 0.015, mu_hat + 0.015, 301)
loglik = np.array([
    np.sum(stats.norm.logpdf(diameter, loc=mu, scale=sigma_hat))
    for mu in mu_candidates
])
relative_loglik = loglik - loglik.max()
cutoff = -stats.chi2.ppf(0.95, df=1) / 2
supported = mu_candidates[relative_loglik >= cutoff]

print(f"尤度比による95%支持範囲（近似）: {supported.min():.5f} ～ {supported.max():.5f} mm")
plt.plot(mu_candidates, relative_loglik, lw=2)
plt.axvline(mu_hat, color="tab:red", ls="--", label="MLE")
plt.axhline(cutoff, color="tab:gray", ls=":", label="95%目安")
plt.fill_between(mu_candidates, relative_loglik, cutoff,
                 where=relative_loglik >= cutoff, alpha=0.2)
plt.title("工程平均の相対対数尤度")
plt.xlabel("平均候補 mu [mm]")
plt.ylabel("相対対数尤度")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

尤度比による95%支持範囲（近似）: 20.01170 ～ 20.02410 mm

png

結果の読み取り

曲線の頂点が最尤推定値で、頂点付近が平坦なら候補を絞りにくく、鋭ければ精度が高いと読めます。支持範囲が管理上許容できる幅より広い場合、過度な微調整を避けて測定数を増やす判断につながります。規格適合の判定区間とは目的が違う点に注意します。

No.046：勾配計算 — 尤度が改善する方向を確認する

実務での意味

勾配は「パラメータを少し動かしたとき目的関数がどれだけ変わるか」です。最適化が止まった理由の診断、独自モデルの実装検証、オンライン推定の更新方向を理解するうえで役立ちます。

分析・モデル化の考え方

$\sigma$ を固定した正規分布の負の対数尤度 $J(\mu)$ の勾配は

\frac{dJ}{d\mu}=\frac{n\mu-\sum_i x_i}{\sigma^2}

です。中心差分 $(J(\mu+h)-J(\mu-h))/(2h)$ と比較し、式とコードを検算します。最尤点では勾配がほぼ0になります。

Pythonで確認する

def nll_mu(mu):
    return -np.sum(stats.norm.logpdf(diameter, loc=mu, scale=sigma_hat))

def analytic_gradient(mu):
    return (len(diameter) * mu - diameter.sum()) / sigma_hat**2

h = 1e-6
check_points = np.array([mu_hat - 0.006, mu_hat, mu_hat + 0.006])
gradient_table = pd.DataFrame({
    "mu候補": check_points,
    "解析勾配": [analytic_gradient(mu) for mu in check_points],
    "数値勾配": [(nll_mu(mu + h) - nll_mu(mu - h)) / (2 * h) for mu in check_points],
})
gradient_table["差の絶対値"] = abs(gradient_table["解析勾配"] - gradient_table["数値勾配"])
gradient_table.round(6)

	mu候補	解析勾配	数値勾配	差の絶対値
0	20.011902	-582.825624	-582.825624	0.000001
1	20.017902	0.000000	-0.000000	0.000000
2	20.023902	582.825624	582.825625	0.000001

結果の読み取り

解析勾配と数値勾配が近ければ、微分式の実装に大きな誤りがないと判断できます。MLEより左では負、右では正となるため、負の対数尤度を下げる方向が中心へ向いています。差分幅が小さすぎると丸め誤差、大きすぎると近似誤差が増えるため、複数の幅で確認します。

No.047：ニュートン法 — 曲率を使って停止発生率を反復推定する

実務での意味

ニュートン法は勾配だけでなく曲率も利用し、解へ反復的に近づきます。日次停止率の例で更新履歴を残すことで、「答え」だけでなく収束過程を監査できる実装にします。

分析・モデル化の考え方

ポアソン対数尤度の1階・2階微分は

\ell'(\lambda)=\frac{\sum_i y_i}{\lambda}-n,\qquad \ell''(\lambda)=-\frac{\sum_i y_i}{\lambda^2}

です。更新式 $\lambda_{t+1}=\lambda_t-\ell'(\lambda_t)/\ell''(\lambda_t)$ を用います。正値制約を破る可能性がある複雑な問題では、ステップ制御や対数変換が必要です。

Pythonで確認する

lam = 0.6
history = []
sum_y, n_days = stoppages.sum(), len(stoppages)
for iteration in range(1, 9):
    score = sum_y / lam - n_days
    curvature = -sum_y / lam**2
    new_lam = lam - score / curvature
    history.append((iteration, lam, score, new_lam))
    if abs(new_lam - lam) < 1e-10:
        lam = new_lam
        break
    lam = new_lam

newton_history = pd.DataFrame(history, columns=["反復", "更新前lambda", "スコア", "更新後lambda"])
display(newton_history.round(6))
print(f"Newton法={lam:.6f}, 解析解（標本平均）={lambda_hat:.6f}")

	反復	更新前lambda	スコア	更新後lambda
0	1	0.600000	126.666667	1.007143
1	2	1.007143	51.205674	1.470891
2	3	1.470891	16.144322	1.782753
3	4	1.782753	2.824173	1.862894
4	5	1.862894	0.121495	1.866659
5	6	1.866659	0.000245	1.866667
6	7	1.866667	0.000000	1.866667

Newton法=1.866667, 解析解（標本平均）=1.866667

結果の読み取り

更新値が標本平均へ収束し、No.043と一致します。単純な例で解析解と照合することは、より複雑な尤度へ進む前の重要なテストです。実務システムでは最大反復回数、許容誤差、非有限値、目的関数の改善、未収束時の扱いをログに残します。

No.048：フィッシャー情報量 — 測定数と推定精度を結び付ける

実務での意味

測定を増やせば精度は上がりますが、費用も増えます。フィッシャー情報量は、観測が未知パラメータについて持つ情報を数量化し、検査設計や必要標本数の議論を支えます。

分析・モデル化の考え方

標準偏差 $\sigma$ が既知の正規平均では

$I_n(\mu)=\frac{n}{\sigma^2},\qquad \mathrm{SE}(\hat\mu)\approx\frac{1}{\sqrt{I_n}}=\frac{\sigma}{\sqrt n}$

です。標準誤差を半分にするには測定数が4倍必要です。独立同分布、モデル適合、既知パラメータという条件を外れると情報量の式も変わります。

Pythonで確認する

sample_sizes = np.array([20, 50, 100, 180, 400, 800])
fisher_info = sample_sizes / sigma_hat**2
standard_error = 1 / np.sqrt(fisher_info)
fisher_table = pd.DataFrame({
    "測定数 n": sample_sizes,
    "Fisher情報量": fisher_info,
    "平均推定の標準誤差 [mm]": standard_error,
    "近似95%半幅 [mm]": 1.96 * standard_error,
})
display(fisher_table.round(6))

plt.plot(sample_sizes, standard_error * 1000, marker="o")
plt.title("測定数と工程平均推定の標準誤差")
plt.xlabel("測定数 n [個]")
plt.ylabel("標準誤差 [μm]")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

	測定数 n	Fisher情報量	平均推定の標準誤差 [mm]	近似95%半幅 [mm]
0	20	10793.067110	0.009626	0.018866
1	50	26982.667774	0.006088	0.011932
2	100	53965.335549	0.004305	0.008437
3	180	97137.603988	0.003209	0.006289
4	400	215861.342195	0.002152	0.004219
5	800	431722.684389	0.001522	0.002983

png

結果の読み取り

測定数の増加に対して標準誤差は $1/\sqrt n$ でしか減らず、限界効果があります。必要精度から標本数を逆算し、測定費用や意思決定損失と比較するのが実務的です。自己相関がある連続測定では実効標本数が小さくなるため、この表をそのまま使えません。

No.049：ブートストラップ — 工程平均の不確かさを再標本化で評価する

実務での意味

複雑なKPIや理論分布を置きにくい統計量では、標準誤差の式を導くのが困難です。ブートストラップは観測標本から復元抽出を繰り返し、推定量のばらつきや区間を評価します。

分析・モデル化の考え方

元標本から同じサイズの標本を復元抽出し、各標本で平均を計算します。ここでは3,000回のパーセンタイル法で95%区間を求めます。観測標本が母集団を代表し、観測単位が交換可能であることが前提です。時系列やロット構造にはブロック単位の再標本化が必要です。

Pythonで確認する

n_boot = 3000
bootstrap_samples = rng.choice(diameter, size=(n_boot, len(diameter)), replace=True)
bootstrap_means = bootstrap_samples.mean(axis=1)
boot_ci = np.quantile(bootstrap_means, [0.025, 0.975])
theory_ci = mu_hat + np.array([-1, 1]) * stats.t.ppf(0.975, len(diameter)-1) * diameter.std(ddof=1) / np.sqrt(len(diameter))

ci_table = pd.DataFrame({
    "方法": ["Bootstrap percentile", "t近似"],
    "下限 [mm]": [boot_ci[0], theory_ci[0]],
    "上限 [mm]": [boot_ci[1], theory_ci[1]],
})
display(ci_table.round(5))

plt.hist(bootstrap_means, bins=30, density=True, alpha=0.7, edgecolor="white")
plt.axvline(boot_ci[0], color="tab:red", ls="--", label="Bootstrap 95%区間")
plt.axvline(boot_ci[1], color="tab:red", ls="--")
plt.axvline(mu_hat, color="black", lw=1.8, label="観測標本平均")
plt.title("ブートストラップ標本平均の分布")
plt.xlabel("再標本化した平均 [mm]")
plt.ylabel("確率密度")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	方法	下限 [mm]	上限 [mm]
0	Bootstrap percentile	20.01178	20.02422
1	t近似	20.01155	20.02425

png

結果の読み取り

ブートストラップ区間とt近似区間が近ければ、このデータでは方法選択による差が小さいと読めます。ただし、区間は工程平均の不確かさであり、個々の製品寸法が入る範囲ではありません。標本に含まれない異常モードや将来の工程変更は再現できないため、データ収集設計の代替にはなりません。

No.050：EMアルゴリズム — 設備IDのない混合工程を推定する

実務での意味

設備IDの欠損や合流工程により、異なる状態の測定値が混ざることがあります。全体平均だけでは、正常な2設備の差を「大きなばらつき」と誤認する可能性があります。EMアルゴリズムは、所属が見えないデータから混合比・平均・標準偏差を反復推定します。

分析・モデル化の考え方

2成分正規混合モデル

$p(x)=\sum_{k=1}^{2}\pi_k\,N(x\mid\mu_k,\sigma_k^2)$

を仮定します。E-stepで所属確率（responsibility）を計算し、M-stepでその確率を重みとして $\pi_k,\mu_k,\sigma_k$ を更新します。混合数は自動決定されず、局所解やラベル入替もあるため、複数初期値と業務知識が必要です。

Pythonで確認する

def em_two_normals(x, max_iter=200, tol=1e-8):
    means = np.quantile(x, [0.3, 0.7]).astype(float)
    sigmas = np.array([x.std(), x.std()])
    weights = np.array([0.5, 0.5])
    loglik_history = []

    for _ in range(max_iter):
        weighted_pdf = np.column_stack([
            weights[k] * stats.norm.pdf(x, means[k], sigmas[k]) for k in range(2)
        ])
        denominator = weighted_pdf.sum(axis=1, keepdims=True)
        responsibilities = weighted_pdf / denominator
        nk = responsibilities.sum(axis=0)
        weights = nk / len(x)
        means = (responsibilities * x[:, None]).sum(axis=0) / nk
        sigmas = np.sqrt(
            (responsibilities * (x[:, None] - means)**2).sum(axis=0) / nk
        )
        loglik = np.log(denominator[:, 0]).sum()
        loglik_history.append(loglik)
        if len(loglik_history) > 1 and abs(loglik_history[-1] - loglik_history[-2]) < tol:
            break

    order = np.argsort(means)
    return weights[order], means[order], sigmas[order], np.array(loglik_history), responsibilities[:, order]

em_weights, em_means, em_sigmas, em_loglik, responsibilities = em_two_normals(torque)
em_result = pd.DataFrame({
    "推定成分": ["低トルク側", "高トルク側"],
    "混合比": em_weights,
    "平均 [N·m]": em_means,
    "標準偏差 [N·m]": em_sigmas,
})
display(em_result.round(3))
print(f"反復回数: {len(em_loglik)}, 最終対数尤度: {em_loglik[-1]:.2f}")

x_grid = np.linspace(torque.min() - 0.5, torque.max() + 0.5, 400)
mixture_pdf = sum(
    em_weights[k] * stats.norm.pdf(x_grid, em_means[k], em_sigmas[k]) for k in range(2)
)
plt.hist(torque, bins=24, density=True, alpha=0.5, edgecolor="white", label="設備IDなし観測")
plt.plot(x_grid, mixture_pdf, color="black", lw=2.2, label="推定混合分布")
for k, color in enumerate(["tab:blue", "tab:orange"]):
    plt.plot(x_grid, em_weights[k] * stats.norm.pdf(x_grid, em_means[k], em_sigmas[k]),
             color=color, ls="--", label=f"成分{k+1}")
plt.title("EMアルゴリズムによる締付トルクの混合分布推定")
plt.xlabel("締付トルク [N·m]")
plt.ylabel("確率密度")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	推定成分	混合比	平均 [N·m]	標準偏差 [N·m]
0	低トルク側	0.657	48.692	0.731
1	高トルク側	0.343	52.592	1.071

反復回数: 34, 最終対数尤度: -589.01

png

結果の読み取り

全体を単一分布として扱うより、低トルク側と高トルク側の2成分として工程状態を説明できます。ただし、成分が実際の設備に対応するという結論は、保全記録や設備IDを復元して検証する必要があります。EMは因果を特定せず、初期値や成分数に依存します。まずトレーサビリティ改善を優先し、EMは欠損期間の補助分析として使います。

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

推定値と不確かさをセットで報告する：平均や率の一点だけでは、調整・保全・追加検査の優先順位を決められません。
データ型に発生機構を対応させる：連続量、合否、件数で尤度が異なり、必要な仮定も変わります。
解析解と数値解を相互検算する：単純なケースで一致を確認すると、複雑な実装の信頼性を高められます。
測定数は必要精度から決める：慣例的な抜取数ではなく、許容する推定誤差と意思決定コストから設計します。
混合分布はデータ管理上の問題を示すことがある：高度な推定で補う前に、設備・材料・作業条件の識別子を整備します。

統計モデルの出力は自動的な処置命令ではありません。工程知識、規格、損失、変更管理と結び付けて初めて意思決定に使えます。

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

目的と意思決定の定義：何を推定し、どの閾値で誰が何を判断するかを明文化する
測定システムの確認：校正、分解能、繰返し性・再現性を評価し、測定誤差を把握する
データ契約と層別キー：設備、金型、材料ロット、シフト、検査員、時刻を欠損なく連携する
モデル診断：独立性、定常性、分布適合、過分散、外れ値、混合状態を確認する
検証と監視：期間を分けた再現性確認、推定値・区間・収束ログ・データ品質の継続監視を行う
運用設計：異常時のエスカレーション、再学習、承認、ロールバック、監査証跡を定める

小規模な対象工程で、現行判断との並行運用から始めるのが安全です。推定精度だけでなく、見逃し・過剰調整・検査工数を含む業務KPIで評価します。

まとめ

No.041〜No.050では、正規・二項・ポアソンモデルの最尤推定を出発点に、SciPyによる数値最適化、尤度の形、勾配、ニュートン法、フィッシャー情報量、ブートストラップ、EMアルゴリズムまで実装しました。

中心となる考え方は、観測値に最も適合するパラメータを求めるだけでなく、仮定・推定精度・収束・データ生成過程を同時に確認することです。製造現場への導入では、モデル精度の改善と同じくらい、測定系、層別キー、判断ルール、監視体制の整備が重要です。

法人向けのご相談

数理工房では、製造業における品質データ分析、統計モデル構築、検査設計、設備保全分析、Python研修、PoCから運用定着までをご支援しています。「平均値の集計から一歩進めたい」「推定結果を現場の判断基準へ落とし込みたい」「設備ラベルやデータ品質に課題がある」といった段階からご相談いただけます。

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

観測データから工程の真の姿を推定する：製造業の統計的推定10本ノック

観測データから工程の真の姿を推定する：製造業の統計的推定10本ノック

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.041：最尤推定（正規分布）— 外径の工程中心とばらつきを推定する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.042：最尤推定（二項分布）— ロット横断の不適合率を推定する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.043：最尤推定（ポアソン分布）— 日次の突発停止発生率を推定する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.044：scipy.optimizeによる最適化 — 制約付きの数値推定を実装する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.045：尤度関数の可視化 — 推定値の周辺にある不確かさを見る

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.046：勾配計算 — 尤度が改善する方向を確認する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.047：ニュートン法 — 曲率を使って停止発生率を反復推定する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.048：フィッシャー情報量 — 測定数と推定精度を結び付ける

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.049：ブートストラップ — 工程平均の不確かさを再標本化で評価する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.050：EMアルゴリズム — 設備IDのない混合工程を推定する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

まとめ

法人向けのご相談

Related Posts

不確実性を設備・在庫・生産能力の判断へつなぐ：製造業シミュレーション10本ノック

需要変動を読み解き、生産計画へつなぐ：製造業の時系列分析10本ノック

不確かさを更新して改善案を選ぶ：製造業のベイズ統計10本ノック

製造業の回帰分析をPythonで実践｜線形回帰・GLM・GAM 10本ノック

改善効果を「偶然」で終わらせない：製造業の仮説検定10本ノック