100本ノック / 確率統計 / 確率・統計理論100本ノック

製造業の不良・設備停止を離散分布で予測する｜Python実践10本ノック

2020-04-10

不良・停止は何件起こるか――離散分布で製造現場の検査・保全計画を立てる

概要

製造現場では、不良の有無、一定数の検査で見つかる不良数、次の異常までの間隔、設備アラームの件数など、整数で数える事象を日々扱います。本記事では、架空の電子部品工場を題材に、ベルヌーイ分布からポアソン過程まで10種類の離散分布を、品質保証・受入検査・保全要員配置の判断へ結び付けます。

対象は「確率・統計理論100本ノック」の No.041〜No.050 です。公式を使うだけでなく、「試行回数を固定するのか」「発生回数まで待つのか」「母集団から戻さずに抜き取るのか」という前提の違いをPythonで確かめます。

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

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

今回の題材は、センサー基板を1日400個生産する架空の工場です。品質保証部門と設備保全部門は、限られた人員で次の問いに答えなければなりません。

1個の製品が不良になる確率を、どのように最小単位として表すか
1日の不良個数や「5個以上発生する日」の確率をどう見積もるか
次の不良、または5件目の不良までに何個検査することになるか
有限ロットからの抜取検査では、二項分布をそのまま使ってよいか
設備アラームの件数と、時間内に対応能力を超える確率をどう評価するか
不良原因別の要員・部品・解析工数をどう配分するか
まれな不良の二項分布を、いつポアソン分布で近似できるか
件数だけでなく、発生時刻を含む保全計画をどう考えるか

離散分布は、現象に合った前提を選ぶことで、平均件数だけでは見えない超過確率や待ち時間を定量化します。

現場でよくある状況

現場では「不良率1.2%」「アラームは1時間平均0.65件」といった平均値が共有されます。しかし、日々の運用では次のような行き違いが生じます。

平均不良数だけを見て、上振れ日に必要な検査員を確保していない
抜取検査を、母集団へ戻して抽出する二項モデルで近似し、有限ロット補正を意識していない
「次の1件まで」と「5件発生するまで」を同じ分布で扱っている
不良原因別件数を独立な二項分布として計算し、合計件数が一致しない
1日件数がポアソン分布らしく見えるだけで、時間的な独立性や一定発生率まで成立したと判断する

分布名を覚えることより、何を固定し、何を確率変数とするかを決めることが重要です。

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

第一に、似た数式でも標本設計が異なります。二項分布は一定確率の独立試行を固定回数繰り返すモデルですが、超幾何分布は有限母集団から戻さずに抜き取るモデルです。抜取率が高いほど両者の差は無視できません。

第二に、管理したいKPIによって確率変数が変わります。「400個中の不良数」は二項分布、「最初の不良までの検査数」は幾何分布、「5件目までの検査数」は負の二項分布です。

第三に、ポアソン分布とポアソン過程には追加の仮定があります。一定時間内の件数がポアソン分布でも、季節性、設備状態、前回停止の影響があれば、一定発生率・独立増分という過程の仮定は崩れます。実務では適合度だけでなく、時間帯別推移や連続発生も確認する必要があります。

今回扱うノックの全体像

No.	テーマ	製造業での問い
041	ベルヌーイ分布	個品の良品・不良をどう表すか
042	二項分布	固定個数を検査したとき不良は何個か
043	幾何分布	最初の不良まで何個検査するか
044	負の二項分布	5件目の不良まで検査数はいくつか
045	超幾何分布	有限ロットの抜取受入確率はいくらか
046	ポアソン分布	一定時間の設備アラームは何件か
047	多項分布	不良総数を原因別にどう配分するか
048	カテゴリ分布	1件の不良原因をどう確率表現するか
049	ポアソン極限定理	まれな不良の二項分布を近似できるか
050	ポアソン過程	アラームの件数と発生時刻をどう扱うか

No.041を個品判定の最小単位とし、固定回数、待ち時間、非復元抽出、多分類、まれな事象、時間軸へ順に拡張します。

Python 環境の準備

NumPyで架空データとシミュレーションを生成し、pandasで表を作り、SciPyで理論確率を計算し、matplotlibで可視化します。外部データには依存せず、再実行して同じ結果になるよう乱数 seed を固定します。

import sys
import platform

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
from scipy.special import comb
from scipy.stats import binom, geom, hypergeom, poisson
import japanize_matplotlib  # noqa: F401  日本語フォント設定
from IPython.display import display

SEED = 20260711
rng = np.random.default_rng(SEED)

pd.set_option("display.max_columns", 20)
pd.set_option("display.float_format", lambda x: f"{x:,.4f}")
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.unicode_minus"] = False

print(f"Python: {sys.version.split()[0]}")
print(f"OS: {platform.system()} {platform.release()}")
print(f"NumPy: {np.__version__}")
print(f"pandas: {pd.__version__}")
print(f"SciPy: {scipy.__version__}")
print(f"matplotlib: {matplotlib.__version__}")
print(f"乱数 seed: {SEED}")

Python: 3.13.1
OS: Darwin 25.3.0
NumPy: 2.5.1
pandas: 3.0.3
SciPy: 1.18.0
matplotlib: 3.11.0
乱数 seed: 20260711

架空データの作成

60稼働日について、日産400個の不良数、8時間内の設備アラーム数、不良原因別件数を生成します。モデル上の基準値は次のとおりです。

個品不良確率： $p=0.012$ （1.2%）
設備アラーム発生率：1時間あたり $\lambda=0.65$ 件
不良原因構成：はんだ45%、部品25%、実装ずれ20%、その他10%
各日の試行は同一条件かつ独立と仮定する

最後の仮定は教材上の単純化です。実務データでは、設備号機、品種、材料ロット、勤務帯、レシピ変更、保全後の経過時間などで確率が変化していないかを先に確認します。

n_days = 60
units_per_day = 400
defect_probability = 0.012
hours_per_day = 8
alarm_rate_per_hour = 0.65
cause_names = np.array(["はんだ", "部品", "実装ずれ", "その他"])
cause_probabilities = np.array([0.45, 0.25, 0.20, 0.10])

daily_defects = rng.binomial(units_per_day, defect_probability, size=n_days)
daily_alarms = rng.poisson(hours_per_day * alarm_rate_per_hour, size=n_days)
daily_causes = np.vstack([
    rng.multinomial(count, cause_probabilities) for count in daily_defects
])

df = pd.DataFrame({
    "稼働日": pd.date_range("2026-04-01", periods=n_days, freq="D"),
    "生産数": units_per_day,
    "不良数": daily_defects,
    "不良率": daily_defects / units_per_day,
    "設備アラーム数": daily_alarms,
})
for i, cause in enumerate(cause_names):
    df[f"原因_{cause}"] = daily_causes[:, i]

display(df.head(10))
summary = pd.DataFrame({
    "指標": ["日次不良数", "日次不良率", "日次設備アラーム数"],
    "実績平均": [df["不良数"].mean(), df["不良率"].mean(), df["設備アラーム数"].mean()],
    "モデル期待値": [units_per_day * defect_probability, defect_probability,
                   hours_per_day * alarm_rate_per_hour],
})
display(summary.round(4))
print(f"欠損数: {int(df.isna().sum().sum())} / 総セル数: {df.size:,}")

	稼働日	生産数	不良数	不良率	設備アラーム数	原因_はんだ	原因_部品	原因_実装ずれ	原因_その他
0	2026-04-01	400	3	0.0075	5	2	0	1	0
1	2026-04-02	400	8	0.0200	7	5	0	2	1
2	2026-04-03	400	5	0.0125	3	1	2	0	2
3	2026-04-04	400	6	0.0150	5	4	0	0	2
4	2026-04-05	400	7	0.0175	5	1	3	3	0
5	2026-04-06	400	3	0.0075	9	2	1	0	0
6	2026-04-07	400	2	0.0050	3	0	2	0	0
7	2026-04-08	400	8	0.0200	2	3	2	3	0
8	2026-04-09	400	5	0.0125	7	0	2	3	0
9	2026-04-10	400	4	0.0100	9	1	0	3	0

	指標	実績平均	モデル期待値
0	日次不良数	4.8333	4.8000
1	日次不良率	0.0121	0.0120
2	日次設備アラーム数	5.1167	5.2000

欠損数: 0 / 総セル数: 540

No.041：ベルヌーイ分布

実務での意味

ベルヌーイ分布は、1個の検査結果を「不良=1、良品=0」の2値で表します。個品単位の合否判定、設備の停止有無、納期遵守の成否など、二者択一のKPIの最小モデルです。

分析・モデル化の考え方

$X\in\{0,1\}$ 、不良確率を $p$ とすると、確率質量関数、期待値、分散は

P(X=x)=p^x(1-p)^{1-x},\qquad E[X]=p,\qquad V[X]=p(1-p)

です。0/1データの標本平均は、そのまま不良率になります。ただし、製品ごとに $p$ が違う場合、全データを同じベルヌーイ試行とみなす前提は崩れます。

Pythonで確認する

bernoulli_table = pd.DataFrame({
    "判定_x": [0, 1],
    "意味": ["良品", "不良"],
    "確率_P(X=x)": [1 - defect_probability, defect_probability],
})
display(bernoulli_table)
print(f"理論期待値 E[X] = p: {defect_probability:.4f}")
print(f"理論分散 V[X] = p(1-p): {defect_probability * (1-defect_probability):.6f}")

plt.bar(bernoulli_table["意味"], bernoulli_table["確率_P(X=x)"], color=["#4C78A8", "#E45756"])
plt.title("個品検査のベルヌーイ分布")
plt.xlabel("検査判定")
plt.ylabel("確率")
plt.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	判定_x	意味	確率_P(X=x)
0	0	良品	0.9880
1	1	不良	0.0120

理論期待値 E[X] = p: 0.0120
理論分散 V[X] = p(1-p): 0.011856

png

結果の読み取り

1個の検査では良品か不良のどちらかしか観測されませんが、0/1を多数蓄積して平均すれば不良率を推定できます。不良確率1.2%は個品の不良を予言する値ではなく、同条件で多数生産したときの長期比率です。設備・品種・期間を混ぜる前に、同じ $p$ とみなせる管理単位を定義する必要があります。

No.042：二項分布

実務での意味

二項分布は、固定した生産数・検査数の中で不良が何個発生するかを表します。平均不良数だけでなく、再検査や応援要員が必要になる上振れ確率を計算できます。

分析・モデル化の考え方

独立なベルヌーイ試行を $n$ 回行い、不良数を $X$ とすると

X\sim\operatorname{Bin}(n,p),\quad P(X=k)=\binom{n}{k}p^k(1-p)^{n-k}

で、 $E[X]=np$ 、 $V[X]=np(1-p)$ です。日産400個、 $p=0.012$ とし、5個以上の不良が出る確率を求めます。

Pythonで確認する

k = np.arange(0, 16)
binom_pmf = binom.pmf(k, units_per_day, defect_probability)
tail_probability = binom.sf(4, units_per_day, defect_probability)

binom_table = pd.DataFrame({"不良数_k": k, "確率": binom_pmf})
display(binom_table.round(5))
print(f"期待不良数 np: {units_per_day * defect_probability:.2f} 個/日")
print(f"標準偏差: {np.sqrt(units_per_day * defect_probability * (1-defect_probability)):.2f} 個")
print(f"5個以上となる確率 P(X>=5): {tail_probability:.2%}")

plt.bar(k, binom_pmf, color="#4C78A8")
plt.axvline(4.5, color="#E45756", linestyle="--", label="5個以上")
plt.title("1日400個生産時の不良個数（二項分布）")
plt.xlabel("1日の不良個数")
plt.ylabel("確率")
plt.grid(axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	不良数_k	確率
0	0	0.0080
1	1	0.0388
2	2	0.0941
3	3	0.1516
4	4	0.1828
5	5	0.1758
6	6	0.1406
7	7	0.0961
8	8	0.0573
9	9	0.0303
10	10	0.0144
11	11	0.0062
12	12	0.0024
13	13	0.0009
14	14	0.0003
15	15	0.0001

期待不良数 np: 4.80 個/日
標準偏差: 2.18 個
5個以上となる確率 P(X>=5): 52.46%

png

結果の読み取り

期待値は4.8個ですが、5個以上になる確率は約半分あります。「5個は平均より多いから異常」とは言えません。日常管理の警戒線は平均との単純比較ではなく、許容できる超過頻度や損失額から設定します。また、連続する日で不良確率が変わるなら、単一の二項分布より品種・設備状態別のモデルが必要です。

No.043：幾何分布

実務での意味

幾何分布は、最初の不良が見つかるまでに何個検査するかを表します。連続良品数の目安、巡回検査間隔、異常を発見するまでの検査負荷を考える入口になります。

分析・モデル化の考え方

$T$ を最初の不良を含む検査個数とすると

P(T=t)=(1-p)^{t-1}p,\quad t=1,2,\ldots

で、 $E[T]=1/p$ です。幾何分布には「これまで良品が何個続いても、次の1個の不良確率は $p$ のまま」という無記憶性があります。工程劣化が進む現場では、この仮定は慎重に確認します。

Pythonで確認する

t = np.arange(1, 401)
geom_pmf = geom.pmf(t, defect_probability)
q50, q90, q95 = geom.ppf([0.50, 0.90, 0.95], defect_probability).astype(int)

geom_summary = pd.DataFrame({
    "指標": ["平均", "中央値", "90%分位点", "95%分位点"],
    "検査個数": [1 / defect_probability, q50, q90, q95],
})
display(geom_summary.round(1))

plt.plot(t, geom.cdf(t, defect_probability), color="#4C78A8", label="累積確率")
plt.axhline(0.90, color="#E45756", linestyle="--", label="90%")
plt.axvline(q90, color="#E45756", linestyle=":", label=f"{q90}個")
plt.title("最初の不良が見つかるまでの検査個数")
plt.xlabel("最初の不良を含む検査個数")
plt.ylabel("累積確率 P(T <= t)")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	指標	検査個数
0	平均	83.3000
1	中央値	58.0000
2	90%分位点	191.0000
3	95%分位点	249.0000

png

結果の読み取り

平均は約83個ですが、分布は右に長く、90%の確率で最初の不良が見つかる検査個数は平均よりかなり大きくなります。平均間隔だけで巡回検査を設計すると、発見確率を過大評価し得ます。検査間隔は「何%の確率で異常を捕捉したいか」というサービス水準から決めます。

No.044：負の二項分布

実務での意味

負の二項分布は、所定件数の不良が見つかるまでに必要な検査個数を表します。たとえば、原因解析用に不良品5個を確保するまでの検査工数や期間を見積もれます。

分析・モデル化の考え方

$T_r$ を $r$ 件目の不良を含む検査個数とすると

P(T_r=t)=\binom{t-1}{r-1}p^r(1-p)^{t-r},\quad t=r,r+1,\ldots

で、 $E[T_r]=r/p$ 、 $V[T_r]=r(1-p)/p^2$ です。SciPyの nbinom は「 $r$ 件目までの良品数」を返すため、総検査数に直すときは $r$ を足す点に注意します。

Pythonで確認する

r = 5
total_inspections = np.arange(r, 1201)
good_before_rth_defect = total_inspections - r
nb_pmf = scipy.stats.nbinom.pmf(good_before_rth_defect, r, defect_probability)
nb_quantiles = scipy.stats.nbinom.ppf([0.50, 0.90, 0.95], r, defect_probability) + r

nb_summary = pd.DataFrame({
    "指標": ["平均", "中央値", "90%分位点", "95%分位点"],
    "総検査個数": [r / defect_probability, *nb_quantiles],
})
display(nb_summary.round(1))

plt.plot(total_inspections, scipy.stats.nbinom.cdf(good_before_rth_defect, r, defect_probability),
         color="#59A14F")
plt.axhline(0.90, color="#E45756", linestyle="--", label="90%")
plt.axvline(nb_quantiles[1], color="#E45756", linestyle=":", label=f"{nb_quantiles[1]:.0f}個")
plt.title("5件目の不良を確保するまでの総検査個数")
plt.xlabel("総検査個数")
plt.ylabel("累積確率")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	指標	総検査個数
0	平均	416.7000
1	中央値	389.0000
2	90%分位点	665.0000
3	95%分位点	761.0000

png

結果の読み取り

5件確保までの平均は約417個ですが、90%の確度で確保するにはさらに多い検査能力が必要です。原因解析の着手日を平均だけで約束すると遅延しやすいため、分位点で計画します。なお、不良が集中する工程では独立・一定確率の仮定が崩れ、実際の分散はこのモデルより大きくなることがあります。

No.045：超幾何分布

実務での意味

超幾何分布は、有限ロットから製品を戻さずに抜き取る受入検査に適します。ロット内不良数が既知だと仮定したとき、抜取標本で不良を何個見つけるかを表します。

分析・モデル化の考え方

ロット総数 $N$ 、不良数 $D$ 、抜取数 $n$ 、発見不良数 $X$ に対して

P(X=k)=\frac{\binom{D}{k}\binom{N-D}{n-k}}{\binom{N}{n}}

です。ここでは $N=500$ 、 $D=8$ 、 $n=50$ とし、不良1個以下なら受入れる規則の受入確率を計算します。これは特定品質のロットを受け入れる確率であり、抜取方式の良否は複数の $D$ に対するOC曲線で評価します。

Pythonで確認する

N, D, sample_n, acceptance_c = 500, 8, 50, 1
x_h = np.arange(0, min(D, sample_n) + 1)
h_pmf = hypergeom.pmf(x_h, N, D, sample_n)
acceptance_probability = hypergeom.cdf(acceptance_c, N, D, sample_n)

display(pd.DataFrame({"抜取標本の不良数": x_h, "確率": h_pmf}).round(6))
print(f"不良1個以下で受入れる確率: {acceptance_probability:.2%}")

defect_counts = np.arange(0, 51)
oc_curve = hypergeom.cdf(acceptance_c, N, defect_counts, sample_n)
plt.plot(defect_counts / N * 100, oc_curve, color="#4C78A8")
plt.axvline(D / N * 100, color="#E45756", linestyle="--", label=f"不良率 {D/N:.1%}")
plt.scatter([D / N * 100], [acceptance_probability], color="#E45756", zorder=3)
plt.title("抜取検査方式のOC曲線（n=50、合格判定1個以下）")
plt.xlabel("ロット不良率（%）")
plt.ylabel("ロット受入確率")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	抜取標本の不良数	確率
0	0	0.4278
1	1	0.3862
2	2	0.1492
3	3	0.0322
4	4	0.0042
5	5	0.0003
6	6	0.0000
7	7	0.0000
8	8	0.0000

不良1個以下で受入れる確率: 81.40%

png

結果の読み取り

不良が8個含まれるロットでも、抜取結果が1個以下となり受け入れられる可能性があります。抜取検査は全数保証ではなく、良いロットを誤って拒否する生産者危険と、悪いロットを受け入れる消費者危険の配分です。許容品質水準だけでなく、重大不良の損失と検査費用を合わせて $n$ と合格判定個数を設計します。

No.046：ポアソン分布

実務での意味

ポアソン分布は、一定時間・面積・長さの中で発生する比較的まれな件数を表します。設備アラーム、異物欠点、問い合わせ、微小停止などの対応能力計画に使えます。

分析・モデル化の考え方

期間内の平均発生数を $\lambda$ とすると

P(X=k)=e^{-\lambda}\frac{\lambda^k}{k!},\qquad E[X]=V[X]=\lambda

です。本例では8時間あたり $\lambda=0.65\times8=5.2$ 件です。平均と分散が等しいことが特徴ですが、実績分散が大きければ発生率の変動や集中発生を疑います。

Pythonで確認する

lambda_day = alarm_rate_per_hour * hours_per_day
alarm_k = np.arange(0, 16)
alarm_pmf = poisson.pmf(alarm_k, lambda_day)
over_capacity_probability = poisson.sf(8, lambda_day)

poisson_check = pd.DataFrame({
    "項目": ["理論平均", "理論分散", "実績平均", "実績分散(ddof=1)", "9件以上の確率"],
    "値": [lambda_day, lambda_day, df["設備アラーム数"].mean(),
          df["設備アラーム数"].var(ddof=1), over_capacity_probability],
})
display(poisson_check.round(4))

plt.bar(alarm_k, alarm_pmf, color="#F28E2B")
plt.axvline(8.5, color="#E45756", linestyle="--", label="対応能力超過：9件以上")
plt.title("8時間シフトの設備アラーム件数（ポアソン分布）")
plt.xlabel("アラーム件数")
plt.ylabel("確率")
plt.grid(axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	項目	値
0	理論平均	5.2000
1	理論分散	5.2000
2	実績平均	5.1167
3	実績分散(ddof=1)	5.0879
4	9件以上の確率	0.0819

png

結果の読み取り

1シフト平均5.2件でも、9件以上となる日は一定確率で発生します。保全要員の能力を平均件数だけに合わせると、超過日に復旧待ちが累積します。一方、平均と分散の実績差は60日程度では揺れるため、直ちにモデル不適合とは断定せず、期間追加、設備別層別、連続発生の確認を行います。

No.047：多項分布

実務での意味

多項分布は、固定された総件数を、相互排他的な複数カテゴリへ配分するモデルです。不良総数を原因別に分け、解析担当者や交換部品の必要量を見積もる場面に対応します。

分析・モデル化の考え方

$n$ 件を $K$ カテゴリへ分類し、件数ベクトルを $\mathbf{X}=(X_1,\ldots,X_K)$ 、確率を $\mathbf{p}$ とすると

P(\mathbf{X}=\mathbf{x})=\frac{n!}{\prod_{j=1}^K x_j!}\prod_{j=1}^Kp_j^{x_j}, \quad \sum_jx_j=n,\quad \sum_jp_j=1

です。 $E[X_j]=np_j$ であり、総数が固定されるためカテゴリ件数間には負の共分散 $\operatorname{Cov}(X_i,X_j)=-np_ip_j$ があります。

Pythonで確認する

monthly_total = 120
n_simulations = 20_000
simulated_causes = rng.multinomial(monthly_total, cause_probabilities, size=n_simulations)

multinomial_summary = pd.DataFrame({
    "原因": cause_names,
    "原因確率": cause_probabilities,
    "理論期待件数": monthly_total * cause_probabilities,
    "シミュレーション平均": simulated_causes.mean(axis=0),
    "90%分位点": np.quantile(simulated_causes, 0.90, axis=0),
})
display(multinomial_summary.round(2))
print(f"各シミュレーションで原因件数の合計が{monthly_total}件: "
      f"{np.all(simulated_causes.sum(axis=1) == monthly_total)}")
print(f"はんだ件数と部品件数の経験共分散: {np.cov(simulated_causes[:, 0], simulated_causes[:, 1], ddof=0)[0,1]:.2f}")

xpos = np.arange(len(cause_names))
plt.bar(xpos - 0.18, monthly_total * cause_probabilities, width=0.36, label="理論期待件数", color="#4C78A8")
plt.bar(xpos + 0.18, np.quantile(simulated_causes, 0.90, axis=0), width=0.36,
        label="90%分位点", color="#F28E2B")
plt.xticks(xpos, cause_names)
plt.title("月120件の不良原因別件数")
plt.xlabel("不良原因")
plt.ylabel("月間件数")
plt.grid(axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	原因	原因確率	理論期待件数	シミュレーション平均	90%分位点
0	はんだ	0.4500	54.0000	53.9400	61.0000
1	部品	0.2500	30.0000	30.0100	36.0000
2	実装ずれ	0.2000	24.0000	24.0300	30.0000
3	その他	0.1000	12.0000	12.0200	16.0000

各シミュレーションで原因件数の合計が120件: True
はんだ件数と部品件数の経験共分散: -13.46

png

結果の読み取り

期待件数は人員・部品の基準量、90%分位点は上振れに備える容量の候補です。各原因件数を別々の独立二項分布として扱うと合計が120件を超えたり下回ったりします。総件数を固定して原因構成を扱うなら、多項分布の制約とカテゴリ間の依存を保つ必要があります。

No.048：カテゴリ分布

実務での意味

カテゴリ分布は、1件の不良を複数原因のいずれか1つへ分類するモデルです。多項分布が複数件の集計なら、カテゴリ分布はその1件分に相当します。自動原因分類モデルの出力確率も同じ形で表せます。

分析・モデル化の考え方

カテゴリ $j$ を表すone-hotベクトル $\mathbf{x}$ に対して

P(\mathbf{X}=\mathbf{x})=\prod_{j=1}^K p_j^{x_j}, \qquad x_j\in\{0,1\},\quad \sum_jx_j=1

です。カテゴリは数値の大小を持たないため、「はんだ=1、部品=2」の平均値に実務的意味はありません。確率またはone-hot表現で扱います。

Pythonで確認する

n_examples = 12
sampled_causes = rng.choice(cause_names, size=n_examples, p=cause_probabilities)
categorical_df = pd.DataFrame({"不良ID": [f"D{i:03d}" for i in range(1, n_examples + 1)],
                               "原因カテゴリ": sampled_causes})
one_hot = pd.get_dummies(categorical_df["原因カテゴリ"], dtype=int).reindex(columns=cause_names, fill_value=0)
display(pd.concat([categorical_df, one_hot], axis=1))

category_table = pd.DataFrame({"原因カテゴリ": cause_names, "確率": cause_probabilities})
plt.bar(category_table["原因カテゴリ"], category_table["確率"], color="#59A14F")
plt.title("1件の不良原因に対するカテゴリ分布")
plt.xlabel("不良原因")
plt.ylabel("選択確率")
plt.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	不良ID	原因カテゴリ	はんだ	部品	実装ずれ	その他
0	D001	その他	0	0	0	1
1	D002	部品	0	1	0	0
2	D003	部品	0	1	0	0
3	D004	はんだ	1	0	0	0
4	D005	はんだ	1	0	0	0
5	D006	はんだ	1	0	0	0
6	D007	実装ずれ	0	0	1	0
7	D008	はんだ	1	0	0	0
8	D009	実装ずれ	0	0	1	0
9	D010	はんだ	1	0	0	0
10	D011	はんだ	1	0	0	0
11	D012	その他	0	0	0	1

png

結果の読み取り

各行のone-hot列は必ず1つだけ1になり、同時に複数原因へ確定分類しない前提を表します。実務では複合原因や原因不明があり得るため、分類ルールを先に定義します。AI分類の確率を使う場合も、最大確率だけで確定せず、信頼度が低い案件を人手確認へ回す閾値設計が必要です。

No.049：ポアソン極限定理

実務での意味

ポアソン極限定理は、試行回数が大きく成功確率が小さい二項分布を、計算しやすいポアソン分布で近似できる根拠です。大量生産におけるまれな欠点件数の概算に役立ちます。

分析・モデル化の考え方

$n\to\infty$ 、 $p\to0$ 、 $np\to\lambda$ のとき

\operatorname{Bin}(n,p)\ \Rightarrow\ \operatorname{Poisson}(\lambda)

となります。有限の $n,p$ では近似誤差があるため、本例の $n=400,p=0.012,\lambda=4.8$ で確率質量関数と裾確率を比較します。平均 $np$ が同じでも $p$ が大きい場合は近似が悪くなります。

Pythonで確認する

lambda_defects = units_per_day * defect_probability
compare_k = np.arange(0, 16)
exact = binom.pmf(compare_k, units_per_day, defect_probability)
approx = poisson.pmf(compare_k, lambda_defects)
comparison = pd.DataFrame({
    "不良数": compare_k,
    "二項分布_厳密": exact,
    "ポアソン近似": approx,
    "絶対誤差": np.abs(exact - approx),
})
display(comparison.round(6))
print(f"表示範囲での最大絶対誤差: {comparison['絶対誤差'].max():.6f}")
print(f"P(X>=10) 二項分布: {binom.sf(9, units_per_day, defect_probability):.4%}")
print(f"P(X>=10) ポアソン近似: {poisson.sf(9, lambda_defects):.4%}")

plt.plot(compare_k, exact, marker="o", label="二項分布（厳密）", color="#4C78A8")
plt.plot(compare_k, approx, marker="s", linestyle="--", label="ポアソン近似", color="#E45756")
plt.title("まれな不良に対する二項分布とポアソン近似")
plt.xlabel("1日の不良個数")
plt.ylabel("確率")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	不良数	二項分布_厳密	ポアソン近似	絶対誤差
0	0	0.0080	0.0082	0.0002
1	1	0.0388	0.0395	0.0007
2	2	0.0941	0.0948	0.0007
3	3	0.1516	0.1517	0.0001
4	4	0.1828	0.1820	0.0008
5	5	0.1758	0.1747	0.0011
6	6	0.1406	0.1398	0.0008
7	7	0.0961	0.0959	0.0003
8	8	0.0574	0.0575	0.0002
9	9	0.0303	0.0307	0.0003
10	10	0.0144	0.0147	0.0003
11	11	0.0062	0.0064	0.0002
12	12	0.0024	0.0026	0.0001
13	13	0.0009	0.0009	0.0001
14	14	0.0003	0.0003	0.0000
15	15	0.0001	0.0001	0.0000

表示範囲での最大絶対誤差: 0.001093
P(X>=10) 二項分布: 2.4368%
P(X>=10) ポアソン近似: 2.5141%

png

結果の読み取り

この条件では中心部分の形はよく一致し、概算にはポアソン近似が使えます。ただし、受入判定や重大品質事故のように裾確率の小さな差が費用へ直結する場合は、計算可能なら二項分布の厳密値を使います。「 $n$ が大きい」だけでなく「 $p$ が十分小さい」「試行が概ね独立」という条件も確認します。

No.050：ポアソン過程

実務での意味

ポアソン過程は、一定期間の件数だけでなく、設備アラームがいつ発生するかを扱います。時間帯別の保全負荷、次の呼出しまでの余裕、応援要員が必要になる確率を考える基礎です。

分析・モデル化の考え方

発生率を $\lambda$ 件/時とする斉次ポアソン過程 $\{N(t)\}$ では

N(t)\sim\operatorname{Poisson}(\lambda t)

で、互いに重ならない時間区間の件数は独立です。また、到着間隔 $W$ は指数分布に従い、 $P(W>w)=e^{-\lambda w}$ 、 $E[W]=1/\lambda$ です。ここでは1シフトの発生時刻をシミュレーションし、時間別累積件数と理論期待値を比較します。

Pythonで確認する

arrival_rng = np.random.default_rng(SEED + 50)
interarrival_times = arrival_rng.exponential(scale=1 / alarm_rate_per_hour, size=100)
arrival_times = np.cumsum(interarrival_times)
arrival_times = arrival_times[arrival_times <= hours_per_day]

time_grid = np.linspace(0, hours_per_day, 161)
cumulative_counts = np.searchsorted(arrival_times, time_grid, side="right")
process_summary = pd.DataFrame({
    "指標": ["このシフトの発生件数", "理論期待件数", "平均到着間隔(時間)",
           "2時間以内に3件以上の確率"],
    "値": [len(arrival_times), alarm_rate_per_hour * hours_per_day, 1 / alarm_rate_per_hour,
          poisson.sf(2, alarm_rate_per_hour * 2)],
})
display(process_summary.round(4))
display(pd.DataFrame({"アラーム番号": np.arange(1, len(arrival_times) + 1),
                      "発生時刻_シフト開始後h": arrival_times}).round(3))

plt.step(time_grid, cumulative_counts, where="post", label="シミュレーション累積件数", color="#4C78A8")
plt.plot(time_grid, alarm_rate_per_hour * time_grid, linestyle="--", label="理論期待値 λt", color="#E45756")
plt.scatter(arrival_times, np.arange(1, len(arrival_times) + 1), color="#4C78A8", s=28)
plt.title("8時間シフトの設備アラーム到着過程")
plt.xlabel("シフト開始後の時間（h）")
plt.ylabel("累積アラーム件数")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	指標	値
0	このシフトの発生件数	9.0000
1	理論期待件数	5.2000
2	平均到着間隔(時間)	1.5385
3	2時間以内に3件以上の確率	0.1429

	アラーム番号	発生時刻_シフト開始後h
0	1	0.7790
1	2	2.5300
2	3	2.6130
3	4	2.7500
4	5	3.4490
5	6	3.5250
6	7	3.6930
7	8	3.8550
8	9	3.8830

png

結果の読み取り

実現した累積件数は期待直線の周りで階段状に動き、短時間に集中する区間もあれば長い空白もあります。平均到着間隔をそのまま「一定間隔で来る」と解釈してはいけません。実務では時刻データから時間帯別発生率、到着間隔、直前故障との依存、保全後のリセット効果を確認し、必要なら非斉次ポアソン過程や更新過程へ拡張します。

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

10の分布は別々の暗記事項ではなく、製造現場の問いを次の順で具体化する道具です。

観測単位を決める：個品の合否はベルヌーイ、固定個数中の件数は二項分布です。
固定する量を決める：試行回数固定なら二項分布、発生回数固定で試行数を問うなら幾何・負の二項分布です。
標本設計を反映する：有限ロットから非復元で抜き取るなら超幾何分布です。
分類構造を保つ：1件の分類はカテゴリ分布、固定総数の原因別件数は多項分布です。
時間と発生率を区別する：期間内件数はポアソン分布、到着時刻まで扱うならポアソン過程です。
平均ではなく超過確率で能力を決める：検査員数、保全要員、予備部品は、期待値だけでなく許容欠品率や対応遅延率に対応する分位点で設計します。

最も重要なのは、計算結果よりも仮定の監視です。品種混在、設備劣化、連鎖不良、時間帯差があると、一定確率・独立性・一定発生率は崩れます。

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

1. 計数基準と母数を固定する

不良「件数」と不良「個数」、手直し件数と廃棄件数を混在させないよう、判定単位、重複計上、再検査の扱いをデータ辞書にします。分母となる生産数・稼働時間・検査数も同じ粒度で記録します。

2. 時刻と層別キーを残す

設備号機、品種、材料ロット、勤務帯、作業者、レシピ、保全履歴と発生時刻を紐付けます。集計済み件数だけでは、発生確率の変化や集中発生を検証できません。

3. 分布仮定を定期的に診断する

平均と分散、実績度数と理論度数、時系列プロット、到着間隔、自己相関を確認します。過分散があれば、負の二項回帰、階層モデル、混合分布などを検討します。

4. 確率を損失とサービス水準へ接続する

「9件以上の確率」だけでは意思決定になりません。1件を待たせた停止損失、応援要員費、見逃し損失、誤判定損失を整理し、許容超過確率と対応能力を合意します。

5. 運用後の変化を監視する

工程改善や設備更新後は $p$ や $\lambda$ が変わります。モデルの前提期間、更新頻度、責任者、警報時のアクションを標準作業として定めます。

まとめ

ベルヌーイ分布は個品の0/1判定、二項分布は固定試行回数内の件数を表す
幾何分布は最初の発生まで、負の二項分布は所定件数の発生までの試行数を表す
有限ロットの非復元抜取には超幾何分布を使う
ポアソン分布は期間内件数、多項・カテゴリ分布は原因分類、ポアソン過程は発生時刻を扱う
ポアソン近似は便利だが、まれな事象・独立性・近似誤差を確認する
人員や検査能力は平均件数だけでなく、分位点、超過確率、損失額から決める
実務適用では、一定確率・独立性・一定発生率という仮定を、設備・品種・時間帯別に監視する

離散分布を正しく選ぶことは、単なる確率計算ではなく、品質保証の強さと運用コストを同じ土俵で比較するための設計作業です。

法人向けのご相談

数理工房では、製造業のお客様に対して、次のようなご支援を行っています。

不良・停止・故障データの定義整理と分析基盤設計
抜取検査方式、OC曲線、検査能力の定量評価
設備アラームや微小停止の発生モデルと保全要員計画
品種・設備・勤務帯を考慮した階層的な不良率分析
Python notebookを用いた企業研修、分析プロトタイプ、現場実装支援

確率分布を現場KPIへどう結び付けるか、既存データでどこまで判断できるかという段階からご相談いただけます。

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

製造業の不良・設備停止を離散分布で予測する｜Python実践10本ノック

不良・停止は何件起こるか――離散分布で製造現場の検査・保全計画を立てる

概要

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.041：ベルヌーイ分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.042：二項分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.043：幾何分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.044：負の二項分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.045：超幾何分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.046：ポアソン分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.047：多項分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.048：カテゴリ分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.049：ポアソン極限定理

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.050：ポアソン過程

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

1. 計数基準と母数を固定する

2. 時刻と層別キーを残す

3. 分布仮定を定期的に診断する

4. 確率を損失とサービス水準へ接続する

5. 運用後の変化を監視する

まとめ

法人向けのご相談

Related Posts

工程変更を「効いた気がする」で終わらせない――統計的推論で品質投資を判断する

製造業の統計推定をPythonで学ぶ｜抜取検査からEMアルゴリズムまで

標本数を、品質判断の根拠に変える――製造KPIの極限定理10本ノック

複数KPIを同時に読む――製造品質・設備負荷・納期リスクの多変量確率

故障・品質・保全時間をどう見積もるか――製造業で使う連続分布10選