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

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

2020-04-15

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

製造現場では、日々の平均不良率や平均サイクルタイムを見ながら、工程停止、追加検査、設備投資を判断します。しかし、少数データの平均は偶然に振れ、分布を仮定しない安全側評価は保守的になりがちです。本記事では、架空の精密部品工場を題材に、**No.071〜No.080（極限定理）**を、必要標本数、管理閾値、推定精度の判断へ結び付けます。

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

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

舞台は、自動車向け精密シャフトを加工する架空工場です。品質保証部門は抜取検査から工程不良率を推定し、生産管理部門はサイクルタイムと日次損失を監視しています。経営会議では、次の問いへの説明が必要です。

何ロット観測すれば、平均KPIを安定した値として扱えるか
元データが歪んでいても、標本平均を正規近似してよいか
母分散が未知でも標準化や誤差評価ができるか
分布を決め切れない状況で、閾値超過をどこまで安全側に評価できるか
不良率推定値の誤差を、納入数量や品質損失の誤差へどう変換するか

現場でよくある状況

月初の会議で「直近10ロットの平均が悪化した」と工程停止を求める一方、別の帳票では「100ロット平均なら問題ない」と報告されることがあります。また、3%という不良率だけが共有され、標本数が100個なのか1万個なのかが省略されることもあります。

平均値は同じでも標本数が違えば精度は違います。逆に、安全を重視して極端に保守的な上界だけを使うと、過剰検査や過剰在庫を招きます。データ数、変数の範囲、独立性、分散について何が分かっているかを整理し、目的に合う理論を選ぶ必要があります。

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

極限定理の多くは「標本数が十分大きいとき」の性質を述べますが、十分の大きさは、元分布の歪み、裾の重さ、自己相関、求める精度によって変わります。 $n=30$ なら常に正規近似できる、という一律の規則ではありません。

さらに、チェビシェフ、マルコフ、Chernoff、Hoeffding などの不等式は、必要とする仮定と上界の鋭さが異なります。上界は実際の確率そのものではなく、「仮定の下でこれを超えない」という保証です。保証値と予測値を区別し、保守性に伴うコストも評価することが実務上の要点です。

今回扱うノックの全体像

No.	テーマ	製造業での判断
071	大数の法則	累積平均が安定するまでの観測量を確認する
072	中心極限定理	歪んだ個票データから標本平均の誤差を近似する
073	スラツキーの定理	未知の母標準偏差を標本標準偏差で置き換える
074	デルタ法	不良率の誤差を必要投入数の誤差へ変換する
075	チェビシェフの不等式	平均・分散だけで逸脱確率を上から抑える
076	マルコフの不等式	非負の品質損失の高額化リスクを評価する
077	Chernoff境界	二項不良個数の上側確率を指数的に評価する
078	Hoeffding不等式	有界な検査結果の平均誤差を標本数へ結び付ける
079	Concentration Inequality	複数の集中不等式の仮定と保守性を比較する
080	漸近正規性	不良率推定量の分布と信頼区間精度を確認する

Python 環境の準備

NumPy、pandas、SciPy、Matplotlibを利用します。外部データには依存せず、乱数 seed は 20260711 に固定します。理論値とシミュレーションを並べ、近似・上界・実測頻度を混同しない表示にします。

import platform
import sys

import japanize_matplotlib  # noqa: F401  日本語フォントを有効化
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
from IPython.display import display
from scipy.stats import binom, norm

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

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"Platform   : {platform.platform()}")

Python     : 3.13.1
NumPy      : 2.5.1
pandas     : 3.0.3
SciPy      : 1.18.0
Matplotlib : 3.11.0
Platform   : macOS-26.3-arm64-arm-64bit-Mach-O

架空データの作成

2,400ロットについて、検査数、不良個数、平均サイクルタイム、品質損失を生成します。不良率には製品差と日々の揺れを持たせ、品質損失は右に裾の長い分布とします。後続の定理を明確に確認する箇所では、この表とは別に条件を固定した反復シミュレーションも行います。

n_lots = 2_400
products = rng.choice(["製品A", "製品B", "製品C"], n_lots, p=[0.50, 0.30, 0.20])
base_p = pd.Series(products).map({"製品A": 0.018, "製品B": 0.030, "製品C": 0.048}).to_numpy()
inspection_n = rng.choice([80, 100, 120], n_lots, p=[0.2, 0.6, 0.2])
lot_p = np.clip(base_p * rng.lognormal(0, 0.18, n_lots), 0.002, 0.15)
defects = rng.binomial(inspection_n, lot_p)
cycle = rng.lognormal(np.log(61), 0.10, n_lots)
quality_loss = 6_000 + 9_000 * defects + rng.lognormal(np.log(8_000), 0.75, n_lots)

df = pd.DataFrame({
    "ロットID": [f"L{i:04d}" for i in range(1, n_lots + 1)],
    "製品": products,
    "検査数": inspection_n,
    "不良個数": defects,
    "不良率": defects / inspection_n,
    "平均サイクルタイム_sec": cycle,
    "品質損失_円": quality_loss,
})
display(df.head())
display(df.groupby("製品", observed=True).agg(
    ロット数=("ロットID", "size"),
    平均不良率=("不良率", "mean"),
    平均サイクルタイム_sec=("平均サイクルタイム_sec", "mean"),
    平均品質損失_円=("品質損失_円", "mean"),
).round(3))

	ロットID	製品	検査数	不良個数	不良率	平均サイクルタイム_sec	品質損失_円
0	L0001	製品A	100	1	0.01	60.847742	18462.889513
1	L0002	製品C	100	5	0.05	75.243633	66797.020735
2	L0003	製品B	100	3	0.03	62.494952	40788.067078
3	L0004	製品B	100	2	0.02	56.958848	31125.459117
4	L0005	製品C	100	5	0.05	59.491457	63804.395203

	ロット数	平均不良率	平均サイクルタイム_sec	平均品質損失_円
製品
製品A	1228	0.018	61.560	33018.239
製品B	689	0.030	60.998	44068.843
製品C	483	0.048	61.182	59565.582

No.071：大数の法則

実務での意味

大数の法則は、独立で同じ分布に従う観測を重ねると、標本平均が母平均へ近づくことを示します。日次品質損失の累積平均が安定する様子を見れば、少数日の実績だけで年間予算を決める危険を説明できます。

分析・モデル化の考え方

期待値 $E[X]=\mu$ が存在する独立同分布の確率変数 $X_1,\ldots,X_n$ に対し、標本平均

\bar{X}_n=\frac{1}{n}\sum_{i=1}^n X_i

は $n\to\infty$ で $\mu$ へ確率収束します。これは平均の誤差が必ず単調に減るという意味ではなく、途中で上下しながら長期的に安定する性質です。

Pythonで確認する

x = df["品質損失_円"].to_numpy()
running_mean = np.cumsum(x) / np.arange(1, len(x) + 1)
reference_mean = x.mean()
checkpoints = [10, 30, 100, 300, 1_000, 2_400]
lln_table = pd.DataFrame({
    "観測ロット数": checkpoints,
    "累積平均損失_円": [running_mean[i - 1] for i in checkpoints],
    "全期間平均との差_pct": [100 * (running_mean[i - 1] / reference_mean - 1) for i in checkpoints],
})
display(lln_table.round(2))

plt.plot(np.arange(1, len(x) + 1), running_mean, label="累積平均")
plt.axhline(reference_mean, color="tab:red", linestyle="--", label="2,400ロット平均")
plt.title("観測ロット数と累積平均品質損失")
plt.xlabel("観測ロット数")
plt.ylabel("累積平均品質損失（円/ロット）")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	観測ロット数	累積平均損失_円	全期間平均との差_pct
0	10	45265.58	8.99
1	30	41531.88	-0.00
2	100	41717.51	0.44
3	300	41504.03	-0.07
4	1000	40694.68	-2.02
5	2400	41533.34	-0.00

png

結果の読み取り

初期の累積平均は高額損失ロットの影響で動きますが、観測数が増えると2,400ロット平均の周辺へ安定します。ただし、大数の法則は「現在の平均が将来も不変」とは保証しません。設備改造、製品構成、検査基準が変わった場合は同一分布という前提が崩れます。累積平均に加え、層別・時系列・工程変更点を管理する必要があります。

No.072：中心極限定理

実務での意味

個々のサイクルタイムが右に歪んでいても、複数ロットの平均は標本数を増やすと正規分布に近づきます。これにより、週次平均の誤差幅や能力計画の余裕を正規近似で説明できます。

分析・モデル化の考え方

平均 $\mu$ 、有限な分散 $\sigma^2$ を持つ独立同分布標本では、

\frac{\sqrt{n}(\bar{X}_n-\mu)}{\sigma}\xrightarrow{d}N(0,1)

が成り立ちます。標本平均の標準偏差は $\sigma/\sqrt{n}$ なので、誤差を半分にするには原則として標本数を4倍にします。

Pythonで確認する

clt_rng = np.random.default_rng(SEED + 72)
mu_log, sigma_log = np.log(61), 0.28
population_mean = np.exp(mu_log + sigma_log**2 / 2)
population_sd = np.sqrt((np.exp(sigma_log**2) - 1) * np.exp(2 * mu_log + sigma_log**2))
n_reps = 8_000

fig, axes = plt.subplots(1, 3, figsize=(13, 3.8))
clt_rows = []
for ax, n in zip(axes, [5, 30, 100]):
    means = clt_rng.lognormal(mu_log, sigma_log, size=(n_reps, n)).mean(axis=1)
    xs = np.linspace(means.min(), means.max(), 300)
    ax.hist(means, bins=40, density=True, alpha=0.65, color="tab:blue")
    ax.plot(xs, norm.pdf(xs, population_mean, population_sd / np.sqrt(n)), color="tab:red")
    ax.set_title(f"標本平均（n={n}）")
    ax.set_xlabel("平均サイクルタイム（秒）")
    ax.set_ylabel("密度")
    ax.grid(alpha=0.3)
    clt_rows.append([n, means.mean(), means.std(ddof=1), population_sd / np.sqrt(n)])
plt.tight_layout()
plt.show()
display(pd.DataFrame(clt_rows, columns=["n", "標本平均の平均", "標本平均の実測SD", "理論標準誤差"]).round(3))

png

	n	標本平均の平均	標本平均の実測SD	理論標準誤差
0	5	63.403	8.042	8.102
1	30	63.406	3.315	3.308
2	100	63.500	1.829	1.812

結果の読み取り

元の対数正規分布は右に歪んでいますが、標本平均の分布は $n$ の増加とともに左右対称な正規曲線へ近づき、実測SDも理論標準誤差へ合います。小標本や極端に裾の重い分布では近似が遅いことがあります。また、連続ロットに自己相関があれば見かけの標本数ほど情報が増えないため、時系列構造も確認します。

No.073：スラツキーの定理

実務での意味

中心極限定理には母標準偏差 $\sigma$ が現れますが、実務では未知です。スラツキーの定理により、一致推定量である標本標準偏差 $S$ へ置き換えても、大標本では同じ標準正規分布へ近づくと説明できます。

分析・モデル化の考え方

$Z_n\xrightarrow{d}Z$ 、 $S_n\xrightarrow{p}\sigma$ なら、連続な四則演算の下で組み合わせた量も対応する分布へ収束します。したがって、

\frac{\sqrt{n}(\bar X_n-\mu)}{S_n}\xrightarrow{d}N(0,1)

です。有限標本での厳密性を主張する定理ではない点に注意します。

Pythonで確認する

sl_rng = np.random.default_rng(SEED + 73)
n, reps = 80, 10_000
samples = sl_rng.lognormal(mu_log, sigma_log, size=(reps, n))
means = samples.mean(axis=1)
s = samples.std(axis=1, ddof=1)
z_known = np.sqrt(n) * (means - population_mean) / population_sd
z_estimated = np.sqrt(n) * (means - population_mean) / s

summary = pd.DataFrame({
    "標準化": ["母SDを使用", "標本SDを使用", "標準正規"],
    "平均": [z_known.mean(), z_estimated.mean(), 0],
    "標準偏差": [z_known.std(ddof=1), z_estimated.std(ddof=1), 1],
    "2.5%点": [np.quantile(z_known, .025), np.quantile(z_estimated, .025), norm.ppf(.025)],
    "97.5%点": [np.quantile(z_known, .975), np.quantile(z_estimated, .975), norm.ppf(.975)],
})
display(summary.round(3))

xs = np.linspace(-4, 4, 300)
plt.hist(z_estimated, bins=50, density=True, alpha=0.65, label="標本SDで標準化")
plt.plot(xs, norm.pdf(xs), color="tab:red", label="標準正規密度")
plt.title("未知の母標準偏差を標本標準偏差で置き換えた統計量")
plt.xlabel("標準化統計量")
plt.ylabel("密度")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	標準化	平均	標準偏差	2.5%点	97.5%点
0	母SDを使用	-0.008	1.002	-1.95	1.998
1	標本SDを使用	-0.058	1.027	-2.19	1.873
2	標準正規	0.000	1.000	-1.96	1.960

png

結果の読み取り

標本SDで標準化した統計量も、平均0・標準偏差1の標準正規分布に近い形になります。これが、未知のばらつきを実績から推定して標準誤差を作る理論的な支えです。ただし小標本では置換誤差を無視できず、正規母集団の平均なら t 分布を使うなど、有限標本に適した方法を優先します。

No.074：デルタ法

実務での意味

不良率そのものではなく、「良品1,000個を納入するための期待投入数」を報告したい場面があります。デルタ法は、推定した不良率の誤差を、非線形に変換したKPIの誤差へ近似的に伝播させます。

分析・モデル化の考え方

$\sqrt{n}(\hat\theta-\theta)\xrightarrow{d}N(0,V)$ で、 $g$ が微分可能なら、

\sqrt{n}\{g(\hat\theta)-g(\theta)\}\xrightarrow{d} N\left(0,\{g'(\theta)\}^2V\right)

です。ここでは不良率 $p$ に対し $g(p)=1000/(1-p)$ 、 $g'(p)=1000/(1-p)^2$ とします。

Pythonで確認する

delta_rng = np.random.default_rng(SEED + 74)
p_true, n, reps = 0.04, 800, 50_000
p_hat = delta_rng.binomial(n, p_true, reps) / n
required_input = 1_000 / (1 - p_hat)
g_true = 1_000 / (1 - p_true)
se_p = np.sqrt(p_true * (1 - p_true) / n)
delta_se = 1_000 / (1 - p_true) ** 2 * se_p

delta_table = pd.DataFrame({
    "指標": ["変換KPIの平均", "変換KPIの標準偏差", "95%下限", "95%上限"],
    "シミュレーション": [required_input.mean(), required_input.std(ddof=1), *np.quantile(required_input, [.025, .975])],
    "デルタ法": [g_true, delta_se, g_true - 1.96 * delta_se, g_true + 1.96 * delta_se],
})
display(delta_table.round(2))

plt.hist(required_input, bins=45, density=True, alpha=0.65, label="シミュレーション")
xs = np.linspace(required_input.min(), required_input.max(), 300)
plt.plot(xs, norm.pdf(xs, g_true, delta_se), color="tab:red", label="デルタ法の正規近似")
plt.title("不良率推定誤差を必要投入数へ伝播")
plt.xlabel("良品1,000個に対する必要投入数（個）")
plt.ylabel("密度")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	指標	シミュレーション	デルタ法
0	変換KPIの平均	1041.78	1041.67
1	変換KPIの標準偏差	7.50	7.52
2	95%下限	1028.28	1026.93
3	95%上限	1056.80	1056.40

png

結果の読み取り

デルタ法の標準誤差と95%範囲は、反復シミュレーションに近い値です。不良率の点推定だけでなく、必要投入数にも幅を付けることで、材料手配や能力余力を議論できます。変換が急激に曲がる領域、境界付近、小標本では正規近似が崩れやすいため、ブートストラップや直接シミュレーションと比較します。

No.075：チェビシェフの不等式

実務での意味

工程KPIの分布形を特定できなくても、平均と有限な分散が分かれば、平均から大きく外れる確率を上から抑えられます。新工程で分布仮定の根拠が弱い時期の安全側評価に使えます。

分析・モデル化の考え方

平均 $\mu$ 、標準偏差 $\sigma$ を持つ任意の確率変数について、 $k>0$ なら

P(|X-\mu|\ge k\sigma)\le \frac{1}{k^2}

です。分布をほとんど仮定しない代わりに、上界は一般に保守的です。

Pythonで確認する

loss = df["品質損失_円"].to_numpy()
mu_loss, sd_loss = loss.mean(), loss.std(ddof=0)
k_values = np.array([1.5, 2.0, 2.5, 3.0, 4.0])
empirical = np.array([np.mean(np.abs(loss - mu_loss) >= k * sd_loss) for k in k_values])
cheb = 1 / k_values**2
cheb_table = pd.DataFrame({"k": k_values, "実測逸脱率": empirical, "チェビシェフ上界": cheb})
display(cheb_table.round(4))

plt.plot(k_values, empirical, marker="o", label="実測逸脱率")
plt.plot(k_values, cheb, marker="s", label="チェビシェフ上界")
plt.title("平均からk標準偏差以上外れる確率")
plt.xlabel("k（標準偏差の倍数）")
plt.ylabel("確率・上界")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	k	実測逸脱率	チェビシェフ上界
0	1.5	0.0996	0.4444
1	2.0	0.0450	0.2500
2	2.5	0.0208	0.1600
3	3.0	0.0100	0.1111
4	4.0	0.0025	0.0625

png

結果の読み取り

すべての閾値で実測逸脱率はチェビシェフ上界以下ですが、差は大きく、上界を予測値として使うとリスクを過大評価します。この不等式は「最低限の保証」を作る用途に向きます。実務では、十分な履歴が得られた後に分布適合や再標本化を併用し、保守性による追加検査・在庫コストも明示します。

No.076：マルコフの不等式

実務での意味

品質損失、停止時間、手直し工数のような非負KPIでは、平均だけから高額・長時間となる確率の上界を作れます。詳細分布がまだない新ラインで、暫定的なリスク枠を置く際に有効です。

分析・モデル化の考え方

非負確率変数 $X$ と $a>0$ に対し、

P(X\ge a)\le\frac{E[X]}{a}

がマルコフの不等式です。非負性だけを使うため適用範囲は広い一方、上界が1を超える場合は自明な上界1へ切り詰めます。

Pythonで確認する

thresholds = np.array([40_000, 60_000, 80_000, 120_000, 180_000])
markov_emp = np.array([(loss >= a).mean() for a in thresholds])
markov_bound = np.minimum(1, mu_loss / thresholds)
markov_table = pd.DataFrame({
    "損失閾値_円": thresholds,
    "実測超過率": markov_emp,
    "マルコフ上界": markov_bound,
    "上界から見た1000ロット中の最大件数": np.ceil(1_000 * markov_bound).astype(int),
})
display(markov_table.round(4))

plt.plot(thresholds, markov_emp, marker="o", label="実測超過率")
plt.plot(thresholds, markov_bound, marker="s", label="マルコフ上界")
plt.title("品質損失の閾値超過率とマルコフ上界")
plt.xlabel("品質損失の閾値（円/ロット）")
plt.ylabel("超過確率・上界")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	損失閾値_円	実測超過率	マルコフ上界	上界から見た1000ロット中の最大件数
0	40000	0.4483	1.0000	1000
1	60000	0.1721	0.6922	693
2	80000	0.0542	0.5192	520
3	120000	0.0029	0.3461	347
4	180000	0.0000	0.2307	231

png

結果の読み取り

マルコフ上界は実測超過率を確かに上回りますが、かなり保守的です。平均損失だけしか共有されていない段階でも上限を示せることが価値であり、精密な予測をする手法ではありません。KPIに負値があり得る場合はそのまま適用できません。保証予算へ使うなら、分散や分布情報を追加した手法との費用差を比較します。

No.077：Chernoff境界

実務での意味

一定不良率を仮定したロットで、不良個数が期待値を大きく超える確率を指数関数的な上界で評価します。全数シミュレーションをせずに、警報閾値のリスクオーダーを把握できます。

分析・モデル化の考え方

$X\sim\mathrm{Binomial}(n,p)$ 、 $\mu=np$ 、 $\delta>0$ に対する乗法型Chernoff境界は

P\{X\ge(1+\delta)\mu\}\le \left(\frac{e^\delta}{(1+\delta)^{1+\delta}}\right)^\mu

です。指数モーメントを使うため、平均・分散だけの一般的な上界より急速に小さくなります。

Pythonで確認する

n, p = 1_000, 0.02
mu = n * p
thresholds = np.array([25, 30, 35, 40, 45])
deltas = thresholds / mu - 1
exact_tail = binom.sf(thresholds - 1, n, p)
chernoff = (np.exp(deltas) / (1 + deltas) ** (1 + deltas)) ** mu
chernoff_table = pd.DataFrame({
    "警報閾値_不良個数": thresholds,
    "期待値の何倍": thresholds / mu,
    "二項分布の厳密確率": exact_tail,
    "Chernoff上界": chernoff,
})
display(chernoff_table)

plt.semilogy(thresholds, exact_tail, marker="o", label="二項分布の厳密確率")
plt.semilogy(thresholds, chernoff, marker="s", label="Chernoff上界")
plt.title("ロット不良個数の上側確率")
plt.xlabel("警報閾値（不良個数/1,000個）")
plt.ylabel("超過確率・上界（対数軸）")
plt.grid(alpha=0.3, which="both")
plt.legend()
plt.tight_layout()
plt.show()

	警報閾値_不良個数	期待値の何倍	二項分布の厳密確率	Chernoff上界
0	25	1.25	1.545154e-01	0.560689
1	30	1.50	2.069652e-02	0.114870
2	35	1.75	1.326559e-03	0.010188
3	40	2.00	4.339876e-05	0.000441
4	45	2.25	7.717810e-07	0.000010

png

結果の読み取り

警報閾値が期待不良個数20個から離れるほど、厳密確率とChernoff上界は指数的に低下します。上界は厳密確率より大きいものの、極小確率の桁を高速に評価できます。ただし、個体ごとの不良が独立で一定確率という仮定が必要です。設備異常により不良が群発する工程では、二項モデル自体が過小評価になるため、過分散やロット内相関を調べます。

No.078：Hoeffding不等式

実務での意味

良品・不良の0/1データは必ず $[0,1]$ に収まります。Hoeffding不等式を使うと、真の不良率から標本不良率が一定以上ずれる確率を、母分布の形や母分散を使わず標本数へ結び付けられます。

分析・モデル化の考え方

独立な $X_i\in[0,1]$ について、

P(|\bar X-E[\bar X]|\ge\varepsilon)\le2\exp(-2n\varepsilon^2)

です。右辺を許容リスク $\alpha$ 以下にすれば、 $n\ge\log(2/\alpha)/(2\varepsilon^2)$ という分布非依存の標本数設計ができます。

Pythonで確認する

hf_rng = np.random.default_rng(SEED + 78)
p, eps, reps = 0.03, 0.02, 80_000
n_values = np.array([500, 1_000, 2_000, 5_000, 10_000])
empirical = []
for n in n_values:
    phat = hf_rng.binomial(n, p, reps) / n
    empirical.append(np.mean(np.abs(phat - p) >= eps))
hoeffding = np.minimum(1, 2 * np.exp(-2 * n_values * eps**2))
hf_table = pd.DataFrame({"検査数_n": n_values, "実測誤差超過率": empirical, "Hoeffding上界": hoeffding})
display(hf_table)

plt.semilogy(n_values, np.maximum(empirical, 1 / reps), marker="o", label="実測誤差超過率")
plt.semilogy(n_values, hoeffding, marker="s", label="Hoeffding上界")
plt.title("標本不良率が真値から2ポイント以上ずれる確率")
plt.xlabel("検査数 n")
plt.ylabel("確率・上界（対数軸）")
plt.grid(alpha=0.3, which="both")
plt.legend()
plt.tight_layout()
plt.show()

alpha = 0.05
n_required = int(np.ceil(np.log(2 / alpha) / (2 * eps**2)))
print(f"誤差±{eps:.1%}を95%以上で保証するHoeffding標本数: {n_required:,}個")

	検査数_n	実測誤差超過率	Hoeffding上界
0	500	0.010700	1.000000
1	1000	0.000375	0.898658
2	2000	0.000000	0.403793
3	5000	0.000000	0.036631
4	10000	0.000000	0.000671

png

誤差±2.0%を95%以上で保証するHoeffding標本数: 4,612個

結果の読み取り

標本数とともに実測誤差率も上界も低下します。Hoeffding標本数は真の不良率を使わないため、工程立上げ時にも設計できますが、低不良率工程では保守的です。連続生産で自己相関がある、検査対象を無作為に抽出していない、測定誤判定がある場合は、形式的に標本数を満たしても保証の前提が崩れます。

No.079：Concentration Inequality（集中不等式）

実務での意味

集中不等式は、標本平均などが期待値の周囲へどれほど集中するかを保証する不等式群の総称です。同じ品質データでも、平均・分散・値の範囲のどこまで信頼できるかに応じて、チェビシェフ、Hoeffding、Bernsteinなどを使い分けます。

分析・モデル化の考え方

独立なBernoulli変数の平均に対する二側Bernstein型上界を

P(|\bar X-p|\ge\varepsilon)\le 2\exp\left\{-\frac{n\varepsilon^2}{2p(1-p)+2\varepsilon/3}\right\}

とします。分散情報を使うため、範囲だけを使うHoeffdingより鋭くなる場合があります。一方、チェビシェフは独立性や有界性を要求しない代わりに緩い上界です。

Pythonで確認する

p, n = 0.03, 2_000
eps_values = np.array([0.005, 0.010, 0.015, 0.020, 0.025])
ci_rng = np.random.default_rng(SEED + 79)
phat = ci_rng.binomial(n, p, 200_000) / n
emp = np.array([np.mean(np.abs(phat - p) >= e) for e in eps_values])
cheb = np.minimum(1, p * (1 - p) / (n * eps_values**2))
hoeff = np.minimum(1, 2 * np.exp(-2 * n * eps_values**2))
bern = np.minimum(1, 2 * np.exp(-n * eps_values**2 / (2 * p * (1 - p) + 2 * eps_values / 3)))

concentration = pd.DataFrame({
    "許容誤差_eps": eps_values,
    "実測確率": emp,
    "Chebyshev": cheb,
    "Hoeffding": hoeff,
    "Bernstein": bern,
})
display(concentration)

for col, marker in [("実測確率", "o"), ("Chebyshev", "s"), ("Hoeffding", "^"), ("Bernstein", "D")]:
    plt.semilogy(eps_values, np.maximum(concentration[col], 1e-6), marker=marker, label=col)
plt.title("集中不等式の仮定と上界の比較")
plt.xlabel("不良率の許容誤差 ε")
plt.ylabel("誤差超過確率・上界（対数軸）")
plt.grid(alpha=0.3, which="both")
plt.legend()
plt.tight_layout()
plt.show()

	許容誤差_eps	実測確率	Chebyshev	Hoeffding	Bernstein
0	0.005	0.189615	0.582000	1.000000	8.874345e-01
1	0.010	0.009295	0.145500	1.000000	9.162047e-02
2	0.015	0.000160	0.064667	0.813139	2.725528e-03
3	0.020	0.000000	0.036375	0.403793	2.780068e-05
4	0.025	0.000000	0.023280	0.164170	1.121754e-07

png

結果の読み取り

どの上界も実測確率以上ですが、使う情報が増えるほど一般に上界は鋭くなります。チェビシェフは仮定が少ない、Hoeffdingは範囲と独立性を使う、Bernsteinはさらに分散を使う、という交換関係です。最も小さい数字を選ぶのではなく、現場データが仮定を満たすと説明できる範囲で選び、保証値・経験頻度・モデル予測を別列で報告します。

No.080：漸近正規性

実務での意味

多くの推定量は、標本数が増えると真値の周りで正規分布に近づきます。不良率の最尤推定量 $\hat p$ についてこの性質を確認すると、大標本の信頼区間や工場間比較に正規近似を使う根拠と限界が分かります。

分析・モデル化の考え方

Bernoulli標本の最尤推定量は $\hat p=\bar X$ で、

\sqrt{n}(\hat p-p)\xrightarrow{d}N\{0,p(1-p)\}

です。したがって $\hat p$ は近似的に $N\{p,p(1-p)/n\}$ に従います。ただし、低不良率かつ小標本では不良0件が頻発し、対称なWald区間は不適切になり得ます。

Pythonで確認する

an_rng = np.random.default_rng(SEED + 80)
p, reps = 0.03, 50_000
rows = []
fig, axes = plt.subplots(1, 3, figsize=(13, 3.8))
for ax, n in zip(axes, [50, 500, 5_000]):
    phat = an_rng.binomial(n, p, reps) / n
    z = np.sqrt(n) * (phat - p) / np.sqrt(p * (1 - p))
    se_hat = np.sqrt(phat * (1 - phat) / n)
    lower, upper = phat - 1.96 * se_hat, phat + 1.96 * se_hat
    coverage = np.mean((lower <= p) & (p <= upper))
    rows.append([n, np.mean(phat == 0), z.mean(), z.std(ddof=1), coverage])
    ax.hist(z, bins=45, density=True, alpha=0.65)
    xs = np.linspace(-4, 4, 300)
    ax.plot(xs, norm.pdf(xs), color="tab:red")
    ax.set_title(f"標準化した不良率（n={n}）")
    ax.set_xlabel("標準化推定誤差")
    ax.set_ylabel("密度")
    ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

asymptotic_table = pd.DataFrame(rows, columns=["検査数_n", "不良0件率", "標準化誤差の平均", "標準化誤差のSD", "Wald95%区間被覆率"])
display(asymptotic_table.round(4))

png

	検査数_n	不良0件率	標準化誤差の平均	標準化誤差のSD	Wald95%区間被覆率
0	50	0.2181	-0.0040	0.9973	0.7812
1	500	0.0000	0.0046	1.0025	0.9233
2	5000	0.0000	0.0038	0.9972	0.9462

結果の読み取り

$n=50$ では不良0件が多く分布は離散的で、推定標準誤差も0となるためWald区間の被覆率が不足します。 $n$ が増えると標準化誤差は標準正規分布へ近づき、被覆率も名目95%へ近づきます。低不良率・小標本ではWilson区間や正確法を使い、「大標本だから正規近似」と判断する前に期待不良件数 $np$ も確認します。

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

No.071〜No.080から、製造KPIの判断には次の4層が必要です。

安定性を確認する：大数の法則に期待するだけでなく、累積平均、層別、変更点を可視化する
推定誤差を数量化する：中心極限定理、スラツキー、デルタ法で、点推定を標準誤差と意思決定KPIの幅へ変換する
保証と予測を分ける：集中不等式の上界は実確率ではない。仮定が少ないほど保守性が増える
有限標本を確認する：漸近理論を使う前に、歪み、裾、低発生率、自己相関、有効標本数を検証する

標本数は「多い・少ない」ではなく、許容誤差、誤判断コスト、工程構造に対して十分かで決めます。

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

1. 観測単位と分母を固定する

個体、ロット、日、設備のどれを1観測とするか、不良率の分母、再検査品の扱い、欠測・測定誤判定をデータ辞書へ記録します。

2. 独立同分布の仮定を点検する

時系列プロットと自己相関、設備・製品・材料・勤務帯による層別、工程変更履歴を確認します。連続ロットが相関する場合は、有効標本数やブロック単位の評価を検討します。

3. 許容誤差と誤判断コストを合意する

工程停止の誤警報、異常見逃し、追加検査、納期遅延を金額・工数へ換算し、必要な信頼度と標本数を決めます。

4. 上界・近似・実績を分けて運用する

帳票には、理論上界、モデルによる推定確率、過去実績の超過率を別項目で表示し、用いた仮定、対象期間、再計算条件を残します。導入後も被覆率や警報頻度を監視し、閾値を更新します。

まとめ

大数の法則は累積平均の安定を支えるが、工程変化や依存性までは保証しない
中心極限定理により、歪んだ個票でも大標本の平均は正規近似できる
スラツキーの定理とデルタ法は、未知分散の置換とKPI変換後の誤差評価を支える
チェビシェフ、マルコフ、Chernoff、Hoeffdingなどは、異なる仮定の下で確率を上から抑える
集中不等式の値は予測確率ではなく保証上界であり、保守性もコストになる
漸近正規性は有用だが、低不良率・小標本では有限標本向けの方法を選ぶ

極限定理は「データが多ければ大丈夫」という標語ではありません。何を仮定し、どの誤差を許し、どの意思決定を守るのかを明文化して初めて、品質保証と生産管理の道具になります。

法人向けのご相談

数理工房では、製造業の品質管理、生産KPI、抜取検査、工程能力、異常検知について、データ定義の整理から統計モデル、シミュレーション、Python notebookによる社内研修、運用設計まで支援しています。

「検査数の根拠を説明したい」「平均KPIの誤差幅を経営会議で示したい」「既存の管理限界や警報閾値を統計的に再設計したい」といった段階からご相談いただけます。

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

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

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

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.071：大数の法則

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.072：中心極限定理

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.073：スラツキーの定理

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.074：デルタ法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.075：チェビシェフの不等式

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.076：マルコフの不等式

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.077：Chernoff境界

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.078：Hoeffding不等式

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.079：Concentration Inequality（集中不等式）

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.080：漸近正規性

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

1. 観測単位と分母を固定する

2. 独立同分布の仮定を点検する

3. 許容誤差と誤判断コストを合意する

4. 上界・近似・実績を分けて運用する

まとめ

法人向けのご相談

Related Posts

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

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

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

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

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