100本ノック / 確率統計 / 確率・統計Python100本ノック
不確かさを更新して改善案を選ぶ：製造業のベイズ統計10本ノック
不確かさを更新して改善案を選ぶ：製造業のベイズ統計10本ノック
製造現場では、不適合率や停止頻度を限られた観測から推定し、追加検査、工程変更、改善案の展開を決めます。本記事では、架空の精密部品工場を題材に、事前分布と観測データを統合して不確かさを更新する方法を、ベータ・ベルヌーイ、ガンマ・ポアソン、MCMC、PyMC、階層ベイズ、予測分布、ベイズA/Bテストまで一続きで実装します。
ベイズ統計の価値は、単に「平均を推定する」ことではありません。「不適合率が管理上の基準を超える確率」「来週に何件起きるか」「改善案が現行より良い確率」のように、意思決定に直結する問いへ答えられることです。一方、結果は事前分布とモデル仮定に依存するため、感度分析と運用設計が欠かせません。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。 掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
架空の精密バルブ工場で、品質保証部と設備保全部が翌月の施策を決める場面を考えます。品質保証部は新しい洗浄条件でシール不適合が減ったかを判断し、設備保全部は突発停止に備える要員数を見積もります。また、複数工場の実績を比較し、データの少ない工場を過度に評価・非難しない仕組みも必要です。
ここで必要なのは、一つの推定値ではなく、未知量の確率分布です。観測が増えるたびに分布を更新し、その分布から超過確率、予測件数、施策の優越確率を計算します。
現場でよくある状況
- 不適合率を「不適合数 ÷ 検査数」の一点だけで報告し、母数の違いを無視して比較する
- 過去実績や類似工程の知見はあるが、担当者の勘としてしか利用されていない
- データが追加されても月次資料を作り直すだけで、更新規則が標準化されていない
- シミュレーションの収束や有効標本数を確認せず、MCMCの平均値だけを採用する
- 「改善後の方が良さそう」と「投資基準を満たすほど改善した」を区別していない
なぜこの問題は判断が難しいのか
不適合が0件でも、検査が10個なら「真の不適合率が0%」とは言えません。逆に、過去知見を強く置きすぎると、工程変化をデータが示しても更新が遅れます。ベイズの定理は、未知パラメータ に対して
と表されます。 は事前分布、 は尤度、 は事後分布です。計算式だけでなく、事前分布が誰のどの情報を表すか、データ生成過程が妥当か、どの確率を意思決定基準にするかを明文化する必要があります。
今回扱うノックの全体像
| No. | テーマ | 製造業での問い |
|---|---|---|
| 071 | ベータ・ベルヌーイ | シール不適合率をどの範囲とみるか |
| 072 | ガンマ・ポアソン | 1日あたりの設備停止率はどの程度か |
| 073 | ベイズ更新 | ロット追加のたびに判断がどう変わるか |
| 074 | MCMC | 解析解を使わず事後分布を近似できるか |
| 075 | Metropolis-Hastings | 提案幅と受理率は推定品質にどう影響するか |
| 076 | Gibbs Sampling | 平均とばらつきを交互に推定できるか |
| 077 | PyMC入門 | 宣言的なモデルで不適合率を推定できるか |
| 078 | 階層ベイズ | 工場別の少数データをどう安定化するか |
| 079 | ベイズ予測分布 | 次回検査で基準超過する確率は何%か |
| 080 | ベイズA/Bテスト | 新条件を展開する根拠は十分か |
Python 環境の準備
NumPyで乱数・数値計算、pandasで表、SciPyで確率分布、matplotlibで可視化し、No.077ではPyMCを使います。日本語表示には japanize_matplotlib を使用します。外部データとseabornは使いません。乱数seedを固定し、同じ環境なら同じ結果を再現できるようにします。
import platform
import warnings
import numpy as np
import pandas as pd
import scipy
from scipy import stats
from scipy.special import betaln, expit
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib
import pytensor
pytensor.config.cxx = "" # C++ツールチェーンがない環境でも実行できるPython実装を使用
import pymc as pm
from IPython.display import display
warnings.filterwarnings("ignore", category=FutureWarning)
pd.set_option("display.precision", 4)
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.unicode_minus"] = False
SEED = 20260711
rng = np.random.default_rng(SEED)
pd.DataFrame({
"項目": ["Python", "NumPy", "pandas", "SciPy", "matplotlib", "PyMC", "乱数seed"],
"値": [platform.python_version(), np.__version__, pd.__version__, scipy.__version__,
matplotlib.__version__, pm.__version__, SEED],
})| 項目 | 値 | |
|---|---|---|
| 0 | Python | 3.13.1 |
| 1 | NumPy | 2.5.1 |
| 2 | pandas | 3.0.3 |
| 3 | SciPy | 1.18.0 |
| 4 | matplotlib | 3.11.0 |
| 5 | PyMC | 5.26.1 |
| 6 | 乱数seed | 20260711 |
架空データの作成
同じ工場内の複数の意思決定を想定し、(1) 12ロットのシール検査、(2) 35日間の設備停止件数、(3) 充填時間、(4) 6工場の検査実績、(5) 現行条件と新洗浄条件の比較データを生成します。実務では、工程・設備・材料・測定方法が期間中に変わっていないかを確認し、変化がある場合は層別または時変モデルを検討します。
# シール検査：ロットごとに検査数が異なる
lot_n = rng.integers(45, 76, size=12)
lot_defects = rng.binomial(lot_n, 0.038)
inspection_df = pd.DataFrame({
"ロット": [f"L{i:02d}" for i in range(1, 13)],
"検査数": lot_n,
"不適合数": lot_defects,
})
inspection_df["不適合率"] = inspection_df["不適合数"] / inspection_df["検査数"]
# 設備停止、充填時間、工場別検査、A/B比較
daily_stops = rng.poisson(1.35, size=35)
fill_time = rng.normal(42.4, 1.8, size=45)
plant_names = ["東北", "関東", "中部", "関西", "中国", "九州"]
plant_n = np.array([80, 520, 140, 950, 65, 300])
plant_true_p = np.array([0.042, 0.031, 0.050, 0.027, 0.060, 0.036])
plant_k = rng.binomial(plant_n, plant_true_p)
ab_n = {"現行条件": 650, "新洗浄条件": 620}
ab_k = {
"現行条件": rng.binomial(ab_n["現行条件"], 0.046),
"新洗浄条件": rng.binomial(ab_n["新洗浄条件"], 0.028),
}
display(inspection_df)
display(pd.DataFrame({
"データ": ["設備停止", "充填時間", "工場別検査", "A/B比較"],
"観測規模": [f"{len(daily_stops)}日", f"{len(fill_time)}個", f"{len(plant_names)}工場",
f"{sum(ab_n.values())}個"],
}))| ロット | 検査数 | 不適合数 | 不適合率 | |
|---|---|---|---|---|
| 0 | L01 | 71 | 1 | 0.0141 |
| 1 | L02 | 51 | 4 | 0.0784 |
| 2 | L03 | 50 | 2 | 0.0400 |
| 3 | L04 | 73 | 2 | 0.0274 |
| 4 | L05 | 47 | 3 | 0.0638 |
| 5 | L06 | 61 | 2 | 0.0328 |
| 6 | L07 | 64 | 3 | 0.0469 |
| 7 | L08 | 65 | 3 | 0.0462 |
| 8 | L09 | 72 | 3 | 0.0417 |
| 9 | L10 | 71 | 1 | 0.0141 |
| 10 | L11 | 55 | 2 | 0.0364 |
| 11 | L12 | 51 | 3 | 0.0588 |
| データ | 観測規模 | |
|---|---|---|
| 0 | 設備停止 | 35日 |
| 1 | 充填時間 | 45個 |
| 2 | 工場別検査 | 6工場 |
| 3 | A/B比較 | 1270個 |
No.071：ベータ・ベルヌーイ — 不適合率を分布で推定する
実務での意味
個々の製品を適合0・不適合1と表すと、未知の不適合率 をベルヌーイ確率として扱えます。点推定だけでなく、 の信用区間や「管理基準3%を超える確率」を求めることで、追加検査や原因調査の優先度を決められます。
分析・モデル化の考え方
事前分布を 、観測した不適合数を 、検査数を とすると、共役性により
です。ここでは弱い事前分布 を置きます。実務で過去実績を用いる場合は、その根拠となる期間・設備・製品族と、事前分布の実効標本数 を併記します。
Pythonで確認する
n_total = int(inspection_df["検査数"].sum())
k_total = int(inspection_df["不適合数"].sum())
a0, b0 = 1.0, 1.0
a_post, b_post = a0 + k_total, b0 + n_total - k_total
ci71 = stats.beta.ppf([0.025, 0.975], a_post, b_post)
threshold = 0.03
beta_summary = pd.DataFrame({
"指標": ["検査数", "不適合数", "観測率", "事後平均", "95%信用区間下限", "95%信用区間上限", "P(p > 3%)"],
"値": [n_total, k_total, k_total/n_total, a_post/(a_post+b_post), *ci71,
stats.beta.sf(threshold, a_post, b_post)],
})
display(beta_summary)
x = np.linspace(0, 0.09, 500)
plt.plot(x, stats.beta.pdf(x, a0, b0), label="事前分布 Beta(1, 1)")
plt.plot(x, stats.beta.pdf(x, a_post, b_post), label="事後分布", linewidth=2)
plt.axvline(threshold, color="crimson", linestyle="--", label="管理基準 3%")
plt.title("シール不適合率の事前分布と事後分布")
plt.xlabel("不適合率 p")
plt.ylabel("確率密度")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 指標 | 値 | |
|---|---|---|
| 0 | 検査数 | 731.0000 |
| 1 | 不適合数 | 29.0000 |
| 2 | 観測率 | 0.0397 |
| 3 | 事後平均 | 0.0409 |
| 4 | 95%信用区間下限 | 0.0278 |
| 5 | 95%信用区間上限 | 0.0564 |
| 6 | P(p > 3%) | 0.9437 |
結果の読み取り
事後平均は今回の不適合率を平滑化した値で、95%信用区間は「モデルと事前分布の下で、未知の が95%の確率で入る範囲」です。管理基準3%を超える事後確率が高ければ、単純な合否判定だけでなく、追加検査や工程要因の層別を優先する根拠になります。ただし、この確率はロット間で同じ が続くという仮定に依存します。
No.072：ガンマ・ポアソン — 設備停止の発生率を更新する
実務での意味
日ごとの突発停止件数をモデル化すると、保全当番、予備品、復旧支援の必要量を確率で見積もれます。停止の有無だけでなく件数を扱うため、一定期間あたりの発生率 を推定します。
分析・モデル化の考え方
、率パラメータ表記の事前分布を とすると、 日の総件数 を観測した事後分布は
です。SciPyの scale は率の逆数 で指定する点に注意します。
Pythonで確認する
ga, gb = 2.0, 1.5 # 事前平均 1.33件/日
ga_post = ga + daily_stops.sum()
gb_post = gb + len(daily_stops)
lambda_mean = ga_post / gb_post
lambda_ci = stats.gamma.ppf([0.025, 0.975], ga_post, scale=1/gb_post)
display(pd.DataFrame({
"指標": ["観測日数", "総停止件数", "標本平均", "事後平均", "95%信用区間下限", "95%信用区間上限"],
"値": [len(daily_stops), daily_stops.sum(), daily_stops.mean(), lambda_mean, *lambda_ci],
}))
lam_x = np.linspace(0.3, 2.8, 500)
plt.plot(lam_x, stats.gamma.pdf(lam_x, ga, scale=1/gb), label="事前分布")
plt.plot(lam_x, stats.gamma.pdf(lam_x, ga_post, scale=1/gb_post), label="事後分布", linewidth=2)
plt.title("1日あたり設備停止率のベイズ更新")
plt.xlabel("停止発生率 λ（件/日）")
plt.ylabel("確率密度")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 指標 | 値 | |
|---|---|---|
| 0 | 観測日数 | 35.0000 |
| 1 | 総停止件数 | 51.0000 |
| 2 | 標本平均 | 1.4571 |
| 3 | 事後平均 | 1.4521 |
| 4 | 95%信用区間下限 | 1.0877 |
| 5 | 95%信用区間上限 | 1.8682 |
結果の読み取り
標本平均と事後平均が近いのは、35日分の観測が弱い事前情報より強く効いているためです。信用区間は、保全計画を一点で固定することの危険を示します。ただしポアソン分布は、日ごとの発生率が一定で平均と分散が等しいと仮定します。曜日、品種、稼働時間、故障の連鎖で過分散が生じる場合は、負の二項モデルや共変量を含むモデルへ拡張します。
No.073：ベイズ更新 — ロット到着ごとに判断を更新する
実務での意味
月末に全データをまとめるだけでなく、ロット検査が終わるたびに不適合率の分布を更新すれば、追加検査や工程停止の判断を早められます。更新履歴を残すと、いつ、どのデータで判断が変わったかを説明できます。
分析・モデル化の考え方
ベータ分布のパラメータは「不適合側の擬似件数」と「適合側の擬似件数」と解釈できます。前ロットまでの事後分布を次ロットの事前分布にする逐次更新は、全ロットを一括更新する計算と一致します。ここでは事後平均、95%信用区間、 の推移を記録します。
Pythonで確認する
a_seq, b_seq = a0, b0
update_rows = []
for row in inspection_df.itertuples(index=False):
a_seq += row.不適合数
b_seq += row.検査数 - row.不適合数
lo, hi = stats.beta.ppf([0.025, 0.975], a_seq, b_seq)
update_rows.append({
"ロット": row.ロット, "累積検査数": int(a_seq + b_seq - a0 - b0),
"事後平均": a_seq/(a_seq+b_seq), "下限": lo, "上限": hi,
"P(p>3%)": stats.beta.sf(0.03, a_seq, b_seq),
})
update_df = pd.DataFrame(update_rows)
display(update_df)
print("一括更新と逐次更新のパラメータ一致:", (a_seq, b_seq) == (a_post, b_post))
idx = np.arange(1, len(update_df)+1)
plt.plot(idx, update_df["事後平均"], marker="o", label="事後平均")
plt.fill_between(idx, update_df["下限"], update_df["上限"], alpha=0.22, label="95%信用区間")
plt.axhline(0.03, color="crimson", linestyle="--", label="管理基準 3%")
plt.title("ロット追加に伴う不適合率推定の更新")
plt.xlabel("更新したロット数")
plt.ylabel("不適合率")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| ロット | 累積検査数 | 事後平均 | 下限 | 上限 | P(p>3%) | |
|---|---|---|---|---|---|---|
| 0 | L01 | 71 | 0.0274 | 0.0034 | 0.0750 | 0.3600 |
| 1 | L02 | 122 | 0.0484 | 0.0181 | 0.0923 | 0.8344 |
| 2 | L03 | 172 | 0.0460 | 0.0202 | 0.0816 | 0.8492 |
| 3 | L04 | 245 | 0.0405 | 0.0197 | 0.0683 | 0.7928 |
| 4 | L05 | 292 | 0.0442 | 0.0238 | 0.0704 | 0.8938 |
| 5 | L06 | 353 | 0.0423 | 0.0239 | 0.0655 | 0.8836 |
| 6 | L07 | 417 | 0.0430 | 0.0257 | 0.0643 | 0.9172 |
| 7 | L08 | 482 | 0.0434 | 0.0271 | 0.0632 | 0.9395 |
| 8 | L09 | 554 | 0.0432 | 0.0279 | 0.0615 | 0.9499 |
| 9 | L10 | 625 | 0.0399 | 0.0260 | 0.0565 | 0.9060 |
| 10 | L11 | 680 | 0.0396 | 0.0263 | 0.0554 | 0.9097 |
| 11 | L12 | 731 | 0.0409 | 0.0278 | 0.0564 | 0.9437 |
一括更新と逐次更新のパラメータ一致: True
結果の読み取り
初期は信用区間が広く、少数の不適合で推定が大きく動きます。検査数が蓄積すると区間が狭まり、後半の1ロットが与える影響は小さくなります。一括更新との一致は実装検算になります。実務の停止ルールには、単一時点の超過だけでなく、確率が何回連続で閾値を超えたか、停止損失と流出損失のどちらが大きいかも組み込みます。
No.074：MCMC — 事後分布を標本で近似する
実務での意味
共役分布が使えない複雑なモデルでは、事後分布の平均や区間を式で計算できません。Markov Chain Monte Carlo（MCMC）は、事後分布を定常分布に持つ連鎖から標本を作り、設備停止率などの不確かさを近似します。
分析・モデル化の考え方
No.072のガンマ事後分布をあえて解析解を使わず、 上のランダムウォークで近似します。変数変換のヤコビアンを含む対数密度は定数を除き
です。バーンインを除き、トレース、自己相関、Monte Carlo標準誤差を確認します。
Pythonで確認する
def log_target_z(z):
return ga_post * z - gb_post * np.exp(z)
mcmc_rng = np.random.default_rng(SEED + 74)
n_iter, burn = 16000, 3000
z = np.log(lambda_mean)
z_trace = np.empty(n_iter)
accepted = 0
for i in range(n_iter):
proposal = z + mcmc_rng.normal(0, 0.22)
if np.log(mcmc_rng.random()) < log_target_z(proposal) - log_target_z(z):
z = proposal
accepted += 1
z_trace[i] = z
lambda_trace = np.exp(z_trace[burn:])
def autocorr(x, max_lag=40):
centered = x - x.mean()
ac = np.correlate(centered, centered, mode="full")[len(x)-1:len(x)+max_lag]
return ac / ac[0]
acf = autocorr(lambda_trace)
positive_acf = acf[1:][acf[1:] > 0]
ess_approx = len(lambda_trace) / (1 + 2 * positive_acf.sum())
mcmc_summary = pd.DataFrame({
"指標": ["受理率", "MCMC平均", "解析的な事後平均", "MCMC 2.5%", "MCMC 97.5%", "概算ESS", "MC標準誤差"],
"値": [accepted/n_iter, lambda_trace.mean(), lambda_mean,
*np.quantile(lambda_trace, [0.025, 0.975]), ess_approx,
lambda_trace.std(ddof=1)/np.sqrt(ess_approx)],
})
display(mcmc_summary)
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(lambda_trace[:2500], linewidth=0.6)
axes[0].set_title("MCMCトレース（先頭2,500標本）")
axes[0].set_xlabel("反復")
axes[0].set_ylabel("停止発生率 λ")
axes[0].grid(alpha=0.3)
axes[1].bar(np.arange(1, 21), acf[1:21])
axes[1].set_title("標本の自己相関")
axes[1].set_xlabel("ラグ")
axes[1].set_ylabel("自己相関")
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()| 指標 | 値 | |
|---|---|---|
| 0 | 受理率 | 0.5665 |
| 1 | MCMC平均 | 1.4504 |
| 2 | 解析的な事後平均 | 1.4521 |
| 3 | MCMC 2.5% | 1.0835 |
| 4 | MCMC 97.5% | 1.8838 |
| 5 | 概算ESS | 2044.2499 |
| 6 | MC標準誤差 | 0.0045 |
結果の読み取り
MCMC平均と解析的な事後平均が近いことは、この例での実装検算です。標本数が多くても、連続する標本に自己相関があれば独立な情報量は少なく、ESSは総標本数を下回ります。実務では複数チェーン、、ESS、トレース、発散、初期値感度を確認し、Monte Carlo誤差が業務上必要な精度より十分小さいことを確認してから確率を報告します。
No.075：Metropolis-Hastings — 提案幅と受理率を診断する
実務での意味
Metropolis-Hastings（MH）法は、現在値から候補を提案し、事後確率の比で採否を決めるMCMCの基本手法です。提案幅が小さすぎるとほぼ受理されても探索が遅く、大きすぎると棄却が増えて連鎖が動きません。
分析・モデル化の考え方
不適合率を必ず0〜1に保つため、 上で対称な正規提案を使います。ベータ事後分布とヤコビアンを合わせた対数標的密度は
です。複数の提案幅を比べ、受理率だけでなくESSと解析解からの誤差を確認します。
Pythonで確認する
def mh_beta(proposal_sd, seed, n_iter=12000, burn=2000):
local_rng = np.random.default_rng(seed)
z = np.log((a_post/(a_post+b_post)) / (1-a_post/(a_post+b_post)))
out = np.empty(n_iter)
accepted = 0
def log_target(logit_p):
p = expit(logit_p)
return a_post*np.log(p) + b_post*np.log1p(-p)
for i in range(n_iter):
prop = z + local_rng.normal(0, proposal_sd)
if np.log(local_rng.random()) < log_target(prop) - log_target(z):
z = prop
accepted += 1
out[i] = expit(z)
sample = out[burn:]
ac = autocorr(sample, 100)
pos = ac[1:][ac[1:] > 0]
ess = len(sample)/(1+2*pos.sum())
return sample, accepted/n_iter, ess
mh_rows = []
mh_samples = {}
for j, proposal_sd in enumerate([0.05, 0.35, 2.0]):
sample, acc, ess = mh_beta(proposal_sd, SEED + 750 + j)
mh_samples[proposal_sd] = sample
mh_rows.append({"提案幅": proposal_sd, "受理率": acc, "事後平均": sample.mean(),
"解析解との差": sample.mean()-a_post/(a_post+b_post), "概算ESS": ess})
mh_df = pd.DataFrame(mh_rows)
display(mh_df)
for proposal_sd, sample in mh_samples.items():
plt.plot(sample[:1000], linewidth=0.7, label=f"提案幅={proposal_sd}")
plt.title("提案幅によるMetropolis-Hastings連鎖の違い")
plt.xlabel("反復")
plt.ylabel("不適合率 p")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 提案幅 | 受理率 | 事後平均 | 解析解との差 | 概算ESS | |
|---|---|---|---|---|---|
| 0 | 0.05 | 0.9097 | 0.0406 | -0.0003 | 138.5534 |
| 1 | 0.35 | 0.5183 | 0.0412 | 0.0003 | 1814.6552 |
| 2 | 2.00 | 0.1148 | 0.0408 | -0.0001 | 845.5656 |
結果の読み取り
小さい提案幅は受理率が高くても少しずつしか移動せず、ESSが伸びません。大きい提案幅は同じ値に留まる区間が増えます。中間の幅がこの例では効率的ですが、万能な受理率はなく、次元や標的分布で変わります。受理率だけをKPIにせず、ESS/秒、複数チェーンの一致、推定対象のMonte Carlo誤差を併せて管理します。
No.076：Gibbs Sampling — 充填時間の平均とばらつきを交互更新する
実務での意味
充填時間の工程中心 と精度 はどちらも未知です。一方を固定すれば他方の条件付き分布から直接サンプリングできる場合、Gibbs Samplingで交互に更新できます。
分析・モデル化の考え方
に対し、共役なNormal-Gamma事前分布
を置きます。 と を交互に引きます。時間順のドリフトがない独立な正規データという仮定は、実務適用前に管理図や残差で確認します。
Pythonで確認する
gibbs_rng = np.random.default_rng(SEED + 76)
mu0, kappa0, shape0, rate0 = 42.0, 0.5, 2.0, 4.0
n = len(fill_time)
n_gibbs, gibbs_burn = 14000, 2000
mu, tau = fill_time.mean(), 1/fill_time.var()
mu_chain = np.empty(n_gibbs)
sigma_chain = np.empty(n_gibbs)
for i in range(n_gibbs):
# tau | mu, x
shape = shape0 + (n + 1)/2
rate = rate0 + 0.5*np.sum((fill_time-mu)**2) + 0.5*kappa0*(mu-mu0)**2
tau = gibbs_rng.gamma(shape, 1/rate)
# mu | tau, x
kappa_n = kappa0 + n
mu_n = (kappa0*mu0 + n*fill_time.mean())/kappa_n
mu = gibbs_rng.normal(mu_n, np.sqrt(1/(kappa_n*tau)))
mu_chain[i] = mu
sigma_chain[i] = 1/np.sqrt(tau)
mu_sample = mu_chain[gibbs_burn:]
sigma_sample = sigma_chain[gibbs_burn:]
display(pd.DataFrame({
"推定対象": ["平均充填時間 μ", "標準偏差 σ"],
"事後平均": [mu_sample.mean(), sigma_sample.mean()],
"95%下限": [np.quantile(mu_sample, .025), np.quantile(sigma_sample, .025)],
"95%上限": [np.quantile(mu_sample, .975), np.quantile(sigma_sample, .975)],
}))
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(mu_sample[:2000], linewidth=0.6)
axes[0].set_title("平均充填時間 μ のトレース")
axes[0].set_xlabel("反復")
axes[0].set_ylabel("秒")
axes[0].grid(alpha=0.3)
axes[1].scatter(mu_sample[::10], sigma_sample[::10], s=7, alpha=0.25)
axes[1].set_title("μ と σ の同時事後標本")
axes[1].set_xlabel("平均 μ（秒）")
axes[1].set_ylabel("標準偏差 σ（秒）")
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()| 推定対象 | 事後平均 | 95%下限 | 95%上限 | |
|---|---|---|---|---|
| 0 | 平均充填時間 μ | 42.6703 | 42.0901 | 43.2449 |
| 1 | 標準偏差 σ | 1.9306 | 1.5955 | 2.3738 |
結果の読み取り
平均と標準偏差を固定値ではなく同時分布として保持できるため、「平均が上がった可能性」と「ばらつきが増えた可能性」を分けて評価できます。Gibbs法でも標本は独立ではないため診断が必要です。また、正規モデルの標準偏差は短期変動と長期ドリフトを自動では分離しません。時刻、設備、品種を記録し、必要なら階層・時系列モデルへ拡張します。
No.077：PyMC入門 — モデルを宣言して不適合率を推定する
実務での意味
PyMCでは、確率変数と観測の関係をモデルとして宣言し、NUTSなどのMCMCで事後分布を計算できます。共役でない回帰、階層構造、測定誤差を含むモデルへ段階的に拡張しやすいことが利点です。
分析・モデル化の考え方
まずNo.071と同じ 、 をPyMCで表します。既知の解析解と照合できる単純なモデルから始めることで、ライブラリの設定、乱数、要約方法を検証します。target_accept は提案の受理を調整する値で、モデル妥当性の指標ではありません。
Pythonで確認する
with pm.Model() as defect_model:
p = pm.Beta("p", alpha=a0, beta=b0)
observed = pm.Binomial("observed", n=n_total, p=p, observed=k_total)
idata = pm.sample(
draws=1000, tune=1000, chains=2, cores=1,
random_seed=SEED + 77,
target_accept=0.9, progressbar=False,
compute_convergence_checks=True,
)
pymc_sample = idata.posterior["p"].values.reshape(-1)
pymc_summary = pm.stats.summary(idata, var_names=["p"], kind="all", round_to=5)
display(pymc_summary)
print(f"解析的事後平均: {a_post/(a_post+b_post):.5f}")
print(f"PyMC事後平均 : {pymc_sample.mean():.5f}")
plt.hist(pymc_sample, bins=35, density=True, alpha=0.55, label="PyMC標本")
plt.plot(x, stats.beta.pdf(x, a_post, b_post), linewidth=2, label="解析的Beta事後分布")
plt.title("PyMC標本と解析的事後分布の照合")
plt.xlabel("不適合率 p")
plt.ylabel("確率密度")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [p]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 1 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| p | 0.041 | 0.0073 | 0.0282 | 0.0549 | 0.0003 | 0.0001 | 736.415 | 1017.243 | 1.0014 |
解析的事後平均: 0.04093
PyMC事後平均 : 0.04103
結果の読み取り
PyMCの事後平均と解析解がMonte Carlo誤差の範囲で近く、ヒストグラムも重なれば基本実装を検算できます。r_hat は1に近いか、ess_bulk と ess_tail は十分かも確認します。本番ではモデルコード、データ版、PyMC版、seed、チェーン数、診断値を保存し、警告や発散を無視して結果だけ配布しない運用が必要です。
No.078：階層ベイズ — 工場別不適合率を部分的にプーリングする
実務での意味
工場別の不適合率を生の比率で順位付けすると、検査数が少ない工場ほど0%や高率になりやすく、偶然を実力差と誤認します。階層ベイズでは、各工場の情報を残しつつ、全工場に共通する分布から情報を借ります。
分析・モデル化の考え方
工場 について
とします。 は全体水準、 は工場間の類似度です。ここでは のグリッドに事前分布を置き、ベータ二項周辺尤度で重みを計算するため、固定したハイパーパラメータだけを使う経験ベイズではなく、その不確かさも事後分布へ反映します。
Pythonで確認する
m_grid = np.linspace(0.008, 0.10, 90)
kappa_grid = np.geomspace(4, 180, 80)
M, KAPPA = np.meshgrid(m_grid, kappa_grid, indexing="ij")
ALPHA = M*KAPPA
BETA = (1-M)*KAPPA
log_weight = -KAPPA/60 # kappaに緩やかな指数事前分布
for kj, nj in zip(plant_k, plant_n):
log_weight += betaln(ALPHA+kj, BETA+nj-kj) - betaln(ALPHA, BETA)
log_weight -= log_weight.max()
weight = np.exp(log_weight)
weight /= weight.sum()
flat_rng = np.random.default_rng(SEED + 78)
grid_indices = flat_rng.choice(weight.size, size=12000, p=weight.ravel())
alpha_draw = ALPHA.ravel()[grid_indices]
beta_draw = BETA.ravel()[grid_indices]
hier_draws = np.column_stack([
flat_rng.beta(alpha_draw+kj, beta_draw+nj-kj)
for kj, nj in zip(plant_k, plant_n)
])
plant_df = pd.DataFrame({
"工場": plant_names,
"検査数": plant_n,
"不適合数": plant_k,
"生の不適合率": plant_k/plant_n,
"階層事後平均": hier_draws.mean(axis=0),
"95%下限": np.quantile(hier_draws, .025, axis=0),
"95%上限": np.quantile(hier_draws, .975, axis=0),
})
display(plant_df)
pos = np.arange(len(plant_names))
plt.scatter(pos, plant_df["生の不適合率"], marker="x", s=70, label="生の不適合率")
plt.errorbar(pos, plant_df["階層事後平均"],
yerr=[plant_df["階層事後平均"]-plant_df["95%下限"],
plant_df["95%上限"]-plant_df["階層事後平均"]],
fmt="o", capsize=4, label="階層事後平均・95%区間")
plt.xticks(pos, plant_names)
plt.title("工場別不適合率の部分的プーリング")
plt.xlabel("工場")
plt.ylabel("不適合率")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 工場 | 検査数 | 不適合数 | 生の不適合率 | 階層事後平均 | 95%下限 | 95%上限 | |
|---|---|---|---|---|---|---|---|
| 0 | 東北 | 80 | 3 | 0.0375 | 0.0425 | 0.0148 | 0.0828 |
| 1 | 関東 | 520 | 16 | 0.0308 | 0.0327 | 0.0198 | 0.0489 |
| 2 | 中部 | 140 | 9 | 0.0643 | 0.0593 | 0.0305 | 0.0982 |
| 3 | 関西 | 950 | 27 | 0.0284 | 0.0297 | 0.0199 | 0.0411 |
| 4 | 中国 | 65 | 3 | 0.0462 | 0.0479 | 0.0171 | 0.0934 |
| 5 | 九州 | 300 | 15 | 0.0500 | 0.0497 | 0.0298 | 0.0745 |
結果の読み取り
検査数の少ない工場ほど全体水準への縮約が強く、区間も広くなります。これは工場差を消す処理ではなく、少数データの極端値を不確かさ込みで扱う部分的プーリングです。縮約後の順位を人事評価へ直結させず、材料構成、製品難度、測定基準などの比較可能性を確認します。共通分布を置けないほど工程が異なる工場は、別グループに分ける必要があります。
No.079：ベイズ予測分布 — 次回検査の不適合数を予測する
実務での意味
品質保証が知りたいのは真の不適合率だけでなく、「次の200個で何個の不適合が出るか」「5%以上になるリスクは何%か」です。予測分布はパラメータ不確かさと将来の偶然変動を両方含めます。
分析・モデル化の考え方
事後分布 の下で、将来 個の不適合数 はベータ二項分布
に従います。 を事後平均に固定した二項分布より広くなり、推定の不確かさを計画へ反映できます。
Pythonで確認する
future_n = 200
future_k = np.arange(future_n+1)
pred_pmf = stats.betabinom.pmf(future_k, future_n, a_post, b_post)
pred_interval = stats.betabinom.ppf([0.025, 0.975], future_n, a_post, b_post).astype(int)
risk_5pct = stats.betabinom.sf(9, future_n, a_post, b_post) # 10個以上
display(pd.DataFrame({
"指標": ["次回検査数", "予測不適合数の平均", "95%予測区間下限", "95%予測区間上限", "P(10個以上 = 5%以上)"],
"値": [future_n, future_n*a_post/(a_post+b_post), *pred_interval, risk_5pct],
}))
show = future_k <= 25
plt.bar(future_k[show], pred_pmf[show], alpha=0.75)
plt.axvline(10, color="crimson", linestyle="--", label="5%に相当する10個")
plt.title("次回200個検査における不適合数の事後予測分布")
plt.xlabel("将来の不適合数")
plt.ylabel("予測確率")
plt.grid(axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 指標 | 値 | |
|---|---|---|
| 0 | 次回検査数 | 200.0000 |
| 1 | 予測不適合数の平均 | 8.1855 |
| 2 | 95%予測区間下限 | 3.0000 |
| 3 | 95%予測区間上限 | 15.0000 |
| 4 | P(10個以上 = 5%以上) | 0.3164 |
結果の読み取り
予測平均だけでなく95%予測区間を見ると、次回検査で起こり得る振れ幅を把握できます。5%以上となる確率は、追加検査要員や隔離スペースを準備する判断材料です。信用区間は未知の率 の範囲、予測区間は将来の不適合数の範囲であり、用途が異なります。予測の妥当性は、後日実績との較正を継続的に確認します。
No.080：ベイズA/Bテスト — 新洗浄条件の展開を判断する
実務での意味
現行条件Aと新洗浄条件Bを比較するとき、「差があるか」だけでなく、BがAより良い確率、改善幅が投資基準を満たす確率、誤って展開した場合の損失を評価します。これにより、統計的な差と経営上意味のある差を分けられます。
分析・モデル化の考え方
両条件の不適合率に独立な 事前分布を置くと、それぞれの事後分布は解析的に得られます。事後標本から差 を計算し、 と、最小実務差1ポイントに対する を求めます。実際のA/Bでは、品種、材料ロット、設備、期間を無作為化または調整し、交絡を避けます。
Pythonで確認する
ab_rng = np.random.default_rng(SEED + 80)
draws = 100000
p_a = ab_rng.beta(1+ab_k["現行条件"], 1+ab_n["現行条件"]-ab_k["現行条件"], draws)
p_b = ab_rng.beta(1+ab_k["新洗浄条件"], 1+ab_n["新洗浄条件"]-ab_k["新洗浄条件"], draws)
delta = p_a - p_b
ab_table = pd.DataFrame({
"条件": ["現行条件A", "新洗浄条件B"],
"検査数": [ab_n["現行条件"], ab_n["新洗浄条件"]],
"不適合数": [ab_k["現行条件"], ab_k["新洗浄条件"]],
"観測不適合率": [ab_k["現行条件"]/ab_n["現行条件"], ab_k["新洗浄条件"]/ab_n["新洗浄条件"]],
"事後平均": [p_a.mean(), p_b.mean()],
})
display(ab_table)
display(pd.DataFrame({
"判断指標": ["P(BがAより低い)", "P(改善幅が1ポイント超)", "改善幅の事後平均", "改善幅の95%信用区間"],
"値": [f"{np.mean(delta>0):.3%}", f"{np.mean(delta>0.01):.3%}",
f"{delta.mean():.3%}", f"{np.quantile(delta,.025):.3%} 〜 {np.quantile(delta,.975):.3%}"],
}))
plt.hist(delta*100, bins=60, density=True, alpha=0.75)
plt.axvline(0, color="black", linestyle="--", label="差なし")
plt.axvline(1, color="crimson", linestyle="--", label="最小実務差 1ポイント")
plt.title("新洗浄条件による不適合率改善幅の事後分布")
plt.xlabel("改善幅 p(A) - p(B)（パーセントポイント）")
plt.ylabel("確率密度")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 条件 | 検査数 | 不適合数 | 観測不適合率 | 事後平均 | |
|---|---|---|---|---|---|
| 0 | 現行条件A | 650 | 32 | 0.0492 | 0.0507 |
| 1 | 新洗浄条件B | 620 | 19 | 0.0306 | 0.0322 |
| 判断指標 | 値 | |
|---|---|---|
| 0 | P(BがAより低い) | 95.283% |
| 1 | P(改善幅が1ポイント超) | 77.666% |
| 2 | 改善幅の事後平均 | 1.847% |
| 3 | 改善幅の95%信用区間 | -0.316% 〜 4.062% |
結果の読み取り
「BがAより低い確率」と「1ポイントを超えて改善する確率」は異なります。前者が高くても後者が投資基準に届かなければ、全面展開ではなく追加試験が妥当な場合があります。最終判断では、流出損失、条件変更費、タクトや他品質特性への副作用を損失関数に含めます。確率の閾値を結果を見た後に変えず、試験前に意思決定ルールを合意しておくことが重要です。
対象ノックを通して見える実務上の示唆
- 点ではなく確率で報告する：事後平均、信用区間、基準超過確率をセットにすると、追加検査・停止・展開の議論が具体化します。
- 事前分布を監査可能にする：過去データ、類似工程、専門家判断のどれを何件分の重みで入れたかを記録します。
- 予測とパラメータ推定を区別する：将来件数には工程率の不確かさと将来の偶然変動の両方が入ります。
- 少数群を順位付けしすぎない：階層ベイズで工場間の情報を共有し、検査数に応じた縮約と区間を示します。
- MCMC診断を成果物に含める：受理率だけでなく、トレース、複数チェーン、、ESS、発散、MC誤差を保存します。
- 統計差を経済価値へ変換する：A/Bの優越確率に加え、最小実務差と損失関数で意思決定します。
実務導入する場合に必要なこと
- 意思決定の定義：誰が、どの事後確率・予測確率で、追加検査・停止・展開を行うか事前に定める
- データ生成過程の確認：設備、金型、材料、品種、シフト、検査員、時刻を記録し、交換可能性や独立性を点検する
- 事前分布のガバナンス：根拠、更新日、対象範囲、実効標本数、弱情報事前分布との感度差をレビューする
- モデル検証：事後予測チェック、時系列外検証、較正、既知の単純モデルとの照合を行う
- 計算品質の管理：環境、依存ライブラリ、seed、チェーン、診断値、警告、実行ログを保存する
- 業務損失との接続：見逃し、過剰停止、検査、設備投資、品質流出の損失を定量化し、確率から行動へ変換する
- 段階導入：過去データで再現検証し、シャドー運用、限定工程での試行、標準作業化、定期監視へ進める
ベイズモデルは現場判断を自動的に正しくする装置ではありません。データ収集、工程知識、モデル診断、意思決定責任を一体で設計することで初めて価値が生まれます。
まとめ
No.071〜No.080では、ベータ・ベルヌーイとガンマ・ポアソンの共役更新から始め、逐次更新、MCMC、Metropolis-Hastings、Gibbs Sampling、PyMC、階層ベイズ、事後予測、ベイズA/Bテストまで実装しました。
中心となる考え方は、過去知見と観測を透明な規則で統合し、未知量と将来結果の不確かさを、現場が判断できる確率へ翻訳することです。事前分布・尤度・モデル構造・計算診断・損失関数を明示し、結果がどの仮定に依存するかを説明できる状態で運用しましょう。
法人向けのご相談
数理工房では、製造業における品質データ分析、ベイズ統計モデル、階層モデル、異常・故障予測、実験設計、Python研修、PoCから運用定着までをご支援しています。「少数データを活用したい」「不適合率や停止件数を確率で判断したい」「工場間比較を公平にしたい」「モデルの診断や運用基準を整えたい」といった段階からご相談いただけます。
📩 お問い合わせ: surikobo.co.jp/contact まずはお気軽にご相談ください。