100本ノック / ベイズ統計 / データ分析のためのベイズ統計100本ノック

製造ラインの不良率を階層ベイズで安定推定|少数データを意思決定に活かす

工場・ライン間のばらつきを意思決定につなげる階層ベイズ分析

タイトル・概要

同じ製品を造る工場やラインでも、不良率・歩留まり・温度感度には差があります。本稿では、架空の精密部品メーカーを題材に、**全体傾向を共有しつつ各グループの違いも残す「部分プーリング」**を中心に、No.071〜No.080を一つの流れで実装します。

扱う問いは「観測不良率が高い少量生産ラインを、直ちに問題ラインと判断してよいか」です。階層構造を使い、データ量に応じて推定の強さを変え、改善・追加測定・横展開を判断します。

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

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

多拠点・多品種の品質管理では、母数の大きい主力ラインと試作・少量ラインを同じ物差しで順位付けすると、偶然の振れが「実力差」に見えます。本稿では、グループ共通の情報を借りながら各ラインを推定し、限られた改善資源をどこへ配分するかを考えます。

現場でよくある状況

  • 月次の不良率ランキングだけで拠点を評価している
  • 新設ラインは検査数が少なく、0%や10%など極端な値が出る
  • 品種・材料・設備条件が混在し、単純平均では比較できない
  • 「全社平均で置き換える」か「各現場を完全に別扱いする」かの二択になっている

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

観測差には、真の工程差と標本誤差が同時に含まれます。特に検査数が少ない群ほど標本誤差が大きく、点推定だけでは区別できません。さらに、現場差を無視すれば異常を薄め、現場を完全に分離すれば推定が不安定になります。

今回扱うノックの全体像

No.071〜080では、階層構造、3種類のプーリング、率の推定、少数群の安定化、ランダム切片・係数、可視化、導入判断を順に扱います。中心となるベータ・二項モデルは

y_j\mid p_j\sim\mathrm{Binomial}(n_j,p_j),\qquad p_j\sim\mathrm{Beta}(lpha,eta)

です。事後分布は p_j\mid y_j\sim\mathrm{Beta}(lpha+y_j,eta+n_j-y_j) となり、事後平均は

E[p_j\mid y_j]= rac{n_j}{n_j+lpha+eta} rac{y_j}{n_j}+ rac{lpha+eta}{n_j+lpha+eta} rac{lpha}{lpha+eta}

です。すなわち、データが多い群は実績を重く、少ない群は全体傾向を重く見ます。

Python 環境の準備

外部データは使わず、乱数 seed を固定します。scipy は分布と区間推定、matplotlib は可視化に用います。

import sys
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy import stats

SEED = 20260712
rng = np.random.default_rng(SEED)
pd.set_option("display.precision", 4)
print(f"Python {sys.version.split()[0]}")
print(f"NumPy {np.__version__} / pandas {pd.__version__} / matplotlib {matplotlib.__version__}")
Python 3.13.1
NumPy 2.5.1 / pandas 3.0.3 / matplotlib 3.11.0

架空データの作成

8製造ラインの検査数と不良数、6製品の合格数、材料セグメント別の改善反応、ライン別の日次工程データを生成します。差が見えるよう真値を設定していますが、分析には真値を使いません。

lines = [f"Line-{c}" for c in "ABCDEFGH"]
n_inspected = np.array([1200, 900, 650, 400, 180, 80, 35, 12])
true_defect = np.array([0.018, 0.022, 0.027, 0.020, 0.035, 0.025, 0.045, 0.030])
defects = rng.binomial(n_inspected, true_defect)
line_df = pd.DataFrame({"line": lines, "inspected": n_inspected, "defects": defects})
line_df["observed_rate"] = line_df["defects"] / line_df["inspected"]

products = [f"Product-{i}" for i in range(1, 7)]
prod_n = np.array([1000, 700, 350, 160, 60, 20])
true_yield = np.array([0.965, 0.952, 0.970, 0.940, 0.955, 0.930])
passed = rng.binomial(prod_n, true_yield)
product_df = pd.DataFrame({"product": products, "produced": prod_n, "passed": passed})
product_df["observed_yield"] = product_df["passed"] / product_df["produced"]

line_df
line inspected defects observed_rate
0 Line-A 1200 24 0.0200
1 Line-B 900 19 0.0211
2 Line-C 650 22 0.0338
3 Line-D 400 7 0.0175
4 Line-E 180 2 0.0111
5 Line-F 80 1 0.0125
6 Line-G 35 4 0.1143
7 Line-H 12 1 0.0833

No.071:階層ベイズの考え方を理解する

実務での意味

ラインは独立した現場である一方、同じ設備標準・材料規格・品質制度を共有しています。この「違うが無関係ではない」を表現するのが階層モデルです。

分析・モデル化の考え方

ここでは全社の不良率を平均2.5%、事前情報量を40件相当とする Beta(1,39)\mathrm{Beta}(1,39) を共通母集団とします。これは結論を固定する値ではなく、後で感度分析すべき仮定です。

Pythonで確認する

alpha0, beta0 = 1.0, 39.0
line_df["post_alpha"] = alpha0 + line_df["defects"]
line_df["post_beta"] = beta0 + line_df["inspected"] - line_df["defects"]
line_df["hier_rate"] = line_df["post_alpha"] / (line_df["post_alpha"] + line_df["post_beta"])
line_df[["line", "inspected", "defects", "observed_rate", "hier_rate"]]
line inspected defects observed_rate hier_rate
0 Line-A 1200 24 0.0200 0.0202
1 Line-B 900 19 0.0211 0.0213
2 Line-C 650 22 0.0338 0.0333
3 Line-D 400 7 0.0175 0.0182
4 Line-E 180 2 0.0111 0.0136
5 Line-F 80 1 0.0125 0.0167
6 Line-G 35 4 0.1143 0.0667
7 Line-H 12 1 0.0833 0.0385

結果の読み取り

検査数が多いラインでは階層推定値は実績率に近く、少数ラインでは2.5%へ強く寄ります。この動きはデータを隠す補正ではなく、標本誤差の大きさを推定へ反映したものです。

No.072:完全プーリング・非プーリング・部分プーリングを比較する

実務での意味

完全プーリングは全ラインを同一視し、非プーリングは各ラインを別会社のように扱います。部分プーリングは、その中間でデータ量に応じて情報共有します。

分析・モデル化の考え方

完全プーリングは全体不良率、非プーリングは観測率、部分プーリングは先ほどの事後平均です。3案を並べると、評価制度が暗黙に置いている仮定を説明できます。

Pythonで確認する

pooled_rate = line_df["defects"].sum() / line_df["inspected"].sum()
comparison = line_df[["line", "inspected", "observed_rate", "hier_rate"]].copy()
comparison["complete_pooling"] = pooled_rate
comparison = comparison.rename(columns={"observed_rate": "no_pooling", "hier_rate": "partial_pooling"})
comparison[["line", "inspected", "complete_pooling", "no_pooling", "partial_pooling"]]
line inspected complete_pooling no_pooling partial_pooling
0 Line-A 1200 0.0231 0.0200 0.0202
1 Line-B 900 0.0231 0.0211 0.0213
2 Line-C 650 0.0231 0.0338 0.0333
3 Line-D 400 0.0231 0.0175 0.0182
4 Line-E 180 0.0231 0.0111 0.0136
5 Line-F 80 0.0231 0.0125 0.0167
6 Line-G 35 0.0231 0.1143 0.0667
7 Line-H 12 0.0231 0.0833 0.0385
plot_df = comparison.set_index("line")[["complete_pooling", "no_pooling", "partial_pooling"]]
ax = plot_df.plot(kind="bar", figsize=(10, 4))
ax.set_title("Three pooling strategies for defect rates")
ax.set_xlabel("Production line")
ax.set_ylabel("Defect rate")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

png

結果の読み取り

完全プーリングでは改善対象の差が消え、非プーリングでは少数群の極端値が強調されます。部分プーリングは差を残しながら極端な順位変動を抑えるため、一次スクリーニングに適します。

No.073:店舗別CVRを階層ベイズで推定する

実務での意味

原題の店舗別CVRは、製造業では「検査機会に対する不良発生率」と同じ二項データです。拠点・ライン別の品質率へ読み替えられます。

分析・モデル化の考え方

各ラインの事後分布から95%信用区間と「不良率が3%を超える確率」を求めます。点推定ではなく、品質基準超過の確率で優先順位を付けます。

Pythonで確認する

line_df["ci_low"] = stats.beta.ppf(0.025, line_df["post_alpha"], line_df["post_beta"])
line_df["ci_high"] = stats.beta.ppf(0.975, line_df["post_alpha"], line_df["post_beta"])
line_df["prob_over_3pct"] = 1 - stats.beta.cdf(0.03, line_df["post_alpha"], line_df["post_beta"])
line_df[["line", "hier_rate", "ci_low", "ci_high", "prob_over_3pct"]].sort_values("prob_over_3pct", ascending=False)
line hier_rate ci_low ci_high prob_over_3pct
6 Line-G 0.0667 0.0223 0.1327 0.9283
2 Line-C 0.0333 0.0213 0.0479 0.6692
7 Line-H 0.0385 0.0048 0.1045 0.5452
5 Line-F 0.0167 0.0020 0.0459 0.1248
3 Line-D 0.0182 0.0079 0.0326 0.0470
1 Line-B 0.0213 0.0131 0.0314 0.0425
4 Line-E 0.0136 0.0028 0.0326 0.0388
0 Line-A 0.0202 0.0131 0.0287 0.0133

結果の読み取り

prob_over_3pct が高いラインは、単なる高い観測率ではなく「基準超過の確からしさ」が高い候補です。信用区間が広いラインは、断定より追加検査を優先します。

No.074:商品別売上率を階層ベイズで推定する

実務での意味

原題の売上率を製品別歩留まりに読み替えます。製品ごとの母数が異なる場面でも、新製品の少量実績を既存製品群の情報で安定化できます。

分析・モデル化の考え方

歩留まりの共通事前分布を平均95%、40件相当の Beta(38,2)\mathrm{Beta}(38,2) とし、製品別の事後分布を求めます。

Pythonで確認する

yield_a0, yield_b0 = 38.0, 2.0
product_df["post_a"] = yield_a0 + product_df["passed"]
product_df["post_b"] = yield_b0 + product_df["produced"] - product_df["passed"]
product_df["hier_yield"] = product_df["post_a"] / (product_df["post_a"] + product_df["post_b"])
product_df["prob_below_94pct"] = stats.beta.cdf(0.94, product_df["post_a"], product_df["post_b"])
product_df
product produced passed observed_yield post_a post_b hier_yield prob_below_94pct
0 Product-1 1000 965 0.9650 1003.0 37.0 0.9644 1.4545e-04
1 Product-2 700 664 0.9486 702.0 38.0 0.9486 1.4382e-01
2 Product-3 350 346 0.9886 384.0 6.0 0.9846 3.3491e-06
3 Product-4 160 145 0.9062 183.0 17.0 0.9150 9.0861e-01
4 Product-5 60 58 0.9667 96.0 4.0 0.9600 1.4837e-01
5 Product-6 20 18 0.9000 56.0 4.0 0.9333 5.2445e-01

結果の読み取り

少量製品の0件・数件の不合格をそのまま恒常的な製品差とみなさず、既存製品群の水準を加味できます。一方、基準未達確率が高い製品は工程条件や設計公差のレビュー対象です。

No.075:顧客セグメント別の反応率を推定する

実務での意味

顧客反応率は、製造現場では材料サプライヤー区分別の改善処置成功率へ置き換えられます。セグメント別の効果差を、試行数の差込みで比較します。

分析・モデル化の考え方

共通の改善成功率60%、20件相当の事前分布から更新し、成功率が60%を超える確率を算出します。

Pythonで確認する

segments = ["Domestic-A", "Domestic-B", "Overseas-A", "Overseas-B"]
trials = np.array([90, 45, 18, 7])
true_response = np.array([0.72, 0.64, 0.55, 0.70])
success = rng.binomial(trials, true_response)
segment_df = pd.DataFrame({"segment": segments, "trials": trials, "success": success})
segment_df["observed"] = segment_df["success"] / segment_df["trials"]
a_seg, b_seg = 12.0, 8.0
segment_df["posterior"] = (a_seg + segment_df["success"]) / (a_seg + b_seg + segment_df["trials"])
segment_df["prob_over_60pct"] = 1 - stats.beta.cdf(0.60, a_seg + segment_df["success"], b_seg + segment_df["trials"] - segment_df["success"])
segment_df
segment trials success observed posterior prob_over_60pct
0 Domestic-A 90 66 0.7333 0.7091 0.9919
1 Domestic-B 45 33 0.7333 0.6923 0.9421
2 Overseas-A 18 7 0.3889 0.5000 0.1080
3 Overseas-B 7 4 0.5714 0.5926 0.4787

結果の読み取り

件数の少ない区分で高い実績率が出ても、事後確率が十分高いとは限りません。全社展開は期待値だけでなく、最低限の確度や追加試験費用と組み合わせて決めます。

No.076:少数データグループの推定を安定化する

実務での意味

試作ラインや新製品はデータが少ないため、ランキングが毎月大きく変わります。縮約の強さを数値で示すと、現場へ「なぜ補正されたか」を説明できます。

分析・モデル化の考え方

事後平均に占める実績の重みは wj=nj/(nj+40)w_j=n_j/(n_j+40) です。残りの 1wj1-w_j が共通事前情報の重みです。

Pythonで確認する

line_df["data_weight"] = line_df["inspected"] / (line_df["inspected"] + alpha0 + beta0)
line_df["shrinkage"] = line_df["observed_rate"] - line_df["hier_rate"]
line_df[["line", "inspected", "observed_rate", "hier_rate", "data_weight", "shrinkage"]]
line inspected observed_rate hier_rate data_weight shrinkage
0 Line-A 1200 0.0200 0.0202 0.9677 -0.0002
1 Line-B 900 0.0211 0.0213 0.9574 -0.0002
2 Line-C 650 0.0338 0.0333 0.9420 0.0005
3 Line-D 400 0.0175 0.0182 0.9091 -0.0007
4 Line-E 180 0.0111 0.0136 0.8182 -0.0025
5 Line-F 80 0.0125 0.0167 0.6667 -0.0042
6 Line-G 35 0.1143 0.0667 0.4667 0.0476
7 Line-H 12 0.0833 0.0385 0.2308 0.0449
fig, ax = plt.subplots(figsize=(8, 4))
ax.scatter(line_df["inspected"], line_df["data_weight"], s=60)
for _, r in line_df.iterrows(): ax.annotate(r["line"], (r["inspected"], r["data_weight"]), xytext=(4, 4), textcoords="offset points")
ax.set_title("Data volume and weight assigned to observed results")
ax.set_xlabel("Number inspected")
ax.set_ylabel("Weight on observed rate")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

png

結果の読み取り

検査数12件の群は主に全体傾向から推定され、1,200件の群はほぼ自群データで決まります。少数群への判断は「正常」と断定せず、推定区間を狭めるための追加測定計画につなげます。

No.077:ランダム切片モデルを理解する

実務での意味

同じ温度・設備条件でも、ライン固有の基準歩留まりがあります。ランダム切片は、このベースライン差を全社分布から来た偏差として扱います。

分析・モデル化の考え方

日次歩留まり y_{ij}=eta_0+u_j+ arepsilon_{ij} とし、まずライン平均を求めます。既知の群間・群内分散を用いた正規—正規モデルの簡易縮約で uju_j を推定します。

Pythonで確認する

days = 30
intercepts = np.array([0.010, -0.004, 0.006, -0.012, 0.003, 0.015, -0.008, 0.000])
daily = []
for line, u in zip(lines, intercepts):
    for day in range(days):
        daily.append((line, day + 1, 0.95 + u + rng.normal(0, 0.012)))
daily_df = pd.DataFrame(daily, columns=["line", "day", "yield_rate"])
means = daily_df.groupby("line")["yield_rate"].agg(["mean", "count", "std"]).reset_index()
grand = daily_df["yield_rate"].mean()
tau2, sigma2 = 0.010**2, 0.012**2
means["reliability"] = tau2 / (tau2 + sigma2 / means["count"])
means["random_intercept_est"] = means["reliability"] * (means["mean"] - grand)
means
line mean count std reliability random_intercept_est
0 Line-A 0.9613 30 0.0111 0.9542 0.0101
1 Line-B 0.9437 30 0.0121 0.9542 -0.0066
2 Line-C 0.9581 30 0.0079 0.9542 0.0071
3 Line-D 0.9350 30 0.0131 0.9542 -0.0150
4 Line-E 0.9514 30 0.0126 0.9542 0.0006
5 Line-F 0.9620 30 0.0145 0.9542 0.0108
6 Line-G 0.9404 30 0.0106 0.9542 -0.0098
7 Line-H 0.9536 30 0.0149 0.9542 0.0027

結果の読み取り

推定されたランダム切片が正なら、条件調整後も平均より高い基準歩留まりを持つ可能性があります。恒久的な現場差と断定する前に、設備世代・班・材料ロットなど未投入要因を確認します。

No.078:ランダム係数モデルを理解する

実務での意味

温度上昇の影響はラインごとに異なることがあります。共通係数だけでは局所的な弱点を見逃し、ライン別回帰だけではノイズを拾います。

分析・モデル化の考え方

y_{ij}=eta_0+u_{0j}+(eta_1+u_{1j})x_{ij}+ arepsilon_{ij} と考え、ライン別の温度係数を計算後、全体係数へ部分縮約します。ここでは仕組みを明示する簡易推定です。

Pythonで確認する

slope_true = np.array([-0.0010, -0.0018, -0.0007, -0.0022, -0.0013, -0.0005, -0.0020, -0.0011])
rows = []
for line, b0, b1 in zip(lines, intercepts, slope_true):
    temp = rng.normal(0, 2.5, days)
    y = 0.95 + b0 + b1 * temp + rng.normal(0, 0.008, days)
    rows.extend(zip([line]*days, temp, y))
temp_df = pd.DataFrame(rows, columns=["line", "temp_deviation", "yield_rate"])
est = []
for line, g in temp_df.groupby("line"):
    slope, intercept = np.polyfit(g["temp_deviation"], g["yield_rate"], 1)
    est.append((line, intercept, slope))
slope_df = pd.DataFrame(est, columns=["line", "intercept", "raw_slope"])
global_slope = slope_df["raw_slope"].mean()
slope_df["partial_slope"] = 0.70 * slope_df["raw_slope"] + 0.30 * global_slope
slope_df.sort_values("partial_slope")
line intercept raw_slope partial_slope
4 Line-E 0.9519 -0.0022 -0.0020
3 Line-D 0.9408 -0.0021 -0.0019
6 Line-G 0.9419 -0.0017 -0.0016
1 Line-B 0.9458 -0.0016 -0.0016
0 Line-A 0.9600 -0.0013 -0.0013
2 Line-C 0.9554 -0.0011 -0.0012
7 Line-H 0.9479 -0.0009 -0.0010
5 Line-F 0.9628 -0.0007 -0.0009

結果の読み取り

係数がより負のラインほど温度偏差に敏感です。保全・空調投資は単なる平均歩留まりだけでなく、環境変動への感度も含めて優先順位付けできます。

No.079:階層ベイズモデルの推定結果を可視化する

実務での意味

経営会議や品質会議では、平均値だけでなく不確実性と基準値の関係を一枚で伝える必要があります。

分析・モデル化の考え方

ライン別の事後平均と95%信用区間をフォレストプロットにし、管理上の注視水準3%を重ねます。

Pythonで確認する

viz = line_df.sort_values("hier_rate").reset_index(drop=True)
xerr = np.vstack([viz["hier_rate"] - viz["ci_low"], viz["ci_high"] - viz["hier_rate"]])
fig, ax = plt.subplots(figsize=(9, 5))
ax.errorbar(viz["hier_rate"], viz["line"], xerr=xerr, fmt="o", capsize=4)
ax.axvline(0.03, color="crimson", linestyle="--", label="Attention threshold: 3%")
ax.set_title("Posterior defect rates with 95% credible intervals")
ax.set_xlabel("Defect rate")
ax.set_ylabel("Production line")
ax.grid(axis="x", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

png

結果の読み取り

区間全体が基準より上なら改善介入の根拠が強く、区間が広く基準をまたぐ場合は追加検査が合理的です。「点が高い順」ではなく、位置と幅を同時に読みます。

No.080:階層ベイズを実務データ分析に使う判断基準を整理する

実務での意味

階層ベイズは常に必要ではありません。グループ数、群ごとのデータ量、共通性、誤判断コストを確認して採用を決めます。

分析・モデル化の考え方

候補案件を「階層構造」「少数群」「群間差」「誤判断コスト」で採点します。これは統計的証明ではなく、モデル化検討を始める業務チェックです。

Pythonで確認する

use_cases = pd.DataFrame({
    "case": ["Multi-line defect rate", "Single machine stable process", "New product yield", "Supplier improvement test"],
    "hierarchy": [2, 0, 2, 2], "small_groups": [2, 0, 2, 1],
    "between_group_diff": [2, 0, 1, 2], "decision_cost": [2, 1, 2, 2]
})
use_cases["score"] = use_cases.iloc[:, 1:].sum(axis=1)
use_cases["recommendation"] = pd.cut(use_cases["score"], [-1, 2, 5, 8], labels=["Simple model first", "Compare both", "Hierarchical model candidate"])
use_cases.sort_values("score", ascending=False)
case hierarchy small_groups between_group_diff decision_cost score recommendation
0 Multi-line defect rate 2 2 2 2 8 Hierarchical model candidate
2 New product yield 2 2 1 2 7 Hierarchical model candidate
3 Supplier improvement test 2 1 2 2 7 Hierarchical model candidate
1 Single machine stable process 0 0 0 1 1 Simple model first

結果の読み取り

多ライン品質や新製品歩留まりは候補になりやすく、単一設備の十分な時系列には単純モデルが先です。導入時は、単純モデルとの予測性能・解釈性・運用費を比較します。

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

  1. 順位ではなく確率で判断する:基準超過確率と信用区間を改善会議の共通指標にする。
  2. 少数群は「良い・悪い」より「情報不足」:追加検査の価値を明示する。
  3. 共通化と個別化を両立する:全社標準を共有しつつ、現場固有差を残す。
  4. 切片だけでなく感度差を見る:環境変動に弱いラインを予防保全へつなげる。

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

  • ライン、設備、品種、材料ロットなど階層キーを欠損なく記録する
  • 事前分布の根拠を過去期間・専門家判断・弱情報事前分布に分けて文書化する
  • 生産量差、検査方式変更、選別後データなどのバイアスを点検する
  • 事前分布の感度分析、事後予測チェック、期間外検証を実施する
  • 「改善」「追加測定」「監視継続」に対応する確率閾値と損失を合意する
  • モデル更新頻度、責任者、監査ログ、現場への説明資料を運用設計に含める

まとめ

階層ベイズの価値は、複雑なモデルそのものではなく、データ量の違いを公平に扱い、全体知と現場差を同時に意思決定へ載せることにあります。まずは共役モデルなど説明可能な小規模実装から始め、単純集計に対する改善を検証してからモデルを拡張するのが現実的です。

法人向けのご相談

数理工房では、品質データの階層設計、ベイズモデルのPoC、意思決定ルール設計、現場説明用レポートまで、業務課題に合わせてご支援します。

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