100本ノック / 確率統計 / 確率・統計理論100本ノック
工程変更を「効いた気がする」で終わらせない――統計的推論で品質投資を判断する
工程変更を「効いた気がする」で終わらせない――統計的推論で品質投資を判断する
概要
製造現場では、新しい設備条件や材料へ切り替えた後に不良率が下がっても、それが施策の効果なのか偶然の揺れなのかを見分ける必要があります。本記事では、架空の精密部品工場における工程変更パイロットを題材に、仮説検定、尤度比・Wald・Score検定、p値、信頼区間、ベイズ推定・予測、AIC・BIC、統計学と機械学習の接点を、量産移行の意思決定へ結び付けます。
対象は「確率・統計理論100本ノック」の No.091〜No.100 です。単に「有意差あり／なし」を判定するのではなく、効果量、不確実性、将来不良数、モデルの複雑さ、予測性能をどう組み合わせるかをPythonで確認します。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
今回の題材は、寸法不良を減らすために加工条件を変更した精密部品工場です。従来条件と新条件を並行運転し、各ロットから100個を検査しました。品質責任者は次の問いに答える必要があります。
- 観測された不良率低下は、偶然の範囲を超えているか
- 統計的に差があっても、投資に値する大きさか
- 検定方法を変えると結論が変わらないか
- 次月2,000個では、不良が何個発生し得るか
- 温度、加工速度、材料供給元を考慮しても新条件の効果は残るか
- 説明を重視する統計モデルと、予測を重視する機械学習をどう使い分けるか
統計的推論は、標本から母集団について断定する技術ではありません。仮定とデータのもとで、判断に残る不確実性を数値化する技術です。
現場でよくある状況
工程変更の評価では、次のような行き違いが起こります。
- 不良率が下がったという一点だけで量産移行を決める
- p値が0.05未満なら効果が大きい、と解釈する
- p値が0.05以上なら「効果がない」と結論する
- 同じデータを何度も切り分け、有意だった比較だけを報告する
- 温度や材料構成が違うのに、単純な新旧比較だけで因果効果とみなす
- 学習データへの当てはまりだけで予測モデルを選ぶ
これらを避けるには、試験前の仮説、評価指標、許容誤り、実務上の最小効果、分析方法、停止条件を決めておく必要があります。
なぜこの問題は判断が難しいのか
第一に、標本不良率は毎回揺れます。母不良率が同じでも検査個数が有限なら新旧の差はゼロになりません。したがって「差があること」と「偶然では説明しにくいこと」を分けます。
第二に、検定には2種類の誤りがあります。実際には効果がないのに採用する第一種の誤りと、実際には効果があるのに見送る第二種の誤りです。品質流出の損失と改善機会の損失が非対称なら、慣例的な有意水準だけでなく事業損失も考える必要があります。
第三に、同じデータでも問いが違えば評価量が違います。頻度論は繰り返し標本抽出に基づく誤差管理、ベイズ推定は事前情報とデータを統合した確率更新、機械学習は未知データへの予測性能を中心に据えます。目的を混ぜず、必要に応じて組み合わせることが重要です。
今回扱うノックの全体像
| No. | テーマ | 製造業での問い |
|---|---|---|
| 091 | 仮説検定とは | 新条件の不良率低下は偶然で説明できるか |
| 092 | 尤度比検定 | 新旧で別の不良率を持つモデルは必要か |
| 093 | Wald検定 | 推定した不良率差は標準誤差の何倍か |
| 094 | Score検定 | 帰無仮説のもとで観測差はどれほど極端か |
| 095 | p値の意味 | p値が表すもの・表さないものは何か |
| 096 | 信頼区間 | 改善幅として整合的な範囲はどこか |
| 097 | ベイズ推定 | 新条件の不良率と改善確率をどう更新するか |
| 098 | ベイズ予測分布 | 次月の不良個数をどの範囲で見込むか |
| 099 | モデル選択（AIC・BIC） | 条件差を説明する複雑さはどこまで必要か |
| 100 | 統計学と機械学習の接点 | 効果説明と個体予測をどう役割分担するか |
No.091〜No.096で頻度論的推論、No.097〜No.098でベイズ推論、No.099〜No.100で説明モデルと予測モデルの選択へ進みます。
Python 環境の準備
NumPyで乱数生成と数値計算、pandasで集計、SciPyで確率分布と最適化、matplotlibで可視化、scikit-learnで未知データへの予測評価を行います。グラフの日本語表示にはjapanize-matplotlibを使います。再実行して同じ結果になるよう乱数seedは固定します。
import sys
import numpy as np
import pandas as pd
import scipy
from scipy.optimize import minimize
from scipy.special import expit, betaln, gammaln
from scipy.stats import beta, chi2, norm
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib
import sklearn
from IPython.display import display
pd.set_option("display.max_columns", 20)
pd.set_option("display.precision", 4)
plt.rcParams["figure.figsize"] = (8, 4.5)
print(f"Python : {sys.version.split()[0]}")
print(f"NumPy : {np.__version__}")
print(f"pandas : {pd.__version__}")
print(f"SciPy : {scipy.__version__}")
print(f"matplotlib : {matplotlib.__version__}")
print(f"scikit-learn: {sklearn.__version__}")Python : 3.13.1
NumPy : 2.5.1
pandas : 3.0.3
SciPy : 1.18.0
matplotlib : 3.11.0
scikit-learn: 1.9.0架空データの作成
180ロットについて従来条件・新条件を並行運転し、各ロット100個を検査したデータを生成します。
- ロットごとに室温、加工速度、材料供給元が異なる
- 高温、高速加工、供給元Cは不良確率を押し上げる
- 新条件は平均的に不良確率を下げ、高速時の悪化も一部緩和する
- 個体データとロット集計データの両方を保持する
新旧の割り付けは架空の並行試験を模しています。実務では、時期、設備、作業者、品種が偏らない無作為化やブロック化が、分析手法以上に重要です。
SEED = 20260711
rng = np.random.default_rng(SEED)
n_lots, units_per_lot = 180, 100
lot_id = np.arange(1, n_lots + 1)
new_process = rng.integers(0, 2, n_lots)
temperature = rng.normal(25.0, 3.0, n_lots)
speed = np.clip(rng.normal(100.0, 5.0, n_lots), 88, 112)
supplier_c = rng.binomial(1, 0.30, n_lots)
logit_p = (
-3.00
- 0.52 * new_process
+ 0.075 * (temperature - 25)
+ 0.040 * (speed - 100)
+ 0.42 * supplier_c
- 0.025 * new_process * (speed - 100)
)
defect_probability = expit(logit_p)
defect_matrix = rng.binomial(1, defect_probability[:, None], (n_lots, units_per_lot))
lots = pd.DataFrame({
"ロットID": lot_id,
"工程条件": np.where(new_process == 1, "新条件", "従来条件"),
"新条件フラグ": new_process,
"室温_C": temperature,
"加工速度_piece_per_min": speed,
"供給元Cフラグ": supplier_c,
"検査数": units_per_lot,
"不良数": defect_matrix.sum(axis=1),
})
lots["不良率"] = lots["不良数"] / lots["検査数"]
unit_df = lots.loc[lots.index.repeat(units_per_lot), [
"ロットID", "工程条件", "新条件フラグ", "室温_C",
"加工速度_piece_per_min", "供給元Cフラグ"
]].reset_index(drop=True)
unit_df["不良"] = defect_matrix.reshape(-1)
process_summary = lots.groupby("工程条件", sort=False).agg(
ロット数=("ロットID", "size"), 検査数=("検査数", "sum"),
不良数=("不良数", "sum"), 平均室温_C=("室温_C", "mean"),
平均加工速度=("加工速度_piece_per_min", "mean")
)
process_summary["不良率"] = process_summary["不良数"] / process_summary["検査数"]
display(lots.head())
display(process_summary)
fig, ax = plt.subplots()
plot_order = ["従来条件", "新条件"]
lot_plot = [lots.loc[lots["工程条件"] == x, "不良率"] * 100 for x in plot_order]
ax.boxplot(lot_plot, tick_labels=plot_order, showmeans=True)
ax.set_title("工程条件別のロット不良率")
ax.set_xlabel("工程条件")
ax.set_ylabel("ロット不良率（%）")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()| ロットID | 工程条件 | 新条件フラグ | 室温_C | 加工速度_piece_per_min | 供給元Cフラグ | 検査数 | 不良数 | 不良率 | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 新条件 | 1 | 26.1386 | 104.2447 | 1 | 100 | 3 | 0.03 |
| 1 | 2 | 従来条件 | 0 | 25.2721 | 107.1965 | 0 | 100 | 5 | 0.05 |
| 2 | 3 | 従来条件 | 0 | 26.6385 | 103.8079 | 0 | 100 | 5 | 0.05 |
| 3 | 4 | 新条件 | 1 | 21.5988 | 91.5556 | 1 | 100 | 0 | 0.00 |
| 4 | 5 | 従来条件 | 0 | 24.5195 | 106.7708 | 0 | 100 | 8 | 0.08 |
| ロット数 | 検査数 | 不良数 | 平均室温_C | 平均加工速度 | 不良率 | |
|---|---|---|---|---|---|---|
| 工程条件 | ||||||
| 新条件 | 91 | 9100 | 325 | 25.1617 | 100.2791 | 0.0357 |
| 従来条件 | 89 | 8900 | 498 | 24.7362 | 101.0576 | 0.0560 |
No.091：仮説検定とは
実務での意味
仮説検定は、観測された不良率差が「工程条件に差がない」という基準からどれほど外れているかを測ります。量産採用の自動判定器ではなく、偶然変動に対する証拠の強さをそろえた手順で評価する道具です。
分析・モデル化の考え方
従来条件と新条件の母不良率を とし、片側仮説
を置きます。有意水準を とし、帰無仮説のもとで標準化した統計量が左側5%の棄却域に入るかを確認します。方向をデータ確認後に選ぶと誤り率が崩れるため、片側か両側かは試験前に決めます。
Pythonで確認する
agg = unit_df.groupby("工程条件")["不良"].agg(["sum", "count"])
x0, n0 = agg.loc["従来条件", ["sum", "count"]]
x1, n1 = agg.loc["新条件", ["sum", "count"]]
p0_hat, p1_hat = x0 / n0, x1 / n1
pooled = (x0 + x1) / (n0 + n1)
se_null = np.sqrt(pooled * (1 - pooled) * (1 / n0 + 1 / n1))
z_score = (p1_hat - p0_hat) / se_null
p_one_sided = norm.cdf(z_score)
critical = norm.ppf(0.05)
test_091 = pd.DataFrame({
"観測差_新-従来": [p1_hat - p0_hat], "z統計量": [z_score],
"片側p値": [p_one_sided], "5%臨界値": [critical],
"判定": ["帰無仮説を棄却" if z_score < critical else "棄却できない"]
})
display(test_091)
z_grid = np.linspace(-4, 4, 500)
fig, ax = plt.subplots()
ax.plot(z_grid, norm.pdf(z_grid), label="帰無仮説下の標準正規分布")
ax.fill_between(z_grid, 0, norm.pdf(z_grid), where=z_grid <= critical,
alpha=0.3, color="tab:red", label="棄却域（5%）")
ax.axvline(z_score, color="black", linestyle="--", label=f"観測 z={z_score:.2f}")
ax.set_title("片側検定の棄却域と観測統計量")
ax.set_xlabel("z統計量")
ax.set_ylabel("確率密度")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 観測差_新-従来 | z統計量 | 片側p値 | 5%臨界値 | 判定 | |
|---|---|---|---|---|---|
| 0 | -0.0202 | -6.4999 | 4.0192e-11 | -1.6449 | 帰無仮説を棄却 |
結果の読み取り
観測された新条件の不良率は従来条件より低く、z統計量は片側5%の棄却域に入ります。このデータと仮定のもとでは「新旧の母不良率が同じ」では説明しにくい差です。ただし、棄却は新条件が原因だと単独で証明しません。割付の公平性、測定方法、期間差、途中停止の有無を確認したうえで、改善幅と費用を次の判断へつなげます。
No.092：尤度比検定
実務での意味
尤度比検定は、「新旧で共通の不良率しか持たない単純なモデル」と「新旧別の不良率を持つモデル」のデータへの適合度を比較します。工程条件をモデルへ追加する価値があるか、というモデル比較として読めます。
分析・モデル化の考え方
二項モデルの対数尤度を とすると、検定統計量は
です。大標本では のもとで自由度1のカイ二乗分布に近づきます。尤度は「仮説が正しい確率」ではなく、パラメータを固定したときに観測データがどれほど整合するかを表します。
Pythonで確認する
def binomial_loglik(x, n, p):
p = np.clip(p, 1e-12, 1 - 1e-12)
return x * np.log(p) + (n - x) * np.log1p(-p)
ll_null = binomial_loglik(x0 + x1, n0 + n1, pooled)
ll_alt = binomial_loglik(x0, n0, p0_hat) + binomial_loglik(x1, n1, p1_hat)
lr_stat = 2 * (ll_alt - ll_null)
lr_p = chi2.sf(lr_stat, df=1)
lr_table = pd.DataFrame({
"モデル": ["共通不良率", "新旧別不良率"],
"パラメータ数": [1, 2], "最大対数尤度": [ll_null, ll_alt]
})
display(lr_table)
display(pd.DataFrame({"尤度比統計量": [lr_stat], "自由度": [1], "p値": [lr_p]}))
fig, ax = plt.subplots()
g_grid = np.linspace(0, max(12, lr_stat * 1.15), 500)
ax.plot(g_grid, chi2.pdf(g_grid, 1), label=r"$\chi^2(1)$")
ax.axvline(lr_stat, color="tab:red", linestyle="--", label=f"観測 $G^2$={lr_stat:.2f}")
ax.set_title("尤度比統計量の基準分布")
ax.set_xlabel("尤度比統計量 $G^2$")
ax.set_ylabel("確率密度")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| モデル | パラメータ数 | 最大対数尤度 | |
|---|---|---|---|
| 0 | 共通不良率 | 1 | -3342.9874 |
| 1 | 新旧別不良率 | 2 | -3321.7292 |
| 尤度比統計量 | 自由度 | p値 | |
|---|---|---|---|
| 0 | 42.5165 | 1 | 7.0089e-11 |
結果の読み取り
新旧別モデルは共通不良率モデルより対数尤度が高く、尤度比統計量に対応するp値も小さくなります。工程条件という説明変数を追加する統計的根拠があります。ただし、適合度の改善は経済価値とは別です。新条件の設備費、速度、歩留まり、保証損失を含む損益分岐点と合わせて採否を決めます。
No.093：Wald検定
実務での意味
Wald検定は、推定した改善幅がその標準誤差に比べて十分大きいかを評価します。推定値と標準誤差があれば計算しやすいため、回帰モデルの係数評価にも広く使われます。
分析・モデル化の考え方
不良率差 のWald統計量を
とします。分母は制約を置かない推定値で評価します。標本が小さい、比率が0や1に近い、モデルが境界に近い場合は近似が悪くなることがあります。
Pythonで確認する
difference = p1_hat - p0_hat
se_wald = np.sqrt(p1_hat * (1 - p1_hat) / n1 + p0_hat * (1 - p0_hat) / n0)
z_wald = difference / se_wald
p_wald_two = 2 * norm.sf(abs(z_wald))
components = pd.DataFrame({
"項目": ["新条件の分散寄与", "従来条件の分散寄与", "標準誤差"],
"値": [p1_hat * (1 - p1_hat) / n1, p0_hat * (1 - p0_hat) / n0, se_wald]
})
display(components)
display(pd.DataFrame({
"不良率差": [difference], "Wald z": [z_wald], "両側p値": [p_wald_two]
}))
fig, ax = plt.subplots()
labels = ["従来条件", "新条件"]
rates = np.array([p0_hat, p1_hat]) * 100
errors = np.array([
np.sqrt(p0_hat * (1 - p0_hat) / n0),
np.sqrt(p1_hat * (1 - p1_hat) / n1)
]) * 1.96 * 100
ax.errorbar(labels, rates, yerr=errors, fmt="o", capsize=6, markersize=8)
ax.set_title("工程条件別不良率のWald近似区間")
ax.set_xlabel("工程条件")
ax.set_ylabel("不良率（%）")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()| 項目 | 値 | |
|---|---|---|
| 0 | 新条件の分散寄与 | 3.7845e-06 |
| 1 | 従来条件の分散寄与 | 5.9353e-06 |
| 2 | 標準誤差 | 3.1177e-03 |
| 不良率差 | Wald z | 両側p値 | |
|---|---|---|---|
| 0 | -0.0202 | -6.4923 | 8.4532e-11 |
結果の読み取り
推定された不良率差は標準誤差の複数倍あり、両側検定でも差を示す結果になります。Wald検定は推定値を中心に見るため、効果量との対応が分かりやすい一方、希少不良や小標本では不安定です。その場合は正確検定、尤度法、プロファイル尤度区間なども候補にします。
No.094：Score検定
実務での意味
Score検定は、差がないと仮定した位置で、尤度が差のある方向へどれだけ強く傾いているかを調べます。代替モデル全体を複雑に推定せずに変数追加を検討できる場面があります。
分析・モデル化の考え方
一般にスコアは対数尤度の傾き
です。2比率の同等性では、帰無仮説下の共通推定値 を標準誤差に使うz検定がScore検定に対応します。Waldは非制約推定値、Scoreは帰無仮説下の推定値、尤度比は両者の最大尤度の差を見る点が異なります。
Pythonで確認する
score_z = difference / np.sqrt(pooled * (1 - pooled) * (1 / n0 + 1 / n1))
score_chi2 = score_z ** 2
score_p = chi2.sf(score_chi2, 1)
comparison = pd.DataFrame({
"検定": ["尤度比", "Wald", "Score"],
"統計量（カイ二乗尺度）": [lr_stat, z_wald ** 2, score_chi2],
"両側p値": [lr_p, p_wald_two, score_p],
"評価位置": ["制約・非制約の尤度差", "非制約推定値", "帰無仮説下"]
})
display(comparison)
fig, ax = plt.subplots()
ax.bar(comparison["検定"], comparison["統計量（カイ二乗尺度）"],
color=["tab:blue", "tab:orange", "tab:green"])
ax.axhline(chi2.ppf(0.95, 1), color="tab:red", linestyle="--", label="5%臨界値")
ax.set_title("3つの大標本検定の比較")
ax.set_xlabel("検定方法")
ax.set_ylabel("カイ二乗尺度の統計量")
ax.grid(axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 検定 | 統計量（カイ二乗尺度） | 両側p値 | 評価位置 | |
|---|---|---|---|---|
| 0 | 尤度比 | 42.5165 | 7.0089e-11 | 制約・非制約の尤度差 |
| 1 | Wald | 42.1500 | 8.4532e-11 | 非制約推定値 |
| 2 | Score | 42.2485 | 8.0383e-11 | 帰無仮説下 |
結果の読み取り
今回は検査数が多く不良率も境界から離れているため、尤度比・Wald・Scoreの3検定は近い結論になります。一致は常に保証されません。値が食い違う場合は多数決にせず、小標本、境界、モデル誤指定、数値最適化を点検し、データ生成過程に適した方法を選びます。
No.095：p値の意味
実務での意味
p値は「新条件に効果がない確率」でも「今回の結論が間違う確率」でもありません。帰無仮説と分析手順を固定したときに、観測値以上に極端な結果が出る確率です。経営判断では効果量と損失を別に評価します。
分析・モデル化の考え方
統計量を 、観測値を とすると、両側p値は概念的に
です。帰無仮説が正しく、連続な検定が正しく運用されるならp値は概ね一様分布になります。任意の時点で繰り返し覗いて止める、多重比較から小さい値だけ選ぶ、といった運用ではこの性質が崩れます。
Pythonで確認する
rng_p = np.random.default_rng(SEED + 95)
n_sim = 10_000
sim_x0 = rng_p.binomial(n0, pooled, n_sim)
sim_x1 = rng_p.binomial(n1, pooled, n_sim)
sim_d = sim_x1 / n1 - sim_x0 / n0
sim_z = sim_d / se_null
simulation_p = np.mean(np.abs(sim_z) >= abs(score_z))
display(pd.DataFrame({
"理論上の両側p値": [score_p], "帰無分布シミュレーションp値": [simulation_p],
"シミュレーション回数": [n_sim]
}))
fig, ax = plt.subplots()
ax.hist(sim_z, bins=50, density=True, alpha=0.7, color="tab:blue", label="帰無仮説下の模擬z")
ax.axvline(score_z, color="tab:red", linestyle="--", label=f"観測 z={score_z:.2f}")
ax.axvline(-score_z, color="tab:red", linestyle="--")
ax.set_title("p値を作る帰無仮説下の反復分布")
ax.set_xlabel("z統計量")
ax.set_ylabel("確率密度")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 理論上の両側p値 | 帰無分布シミュレーションp値 | シミュレーション回数 | |
|---|---|---|---|
| 0 | 8.0383e-11 | 0.0 | 10000 |
結果の読み取り
帰無仮説下で新旧に同じ不良率を設定した反復では、観測統計量以上に極端な値はほとんど現れません。これが小さいp値の意味です。ただし、p値は改善額や再現確率を直接与えません。試験前登録、効果量、信頼区間、測定品質、追試を併記し、「p値が小さい」だけの承認資料にしないことが重要です。
No.096：信頼区間
実務での意味
信頼区間は、改善幅としてデータと整合的な範囲を示します。点推定だけでなく、悲観側でも投資基準を満たすかを見ることで、量産採用の判断が実務的になります。
分析・モデル化の考え方
不良率差の近似95%信頼区間は
です。頻度論的には、同じ手順を繰り返したときに作られる区間の約95%が真値を覆う、という手続きの性質です。今回得た固定区間に真値が入る確率が95%という意味ではありません。
Pythonで確認する
z975 = norm.ppf(0.975)
ci_low, ci_high = difference - z975 * se_wald, difference + z975 * se_wald
# 月産100,000個、流出不良1個あたり8,000円という架空条件で金額換算
monthly_volume, loss_per_defect = 100_000, 8_000
avoided_defects_range = (-ci_high * monthly_volume, -ci_low * monthly_volume)
avoided_loss_range = tuple(x * loss_per_defect for x in avoided_defects_range)
ci_table = pd.DataFrame({
"指標": ["不良率差（新-従来）", "月間回避不良個数", "月間回避損失_円"],
"点推定": [difference, -difference * monthly_volume,
-difference * monthly_volume * loss_per_defect],
"95%下限": [ci_low, avoided_defects_range[0], avoided_loss_range[0]],
"95%上限": [ci_high, avoided_defects_range[1], avoided_loss_range[1]]
})
display(ci_table)
fig, ax = plt.subplots()
ax.errorbar([difference * 100], [0],
xerr=[[difference * 100 - ci_low * 100], [ci_high * 100 - difference * 100]],
fmt="o", capsize=7, markersize=8)
ax.axvline(0, color="tab:red", linestyle="--", label="差なし")
ax.set_title("新条件と従来条件の不良率差：95%信頼区間")
ax.set_xlabel("不良率差（新条件 − 従来条件、%ポイント）")
ax.set_ylabel("推定結果")
ax.set_yticks([])
ax.grid(axis="x", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 指標 | 点推定 | 95%下限 | 95%上限 | |
|---|---|---|---|---|
| 0 | 不良率差（新-従来） | -2.0241e-02 | -2.6351e-02 | -1.4130e-02 |
| 1 | 月間回避不良個数 | 2.0241e+03 | 1.4130e+03 | 2.6351e+03 |
| 2 | 月間回避損失_円 | 1.6193e+07 | 1.1304e+07 | 2.1081e+07 |
結果の読み取り
不良率差の95%信頼区間はゼロより下側にあり、改善方向と整合します。さらに架空の生産量と損失単価で、回避不良個数と損失額の幅へ翻訳できます。実務ではこの悲観側の効果と導入費を比較します。なお、損失単価や将来の製品構成にも不確実性があるため、金額区間は統計区間だけで完結しません。
No.097：ベイズ推定
実務での意味
ベイズ推定は、不良率についての事前情報を今回の検査結果で更新し、「不良率が何%以下である確率」や「新条件が従来条件より良い確率」を直接計算します。過去試験や類似設備の情報を明示的に統合できます。
分析・モデル化の考え方
不良確率 にBeta事前分布、観測不良数 に二項分布を置くと
となります。ここでは弱情報のJeffreys事前分布 を両条件に使います。事前分布は隠さず、感度分析の対象にします。
Pythonで確認する
a_prior = b_prior = 0.5
post0 = (a_prior + x0, b_prior + n0 - x0)
post1 = (a_prior + x1, b_prior + n1 - x1)
rng_b = np.random.default_rng(SEED + 97)
draws0 = rng_b.beta(*post0, 100_000)
draws1 = rng_b.beta(*post1, 100_000)
posterior_table = pd.DataFrame({
"工程条件": ["従来条件", "新条件"],
"事後平均": [draws0.mean(), draws1.mean()],
"95%信用区間下限": [np.quantile(draws0, 0.025), np.quantile(draws1, 0.025)],
"95%信用区間上限": [np.quantile(draws0, 0.975), np.quantile(draws1, 0.975)]
})
display(posterior_table)
display(pd.DataFrame({
"P(新条件の不良率 < 従来条件)": [np.mean(draws1 < draws0)],
"P(1%ポイント以上改善)": [np.mean(draws0 - draws1 > 0.01)]
}))
p_grid = np.linspace(0.01, 0.08, 500)
fig, ax = plt.subplots()
ax.plot(p_grid * 100, beta.pdf(p_grid, *post0), label="従来条件")
ax.plot(p_grid * 100, beta.pdf(p_grid, *post1), label="新条件")
ax.set_title("工程条件別の不良率の事後分布")
ax.set_xlabel("母不良率（%）")
ax.set_ylabel("事後確率密度")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 工程条件 | 事後平均 | 95%信用区間下限 | 95%信用区間上限 | |
|---|---|---|---|---|
| 0 | 従来条件 | 0.0560 | 0.0513 | 0.0609 |
| 1 | 新条件 | 0.0358 | 0.0320 | 0.0397 |
| P(新条件の不良率 < 従来条件) | P(1%ポイント以上改善) | |
|---|---|---|
| 0 | 1.0 | 0.9994 |
結果の読み取り
新条件の事後分布は従来条件より低い側へ位置し、新条件の母不良率が低い事後確率を直接確認できます。「1%ポイント以上改善」の確率は、実務上の最小効果を満たす確からしさとして読めます。ただし、無作為化不足や測定バイアスは事後分布を狭くしても消えません。モデルとデータ収集の妥当性が前提です。
No.098：ベイズ予測分布
実務での意味
品質部門が知りたいのは母不良率だけでなく、次月に不良が何個発生し得るかです。ベイズ予測分布は、二項変動とパラメータ推定の不確実性をまとめ、検査要員、手直し能力、予備費の設計へつなげます。
分析・モデル化の考え方
将来の不良数を とすると、事後予測分布は
です。Beta-Binomialモデルでは、事後分布から不良率を引き、その不良率で将来の二項乱数を生成すれば予測できます。点推定を固定した二項予測より、パラメータ不確実性の分だけ裾が広くなります。
Pythonで確認する
future_n = 2_000
rng_pred = np.random.default_rng(SEED + 98)
posterior_p = rng_pred.beta(*post1, 100_000)
predictive_defects = rng_pred.binomial(future_n, posterior_p)
plug_in_defects = rng_pred.binomial(future_n, p1_hat, 100_000)
pred_table = pd.DataFrame({
"方法": ["ベイズ事後予測", "点推定固定の二項予測"],
"平均不良数": [predictive_defects.mean(), plug_in_defects.mean()],
"2.5%点": [np.quantile(predictive_defects, 0.025), np.quantile(plug_in_defects, 0.025)],
"97.5%点": [np.quantile(predictive_defects, 0.975), np.quantile(plug_in_defects, 0.975)],
"標準偏差": [predictive_defects.std(), plug_in_defects.std()]
})
display(pred_table)
fig, ax = plt.subplots()
bins = np.arange(min(predictive_defects.min(), plug_in_defects.min()),
max(predictive_defects.max(), plug_in_defects.max()) + 2) - 0.5
ax.hist(plug_in_defects, bins=bins, density=True, alpha=0.55, label="点推定固定")
ax.hist(predictive_defects, bins=bins, density=True, alpha=0.55, label="ベイズ事後予測")
ax.set_title("新条件で次の2,000個を生産したときの不良数予測")
ax.set_xlabel("将来不良数")
ax.set_ylabel("予測確率密度")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 方法 | 平均不良数 | 2.5%点 | 97.5%点 | 標準偏差 | |
|---|---|---|---|---|---|
| 0 | ベイズ事後予測 | 71.5155 | 54.0 | 90.0 | 9.1363 |
| 1 | 点推定固定の二項予測 | 71.4918 | 56.0 | 88.0 | 8.3184 |
結果の読み取り
事後予測の平均は新条件の推定不良率から想定する個数に近い一方、点推定固定の予測より幅が広くなります。これは不良発生の偶然だけでなく、母不良率を有限データから推定した不確実性も含むためです。人員や予算を平均だけで組まず、たとえば97.5%点を高負荷シナリオとして使えます。
No.099：モデル選択（AIC・BIC）
実務での意味
不良率は工程条件だけでなく、室温、加工速度、材料供給元にも左右されます。しかし変数や交互作用を増やせば学習データへの適合は必ず改善します。AIC・BICは適合度と複雑さを同じ尺度で比較します。
分析・モデル化の考え方
最大対数尤度を 、パラメータ数を 、観測数を とすると
です。小さいモデルを選びます。AICは予測損失の近似、BICは候補の中に真のモデルがある等の仮定のもとでより強く複雑さを罰します。異なる目的変数や異なるデータで計算した値は比較できません。
Pythonで確認する
y = unit_df["不良"].to_numpy(dtype=float)
new = unit_df["新条件フラグ"].to_numpy(dtype=float)
temp_c = unit_df["室温_C"].to_numpy() - 25
speed_c = unit_df["加工速度_piece_per_min"].to_numpy() - 100
supplier = unit_df["供給元Cフラグ"].to_numpy(dtype=float)
design_matrices = {
"M1: 工程条件のみ": np.column_stack([np.ones(len(y)), new]),
"M2: 主効果": np.column_stack([np.ones(len(y)), new, temp_c, speed_c, supplier]),
"M3: 主効果+速度交互作用": np.column_stack([
np.ones(len(y)), new, temp_c, speed_c, supplier, new * speed_c
]),
}
def fit_logistic(X, y):
def negative_loglik(beta_coef):
eta = X @ beta_coef
return np.sum(np.logaddexp(0, eta) - y * eta)
result = minimize(negative_loglik, np.zeros(X.shape[1]), method="L-BFGS-B")
return result.x, -result.fun, result.success
model_rows = []
fitted_models = {}
for name, X in design_matrices.items():
coef, loglik, success = fit_logistic(X, y)
k = X.shape[1]
fitted_models[name] = coef
model_rows.append({
"モデル": name, "パラメータ数": k, "最大対数尤度": loglik,
"AIC": 2 * k - 2 * loglik,
"BIC": k * np.log(len(y)) - 2 * loglik,
"収束": success
})
model_comparison = pd.DataFrame(model_rows).sort_values("AIC")
display(model_comparison)
fig, ax = plt.subplots()
x_pos = np.arange(len(model_comparison))
width = 0.36
ax.bar(x_pos - width / 2, model_comparison["AIC"], width, label="AIC")
ax.bar(x_pos + width / 2, model_comparison["BIC"], width, label="BIC")
ax.set_xticks(x_pos, model_comparison["モデル"], rotation=12)
ax.set_title("候補ロジスティック回帰モデルのAIC・BIC")
ax.set_xlabel("候補モデル")
ax.set_ylabel("情報量規準（小さいほど良い）")
ax.grid(axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| モデル | パラメータ数 | 最大対数尤度 | AIC | BIC | 収束 | |
|---|---|---|---|---|---|---|
| 1 | M2: 主効果 | 5 | -3274.5749 | 6559.1499 | 6598.1405 | True |
| 2 | M3: 主効果+速度交互作用 | 6 | -3274.1920 | 6560.3839 | 6607.1727 | True |
| 0 | M1: 工程条件のみ | 2 | -3321.7292 | 6647.4583 | 6663.0546 | True |
結果の読み取り
工程条件だけのモデルより、操業条件を加えたモデルの情報量規準が改善します。交互作用を追加する価値はAICとBICで評価が分かれる場合があり、目的に応じた判断が必要です。AIC・BICの差だけで因果変数を決めず、階層原則、現場知識、残差診断、外部期間での再現性も確認します。
No.100：統計学と機械学習の接点
実務での意味
統計学は工程変更の効果や不確実性の説明に、機械学習は次の個体・ロットのリスク順位付けに強みがあります。品質改善では「なぜ変わったか」と「どれが危ないか」の両方が必要です。
分析・モデル化の考え方
両者に明確な境界があるわけではなく、ロジスティック回帰は統計推論にも機械学習にも使われます。違いは主に評価目的です。
- 推論：係数、効果量、信頼区間、仮定、識別可能性を重視する
- 予測：未使用データの対数損失、Brier score、ROC-AUC、校正を重視する
不良は少ないため、正解率だけでは「すべて良品」とするモデルを高評価しかねません。確率予測の鋭さと校正を分離して確認します。
Pythonで確認する
from sklearn.calibration import calibration_curve
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, brier_score_loss, log_loss, roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
feature_cols = ["新条件フラグ", "室温_C", "加工速度_piece_per_min", "供給元Cフラグ"]
X_ml = unit_df[feature_cols]
y_ml = unit_df["不良"]
X_train, X_test, y_train, y_test = train_test_split(
X_ml, y_ml, test_size=0.30, random_state=SEED, stratify=y_ml
)
predictors = {
"ロジスティック回帰": LogisticRegression(max_iter=1_000),
"決定木（深さ4）": DecisionTreeClassifier(max_depth=4, min_samples_leaf=150, random_state=SEED),
}
metric_rows, probability_predictions = [], {}
for name, estimator in predictors.items():
estimator.fit(X_train, y_train)
probability = estimator.predict_proba(X_test)[:, 1]
probability_predictions[name] = probability
metric_rows.append({
"モデル": name,
"Accuracy": accuracy_score(y_test, probability >= 0.5),
"ROC-AUC": roc_auc_score(y_test, probability),
"Log loss": log_loss(y_test, probability),
"Brier score": brier_score_loss(y_test, probability),
})
metrics = pd.DataFrame(metric_rows)
display(metrics)
print(f"テストデータの不良率: {y_test.mean():.4f}")
fig, ax = plt.subplots()
for name, probability in probability_predictions.items():
observed, predicted = calibration_curve(y_test, probability, n_bins=6, strategy="quantile")
ax.plot(predicted * 100, observed * 100, marker="o", label=name)
max_rate = max(ax.get_xlim()[1], ax.get_ylim()[1])
ax.plot([0, max_rate], [0, max_rate], linestyle="--", color="black", label="完全な校正")
ax.set_title("未知データにおける不良確率の校正")
ax.set_xlabel("予測不良率（%）")
ax.set_ylabel("実測不良率（%）")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| モデル | Accuracy | ROC-AUC | Log loss | Brier score | |
|---|---|---|---|---|---|
| 0 | ロジスティック回帰 | 0.9543 | 0.6156 | 0.1819 | 0.0433 |
| 1 | 決定木（深さ4） | 0.9543 | 0.5799 | 0.1835 | 0.0434 |
テストデータの不良率: 0.0457
結果の読み取り
不良が少ないため両モデルのAccuracyは高く見えますが、それだけではモデルの価値を比較できません。ROC-AUCは順位付け、Log lossとBrier scoreは確率予測、校正曲線は「予測5%の群で本当に約5%不良か」を評価します。効果説明には試験設計と推論モデル、検査優先順位には外部検証済みの予測モデルを使い、最終判断では誤検知・見逃しの費用を閾値へ反映します。
対象ノックを通して見える実務上の示唆
No.091〜No.100を通して、工程変更の評価を次の4層で設計できることが分かります。
- 証拠をそろえる：仮説、片側・両側、有意水準、検査数を試験前に決める
- 大きさを読む：p値だけでなく、効果量と信頼区間・信用区間を示す
- 将来へ翻訳する：予測分布を不良個数、人員、損失額、供給リスクへ変換する
- 目的別にモデルを選ぶ：説明には仮定と不確実性、予測には未知データ性能と校正を重視する
尤度比・Wald・Scoreは同じ問いを異なる位置から評価します。大標本で近い結果になることは安心材料ですが、検定の一致が試験設計上の偏りを修正するわけではありません。また、頻度論とベイズ推論は敵対する選択肢ではなく、誤り率管理と意思決定確率という異なる情報を提供します。
経営会議では「有意差あり」だけでなく、不良率差、区間、最小実務効果を満たす確率、次期不良数の予測範囲、金額換算、前提条件を1枚にまとめると、技術結果を投資判断へ接続しやすくなります。
実務導入する場合に必要なこと
1. 試験計画を分析前に固定する
主評価指標、帰無仮説、片側・両側、有意水準、必要検査数、除外基準、途中解析、最小実務効果を計画書へ記載します。複数KPIや多数条件を試す場合は、多重性を管理します。
2. 比較可能性を工程側で作る
新旧条件を時期、設備、品種、材料、作業者へ偏らせず、可能なら無作為化・ブロック化します。同一の測定器、検査基準、サンプリング方法を使い、測定システム解析も行います。
3. 効果量を損益へつなげる
不良率差を廃棄、手直し、選別、流出保証、納期、能力損失へ換算します。導入費、保全費、速度低下、副作用を含め、悲観・標準・楽観シナリオで投資回収を評価します。
4. モデル仮定と感度を確認する
ロット内相関、過分散、時系列変化、欠測、停止ルール、事前分布を点検します。個体を独立とみなせない場合は、ロットを分析単位にするか階層モデル・クラスタ頑健標準誤差を検討します。
5. 推論と予測の検証を分ける
因果効果の推定には比較設計、予測には時間的に後のデータや別設備での外部検証が必要です。予測モデルは識別性能だけでなく校正、閾値別費用、データ変化も監視します。
6. 意思決定と再評価のルールを運用する
「採用・追加試験・見送り」の条件、承認者、再評価時期、量産後の停止基準を定めます。データ、コード、ライブラリ版、分析計画、承認記録を保存し、同じ結論を再現できるようにします。
まとめ
- 仮説検定は、帰無仮説のもとで観測差がどれほど極端かを評価する
- 尤度比・Wald・Score検定は評価位置が異なり、大標本では近い結果になりやすい
- p値は仮説が正しい確率ではなく、効果の大きさも表さない
- 信頼区間は、改善幅としてデータと整合する範囲を示す
- ベイズ推定は事前情報をデータで更新し、改善確率を直接計算できる
- ベイズ予測分布は、将来不良数の偶然変動と推定不確実性を統合する
- AIC・BICは適合度と複雑さを比較するが、目的や外部妥当性の判断を代替しない
- 統計的推論と機械学習は、効果説明と未知データ予測を役割分担できる
統計的推論の価値は、有意差を付けることではありません。限られた試験データを、誤りの可能性、改善幅、将来リスク、経済価値まで含む意思決定資料へ変換することにあります。
法人向けのご相談
数理工房では、製造業の品質・生産・保全データを対象に、以下のような支援を行っています。
- 工程変更・設備投資の統計的な試験計画と必要サンプル数設計
- p値、信頼区間、ベイズ推定を用いた意思決定ルールの設計
- 品質KPIと廃棄・手直し・保証損失を結ぶ効果金額化
- AIC・BIC、外部検証、校正を含む予測モデル評価
- Python notebookを用いた企業研修と分析プロトタイプ開発
- 分析結果を定期監視・承認フローへ組み込むシステム化
「有意差の読み方を社内で統一したい」「工程試験のサンプル数を設計したい」「PoCの精度を量産判断へつなげたい」といった段階からご相談いただけます。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。