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

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

2020-04-16

抜取検査から工程推定へ：限られたデータで製造条件を判断する

全数を測れない精密部品工程を題材に、標本から母集団を推定する考え方を、推定量の性質、最尤法、情報量、推定精度の限界、混合工程の分離まで一貫して扱います。数式と Python のシミュレーションを往復しながら、「何点測ればよいか」「その推定値をどこまで信頼できるか」「混ざったデータをどう解釈するか」を製造業の意思決定へ接続します。

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

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

量産工程では、製造した全数について寸法を精密測定することは、測定時間・設備能力・費用の面から現実的でない場合があります。そのため、各ロットから一部を抜き取り、平均寸法、ばらつき、不適合率などを推定して、出荷、条件補正、追加検査を判断します。

しかし、標本の数値は抜き取るたびに変わります。標本平均が規格中心に近いというだけで工程全体が安定しているとは限らず、標本分散の計算方法やサンプルサイズによって推定の性質も変わります。本記事では、推定値そのものだけでなく、偏り・ばらつき・情報量・理論上の精度限界を判断材料として扱います。

現場でよくある状況

1ロット10点の慣例が続いているが、必要な推定精度から点数を決めていない
月報の平均値は確認しているが、母集団と標本の範囲・期間・抽出方法が曖昧
標本分散で $n$ と $n-1$ のどちらを使うかが担当者やツールで異なる
推定値の小数点以下を細かく報告する一方、標準誤差を併記していない
複数設備や複数条件の測定値が混在し、単一分布として平均・分散を計算している

こうした状態では、計算が正しくても意思決定を誤ります。まず推定対象を定義し、抽出設計と推定方法を対応させる必要があります。

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

母集団の真の平均や分散は通常観測できません。観測できるのは、抽出された有限個のデータだけです。同じ工程から抜き取っても標本平均は毎回異なり、サンプルサイズが小さいほど大きく振れます。また、設備、シフト、材料ロットなどの構成比が偏れば、ランダム誤差だけでなく系統的な偏りが入ります。

したがって「推定値がいくつか」だけでは不十分です。対象母集団、抽出方法、推定量、標準誤差、モデル仮定を一組で管理し、推定の不確実性が出荷損失や調整費用に対して十分小さいかを判断します。

今回扱うノックの全体像

No.	テーマ	製造業での問い	主な確認内容
081	標本と母集団	抜取データは工程全体を表しているか	母集団・標本・抽出枠
082	標本平均	平均値は何点測れば安定するか	標本分布・標準誤差
083	標本分散	工程ばらつきをどう計算するか	自由度・ $n-1$ 補正
084	不偏推定量	繰り返したとき推定に偏りはないか	バイアスの比較
085	一致推定量	データを増やせば真値へ近づくか	MSE とサンプルサイズ
086	十分統計量	推定に必要な情報を要約できるか	二項モデルと不適合数
087	最尤推定法	観測データを最も説明する条件は何か	尤度・数値最適化
088	フィッシャー情報量	測定数と測定精度は推定力にどう効くか	情報量・標準誤差
089	クラメール・ラオ下界	不偏推定の精度にはどんな限界があるか	下界・効率性
090	EMアルゴリズム	混在した工程状態をどう分離するか	潜在変数・反復推定

Python 環境の準備

外部データには依存せず、NumPy、pandas、SciPy、matplotlib だけを使います。乱数生成器の seed を固定し、同じ環境で表とグラフを再現できるようにします。グラフ内は Markdown 変換先でのフォント差を避けるため英語表記にします。

%matplotlib inline
import platform

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
from IPython.display import display
from scipy import optimize, stats

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

print(f"Python      : {platform.python_version()}")
print(f"NumPy       : {np.__version__}")
print(f"pandas      : {pd.__version__}")
print(f"SciPy       : {scipy.__version__}")

Python      : 3.13.1
NumPy       : 2.5.1
pandas      : 3.0.3
SciPy       : 1.18.0

架空データの作成

精密シャフトを2台の設備（M1、M2）、3シフトで生産する1か月分24,000本の架空データを作ります。管理対象は直径で、目標値は10.000 mm、規格は 9.980〜10.020 mm とします。設備・シフトごとに小さな平均差を持たせ、全数データから層別抜取標本も作ります。

以降、全数データを説明のための「既知の有限母集団」として使います。実務では真の母集団値は分からないため、標本から推定します。また No.090 用に、設備ラベルが欠落した混合測定データも別に生成します。

n_population = 24_000
machines = rng.choice(["M1", "M2"], size=n_population, p=[0.58, 0.42])
shifts = rng.choice(["Day", "Evening", "Night"], size=n_population, p=[0.45, 0.35, 0.20])
machine_offset = pd.Series(machines).map({"M1": -0.0015, "M2": 0.0025}).to_numpy()
shift_offset = pd.Series(shifts).map({"Day": 0.0000, "Evening": 0.0008, "Night": -0.0007}).to_numpy()
sigma = np.where(machines == "M1", 0.0060, 0.0075)
diameter = 10.000 + machine_offset + shift_offset + rng.normal(0, sigma)

population_df = pd.DataFrame({
    "unit_id": np.arange(1, n_population + 1),
    "machine": machines,
    "shift": shifts,
    "diameter_mm": diameter,
})
population_df["nonconforming"] = ~population_df["diameter_mm"].between(9.980, 10.020)

# 各設備×シフトから80本ずつ抽出した層別標本
sample_df = (
    population_df.groupby(["machine", "shift"], group_keys=False)
    .sample(n=80, random_state=SEED)
    .sort_values("unit_id")
    .reset_index(drop=True)
)

# ラベルが欠落した2状態混合データ（No.090）
n_mixed = 700
hidden_state = rng.choice([0, 1], n_mixed, p=[0.64, 0.36])
mixed_values = rng.normal(
    np.where(hidden_state == 0, 9.990, 10.014),
    np.where(hidden_state == 0, 0.0050, 0.0065),
)

display(population_df.head())
display(pd.DataFrame({
    "dataset": ["有限母集団", "層別抜取標本", "ラベルなし混合データ"],
    "rows": [len(population_df), len(sample_df), len(mixed_values)],
}))

	unit_id	machine	shift	diameter_mm	nonconforming
0	1	M1	Night	9.998475	False
1	2	M2	Day	10.015734	False
2	3	M1	Night	9.993877	False
3	4	M2	Night	10.001658	False
4	5	M2	Day	9.999186	False

	dataset	rows
0	有限母集団	24000
1	層別抜取標本	480
2	ラベルなし混合データ	700

No.081：標本と母集団 — 抜取検査の対象を定義する

実務での意味

母集団は「推定したい対象の全体」、標本は「実際に観測した一部」です。たとえば今月の全製品を母集団とするのか、今後も同じ条件で作る製品まで含む概念的な母集団とするのかで、推定結果の適用範囲が変わります。抽出枠から夜勤や特定設備が漏れれば、測定点数が多くても代表性は得られません。

分析・モデル化の考え方

有限母集団の平均を $mu_N=N^{-1}\sum_{i=1}^{N}x_i$ 、標本平均を $\bar{x}=n^{-1}\sum_{i=1}^{n}x_i$ とします。単純無作為抽出では各個体に等しい抽出機会があります。層別抽出では設備・シフトなどの層ごとに抽出し、母集団構成比と標本構成比が違う場合は重み付けします。ここでは全数データと層別標本の分布・構成を比較します。

Pythonで確認する

population_summary = pd.Series({
    "件数": len(population_df),
    "平均直径_mm": population_df["diameter_mm"].mean(),
    "標準偏差_mm": population_df["diameter_mm"].std(ddof=0),
    "不適合率": population_df["nonconforming"].mean(),
})
sample_summary = pd.Series({
    "件数": len(sample_df),
    "平均直径_mm": sample_df["diameter_mm"].mean(),
    "標準偏差_mm": sample_df["diameter_mm"].std(ddof=1),
    "不適合率": sample_df["nonconforming"].mean(),
})
comparison_081 = pd.concat(
    [population_summary.rename("有限母集団"), sample_summary.rename("層別標本")], axis=1
)
display(comparison_081.round(6))

fig, ax = plt.subplots()
ax.hist(population_df["diameter_mm"], bins=45, density=True, alpha=0.45, label="Population")
ax.hist(sample_df["diameter_mm"], bins=25, density=True, alpha=0.55, label="Stratified sample")
ax.set_title("Population and Inspection Sample")
ax.set_xlabel("Diameter (mm)")
ax.set_ylabel("Density")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	有限母集団	層別標本
件数	24000.000000	480.00000
平均直径_mm	10.000264	10.00063
標準偏差_mm	0.007024	0.00690
不適合率	0.006000	0.00625

png

結果の読み取り

層別標本の分布は母集団の中心と広がりを概ね再現しますが、値は完全には一致しません。この差が標本誤差です。今回の抽出は各層80本の均等割付なので、母集団の層構成比とは異なります。工程全体の指標を厳密に推定する際は母集団構成比で重み付けし、設備別比較が目的なら層別値をそのまま報告するなど、目的に合わせます。対象期間、対象設備、除外条件、抽出方法を分析結果と一緒に記録することが重要です。

No.082：標本平均 — 測定点数と平均値の安定性を結び付ける

実務での意味

標本平均は工程中心を表す基本 KPI ですが、少数点の平均は偶然に左右されます。測定点数を増やすと平均値は安定する一方、検査費用と判定時間が増えます。標本平均の標準誤差を使えば、必要精度と測定負荷のトレードオフを定量化できます。

分析・モデル化の考え方

独立同分布の観測 $X_1,\ldots,X_n$ に対し、標本平均は

\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_i,\qquad E[\bar{X}]=\mu,\qquad \mathrm{Var}(\bar{X})=\frac{\sigma^2}{n}

です。標準誤差は $\sigma/\sqrt{n}$ で減少するため、精度を半分にするには原則4倍のデータが必要です。母集団からサイズの異なる標本を繰り返し抽出し、標本平均の分布を確認します。

Pythonで確認する

values = population_df["diameter_mm"].to_numpy()
sample_sizes = [10, 40, 160]
repetitions = 2_000
mean_draws = {
    n: np.array([rng.choice(values, n, replace=False).mean() for _ in range(repetitions)])
    for n in sample_sizes
}
summary_082 = pd.DataFrame([
    {
        "n": n,
        "標本平均の平均": draws.mean(),
        "標本平均の実測SD": draws.std(ddof=1),
        "理論SE近似": values.std(ddof=0) / np.sqrt(n),
    }
    for n, draws in mean_draws.items()
])
display(summary_082.round(7))

fig, ax = plt.subplots()
ax.boxplot([mean_draws[n] for n in sample_sizes], tick_labels=[str(n) for n in sample_sizes])
ax.axhline(values.mean(), color="red", linestyle="--", label="Population mean")
ax.set_title("Sampling Distribution of the Mean")
ax.set_xlabel("Sample size n")
ax.set_ylabel("Sample mean diameter (mm)")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	n	標本平均の平均	標本平均の実測SD	理論SE近似
0	10	10.000246	0.002202	0.002221
1	40	10.000261	0.001125	0.001111
2	160	10.000268	0.000566	0.000555

png

結果の読み取り

$n$ が大きいほど標本平均の箱が狭まり、実測した標準偏差も理論的な標準誤差に近づきます。一方、測定数を4倍にしても標準誤差は約半分にしかなりません。検査点数は慣例ではなく、許容する推定誤差、工程標準偏差、測定費用、判断の損失を使って設計します。連続生産で自己相関がある場合は独立性が崩れ、単純式より実効的な情報量が小さくなる点にも注意が必要です。

No.083：標本分散 — 工程ばらつきを自由度とともに推定する

実務での意味

平均が規格中心でも、ばらつきが大きければ規格外が増えます。工程能力や管理限界の基礎となる分散を過小評価すると、品質リスクを楽観視します。特に小標本では、分母を $n$ とするか $n-1$ とするかの差が無視できません。

分析・モデル化の考え方

母平均も標本から推定するため、偏差 $X_i-\bar{X}$ には独立な情報が $n-1$ 個しかありません。不偏標本分散は

S^2=\frac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar{X})^2

と定義され、独立同分布で分散が有限なら $E[S^2]=\sigma^2$ です。分母 $n$ の分散は正規分布の最尤推定量ですが、有限標本では平均的に小さくなります。

Pythonで確認する

true_variance = values.var(ddof=0)
variance_rows = []
for n in [5, 10, 30, 100]:
    draws = rng.choice(values, size=(4_000, n), replace=True)
    variance_rows.append({
        "n": n,
        "分母nの平均": draws.var(axis=1, ddof=0).mean(),
        "分母n-1の平均": draws.var(axis=1, ddof=1).mean(),
        "母分散": true_variance,
    })
summary_083 = pd.DataFrame(variance_rows)
display(summary_083.round(9))

fig, ax = plt.subplots()
ax.plot(summary_083["n"], summary_083["分母nの平均"], "o-", label="Denominator n")
ax.plot(summary_083["n"], summary_083["分母n-1の平均"], "s-", label="Denominator n-1")
ax.axhline(true_variance, color="black", linestyle="--", label="Population variance")
ax.set_title("Average Variance Estimates")
ax.set_xlabel("Sample size n")
ax.set_ylabel("Estimated variance (mm squared)")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	n	分母nの平均	分母n-1の平均	母分散
0	5	0.000041	0.000051	0.000049
1	10	0.000044	0.000049	0.000049
2	30	0.000048	0.000050	0.000049
3	100	0.000049	0.000049	0.000049

png

結果の読み取り

小標本では分母 $n$ の推定値が母分散を系統的に下回り、 $n-1$ 補正がその偏りを取り除いています。ただし、不偏であることは一回の推定が真値に近いことを保証しません。 $n=5$ の標本分散は繰り返し間の変動が大きいため、ばらつき監視では測定数、合理的な群分け、時系列変化も合わせて確認します。利用ライブラリの ddof 設定をデータ仕様に明記すると、部署間の計算差を防げます。

No.084：不偏推定量 — 長期的な推定の偏りを点検する

実務での意味

推定ルールが毎月わずかに工程ばらつきを過小評価するなら、長期的には検査基準や能力評価が楽観側へずれます。不偏性は、同じ条件で標本抽出を繰り返したとき、推定量の平均が真の値に一致するという性質です。

分析・モデル化の考え方

パラメータ $\theta$ の推定量 $\hat{\theta}$ が

$E[\hat{\theta}]=\theta$

を満たすとき不偏といいます。バイアスは $\mathrm{Bias}(\hat{\theta})=E[\hat{\theta}]-\theta$ です。標本平均と $n-1$ で割る標本分散は不偏ですが、分母 $n$ の分散推定量は $-\sigma^2/n$ のバイアスを持ちます。シミュレーションで3つの推定ルールを比較します。

Pythonで確認する

n = 12
draws = rng.choice(values, size=(8_000, n), replace=True)
estimators = pd.DataFrame({
    "推定対象": ["母平均", "母分散", "母分散"],
    "推定量": ["標本平均", "分母n-1の標本分散", "分母nの標本分散"],
    "真値": [values.mean(), true_variance, true_variance],
    "推定量の期待値近似": [
        draws.mean(axis=1).mean(),
        draws.var(axis=1, ddof=1).mean(),
        draws.var(axis=1, ddof=0).mean(),
    ],
})
estimators["バイアス"] = estimators["推定量の期待値近似"] - estimators["真値"]
display(estimators.round(9))

variance_estimates = [draws.var(axis=1, ddof=1), draws.var(axis=1, ddof=0)]
fig, ax = plt.subplots()
ax.boxplot(variance_estimates, tick_labels=["n-1", "n"], showfliers=False)
ax.axhline(true_variance, color="red", linestyle="--", label="Population variance")
ax.set_title("Bias of Variance Estimators")
ax.set_xlabel("Variance denominator")
ax.set_ylabel("Estimated variance (mm squared)")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	推定対象	推定量	真値	推定量の期待値近似	バイアス
0	母平均	標本平均	10.000264	10.000287	2.290500e-05
1	母分散	分母n-1の標本分散	0.000049	0.000049	-7.500000e-08
2	母分散	分母nの標本分散	0.000049	0.000045	-4.181000e-06

png

結果の読み取り

標本平均と $n-1$ の標本分散はシミュレーション誤差の範囲で真値に一致し、分母 $n$ の推定値には理論どおり負のバイアスがあります。ただし実務では、不偏性だけで推定量を選びません。多少のバイアスを許すことで推定値の振れを大きく減らせる場合もあるため、平均二乗誤差、異常時の損失、説明可能性も比較します。また抽出自体に偏りがあれば、数式上不偏な推定量でも対象母集団には不偏になりません。

No.085：一致推定量 — データ増加が真値への接近につながるか確認する

実務での意味

センサーや自動測定機を導入してデータ量を増やしても、推定方法が適切でなければ真値へ近づくとは限りません。一致性は、サンプルサイズを大きくしたとき推定量が真のパラメータへ確率的に近づく性質です。長期データ蓄積の価値を支える条件でもあります。

分析・モデル化の考え方

任意の $\varepsilon>0$ に対し

P(|\hat{\theta}_n-\theta|>\varepsilon)\to 0\quad(n\to\infty)

なら $\hat{\theta}_n$ は一致推定量です。標本平均では大数の法則により一致性が得られ、平均二乗誤差は独立同分布なら概ね $\sigma^2/n$ で減少します。データ量別に誤差と許容誤差超過率を測ります。

Pythonで確認する

mu = values.mean()
epsilon = 0.001  # 1 micrometer
consistency_rows = []
for n in [5, 10, 25, 50, 100, 250, 500]:
    sample_means = rng.choice(values, size=(3_000, n), replace=True).mean(axis=1)
    consistency_rows.append({
        "n": n,
        "MSE": np.mean((sample_means - mu) ** 2),
        "絶対誤差1μm超の割合": np.mean(np.abs(sample_means - mu) > epsilon),
    })
summary_085 = pd.DataFrame(consistency_rows)
display(summary_085.round(8))

fig, ax = plt.subplots()
ax.loglog(summary_085["n"], summary_085["MSE"], "o-", label="Simulated MSE")
ax.loglog(summary_085["n"], true_variance / summary_085["n"], "--", label="Variance / n")
ax.set_title("Consistency of the Sample Mean")
ax.set_xlabel("Sample size n (log scale)")
ax.set_ylabel("Mean squared error (log scale)")
ax.grid(True, which="both", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	n	MSE	絶対誤差1μm超の割合
0	5	9.720000e-06	0.757000
1	10	4.840000e-06	0.655667
2	25	1.970000e-06	0.468333
3	50	9.400000e-07	0.301000
4	100	5.100000e-07	0.157667
5	250	2.000000e-07	0.023333
6	500	1.000000e-07	0.001000

png

結果の読み取り

サンプルサイズの増加とともに MSE と1 μm超の誤差割合が下がり、標本平均の一致性を確認できます。ただし、同じ設備状態から高頻度に測った相関データ、校正ずれを含む測定値、将来の工程変更に対しては、件数を増やすだけでは誤差が消えません。一致性はモデル仮定の下での性質です。設備・材料・期間をまたいだ代表性、測定系のバイアス、ドリフトを別途管理します。

No.086：十分統計量 — 不適合数に推定情報を集約する

実務での意味

検査結果の生データをすべて会議資料へ並べる必要はありません。推定対象に関する情報を失わず要約できれば、監視、通信、保存、説明を簡潔にできます。ただし「十分」は特定の確率モデルとパラメータに対する性質であり、原因分析に必要な時系列順序や設備情報まで不要になるわけではありません。

分析・モデル化の考え方

各検査が独立に確率 $p$ で不適合となるベルヌーイモデルでは、 $n$ 件中の不適合総数 $T=\sum_i X_i$ の尤度は

L(p\mid\boldsymbol{x})=p^T(1-p)^{n-T}

です。尤度は生の並び順ではなく $(n,T)$ だけに依存するため、因子分解定理により $T$ は $p$ の十分統計量です。同じ件数・同じ不適合数で並び順だけ異なる2系列の尤度を比較します。

Pythonで確認する

sequence_a = np.array([1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
sequence_b = np.array([0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0])
p_grid = np.linspace(0.001, 0.45, 400)

def bernoulli_likelihood(sequence, p):
    total = sequence.sum()
    return p ** total * (1 - p) ** (len(sequence) - total)

lik_a = bernoulli_likelihood(sequence_a, p_grid)
lik_b = bernoulli_likelihood(sequence_b, p_grid)
display(pd.DataFrame({
    "系列": ["A", "B"],
    "検査数n": [len(sequence_a), len(sequence_b)],
    "不適合数T": [sequence_a.sum(), sequence_b.sum()],
    "最尤不適合率T/n": [sequence_a.mean(), sequence_b.mean()],
    "尤度曲線の最大差": [np.max(np.abs(lik_a - lik_b)), np.max(np.abs(lik_a - lik_b))],
}))

fig, ax = plt.subplots()
ax.plot(p_grid, lik_a, label="Sequence A")
ax.plot(p_grid, lik_b, "--", label="Sequence B")
ax.axvline(sequence_a.mean(), color="red", linestyle=":", label="T / n")
ax.set_title("Likelihood Depends on the Defect Count")
ax.set_xlabel("Nonconforming probability p")
ax.set_ylabel("Likelihood")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	系列	検査数n	不適合数T	最尤不適合率T/n	尤度曲線の最大差
0	A	20	3	0.15	0.0
1	B	20	3	0.15	0.0

png

結果の読み取り

2系列は並び順が異なっても $(n,T)=(20,3)$ が同じなので、 $p$ に関する尤度曲線は完全に重なります。不適合率の推定だけなら検査数と不適合数へ情報を集約できます。一方、3件が連続発生したか、特定設備・材料に集中したかは原因究明に重要です。集計 KPI 用の十分統計量と、トレーサビリティ用の生データを目的別に保持し、モデルの独立・同一確率という仮定も監視します。

No.087：最尤推定法 — 観測寸法を最も説明する工程条件を求める

実務での意味

工程平均や標準偏差をデータから一貫した規則で求めることは、能力評価、条件補正、異常検知の基礎です。最尤推定法は、観測したデータが最も起こりやすくなるパラメータを選びます。多くの統計モデルや機械学習モデルに共通する推定原理です。

分析・モデル化の考え方

寸法 $x_i$ が独立に正規分布 $N(\mu,\sigma^2)$ に従うとすると、対数尤度は

\ell(\mu,\sigma)=-n\log\sigma-\frac{n}{2}\log(2\pi) -\frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\mu)^2

です。最大化すると $\hat{\mu}=\bar{x}$ 、 $\hat{\sigma}^2=n^{-1}\sum_i(x_i-\bar{x})^2$ が得られます。解析解と数値最適化の結果を照合します。

Pythonで確認する

m1_data = sample_df.loc[sample_df["machine"] == "M1", "diameter_mm"].to_numpy()

def normal_negative_loglik(params, x):
    mu_value, log_sigma = params
    sigma_value = np.exp(log_sigma)
    return -np.sum(stats.norm.logpdf(x, loc=mu_value, scale=sigma_value))

result = optimize.minimize(
    normal_negative_loglik,
    x0=np.array([m1_data.mean(), np.log(m1_data.std(ddof=1))]),
    args=(m1_data,),
    method="BFGS",
)
mu_mle, sigma_mle = result.x[0], np.exp(result.x[1])
display(pd.DataFrame({
    "method": ["Analytic MLE", "Numerical MLE"],
    "mu": [m1_data.mean(), mu_mle],
    "sigma": [m1_data.std(ddof=0), sigma_mle],
    "converged": [True, result.success],
}).round(8))

mu_grid = np.linspace(mu_mle - 0.0025, mu_mle + 0.0025, 300)
relative_loglik = np.array([
    -normal_negative_loglik((mu_value, np.log(sigma_mle)), m1_data) for mu_value in mu_grid
])
relative_loglik -= relative_loglik.max()
fig, ax = plt.subplots()
ax.plot(mu_grid, relative_loglik)
ax.axvline(mu_mle, color="red", linestyle="--", label="MLE")
ax.set_title("Profile of the Normal Log-Likelihood")
ax.set_xlabel("Process mean mu (mm)")
ax.set_ylabel("Relative log-likelihood")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	method	mu	sigma	converged
0	Analytic MLE	9.99879	0.005808	True
1	Numerical MLE	9.99879	0.005808	True

png

結果の読み取り

解析解と数値解が一致し、対数尤度は推定平均で最大になります。最尤推定は複雑なモデルへ拡張しやすい一方、正規性、独立性、単一工程という仮定が誤っていれば、得られる値の意味も変わります。また分散の最尤推定量は分母 $n$ で、No.083 の不偏分散とは目的と性質が異なります。推定方法名、尤度、収束判定、初期値、適合度診断を記録し、最適化成功とモデル妥当性を区別します。

No.088：フィッシャー情報量 — 測定数と精度から推定力を評価する

実務での意味

同じ100点でも、高精度測定器で安定工程を測る場合と、測定誤差が大きい環境で変動工程を測る場合では、工程平均について得られる情報が異なります。フィッシャー情報量は、尤度が未知パラメータの周囲でどれだけ鋭く変化するかを表し、検査数・測定精度・推定精度を結び付けます。

分析・モデル化の考え方

スコア $\partial\ell/\partial\theta$ の二乗の期待値をフィッシャー情報量

I(\theta)=E\left[\left(\frac{\partial}{\partial\theta}\log f(X;\theta)\right)^2\right]

と定義します。分散 $\sigma^2$ が既知の正規分布で平均 $\mu$ を推定する場合、1観測の情報量は $1/\sigma^2$ 、 $n$ 観測では $I_n(\mu)=n/\sigma^2$ です。情報量の逆平方根 $1/\sqrt{I_n}=\sigma/\sqrt{n}$ が推定精度の尺度になります。

Pythonで確認する

design_rows = []
for measurement_system, sigma_value in [("High precision", 0.006), ("Standard", 0.009)]:
    for n in [10, 25, 50, 100, 200]:
        information = n / sigma_value**2
        design_rows.append({
            "measurement_system": measurement_system,
            "n": n,
            "sigma_mm": sigma_value,
            "Fisher_information": information,
            "expected_SE_mm": 1 / np.sqrt(information),
        })
information_df = pd.DataFrame(design_rows)
display(information_df.round(7))

fig, ax = plt.subplots()
for name, group in information_df.groupby("measurement_system"):
    ax.plot(group["n"], group["expected_SE_mm"] * 1000, "o-", label=name)
ax.set_title("Fisher Information and Expected Precision")
ax.set_xlabel("Number of measurements n")
ax.set_ylabel("Expected standard error (micrometers)")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	measurement_system	n	sigma_mm	Fisher_information	expected_SE_mm
0	High precision	10	0.006	2.777778e+05	0.001897
1	High precision	25	0.006	6.944444e+05	0.001200
2	High precision	50	0.006	1.388889e+06	0.000848
3	High precision	100	0.006	2.777778e+06	0.000600
4	High precision	200	0.006	5.555556e+06	0.000424
5	Standard	10	0.009	1.234568e+05	0.002846
6	Standard	25	0.009	3.086420e+05	0.001800
7	Standard	50	0.009	6.172840e+05	0.001273
8	Standard	100	0.009	1.234568e+06	0.000900
9	Standard	200	0.009	2.469136e+06	0.000636

png

結果の読み取り

測定数を増やすほど期待標準誤差は下がりますが、平方根則のため限界効果は逓減します。また標準偏差の小さい高精度系は、同じ測定数でも多くの情報を持ちます。これは検査点数の追加と測定系改善を同じ精度指標で比較できることを意味します。ただし、式の $\sigma$ には工程変動と測定誤差が混ざり得ます。Gage R&R などで測定系変動を分離し、費用、サイクルタイム、必要精度を含めて検査設計を選びます。

No.089：クラメール・ラオ下界 — 不偏推定の精度限界を知る

実務での意味

推定アルゴリズムを複雑にすれば無制限に精度が上がるわけではありません。一定のモデルとデータ量の下では、不偏推定量の分散に理論的な下限があります。下限を知ると、手法改善で得られる余地と、測定数・測定精度を改善すべき領域を分けられます。

分析・モデル化の考え方

正則条件の下で、不偏推定量 $\hat{\theta}$ にはクラメール・ラオ下界

\mathrm{Var}(\hat{\theta})\geq\frac{1}{I_n(\theta)}

が成り立ちます。分散既知の正規分布の平均では下界は $\sigma^2/n$ で、標本平均はこの下界を達成する効率的な推定量です。標本平均と標本中央値の分散をシミュレーションで比較します。

Pythonで確認する

sigma_known = 0.007
crlb_rows = []
for n in [10, 30, 100, 300]:
    normal_draws = rng.normal(10.0, sigma_known, size=(8_000, n))
    mean_estimator = normal_draws.mean(axis=1)
    median_estimator = np.median(normal_draws, axis=1)
    lower_bound = sigma_known**2 / n
    crlb_rows.append({
        "n": n,
        "CRLB": lower_bound,
        "標本平均の分散": mean_estimator.var(ddof=1),
        "標本中央値の分散": median_estimator.var(ddof=1),
        "平均の効率_CRLB/分散": lower_bound / mean_estimator.var(ddof=1),
    })
crlb_df = pd.DataFrame(crlb_rows)
display(crlb_df.round(10))

fig, ax = plt.subplots()
ax.loglog(crlb_df["n"], crlb_df["CRLB"], "k--", label="Cramer-Rao lower bound")
ax.loglog(crlb_df["n"], crlb_df["標本平均の分散"], "o-", label="Sample mean")
ax.loglog(crlb_df["n"], crlb_df["標本中央値の分散"], "s-", label="Sample median")
ax.set_title("Estimator Variance and the Cramer-Rao Bound")
ax.set_xlabel("Sample size n (log scale)")
ax.set_ylabel("Estimator variance (log scale)")
ax.grid(True, which="both", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	n	CRLB	標本平均の分散	標本中央値の分散	平均の効率_CRLB/分散
0	10	4.900000e-06	4.922100e-06	6.895100e-06	0.995514
1	30	1.633300e-06	1.673900e-06	2.490800e-06	0.975772
2	100	4.900000e-07	4.748000e-07	7.492000e-07	1.031930
3	300	1.633000e-07	1.656000e-07	2.566000e-07	0.986608

png

結果の読み取り

標本平均の分散は下界にほぼ一致し、正規モデルの平均推定では効率的です。標本中央値の分散は大きく、正規分布という仮定の下では効率が低くなります。一方、外れ値や裾の重い分布では中央値の頑健性が意思決定上有利な場合があります。クラメール・ラオ下界は不偏性、モデルの正則性、パラメータ設定などの条件に依存するため、あらゆる推定量への万能な順位付けではありません。現場分布と損失関数に照らして使います。

No.090：EMアルゴリズム — ラベルのない混合工程を反復推定する

実務での意味

設備IDや段取り条件の記録が欠落すると、異なる工程状態の測定値が一つの分布に混ざります。全体平均だけでは、中心の異なる2状態や状態別のばらつきを見落とします。EMアルゴリズムは、どの状態から生成されたかという潜在ラベルを確率的に補いながら、各状態のパラメータを推定します。

分析・モデル化の考え方

2成分正規混合モデルを

p(x)=\sum_{k=1}^{2}\pi_k\,\mathcal{N}(x\mid\mu_k,\sigma_k^2)

とします。Eステップでは現在のパラメータから各観測が成分 $k$ に属する事後確率（責任度） $r_{ik}$ を計算します。Mステップでは責任度を重みとして $\pi_k,\mu_k,\sigma_k^2$ を更新します。観測対数尤度の増加が十分小さくなるまで反復します。

Pythonで確認する

def gaussian_mixture_em(x, n_components=2, max_iter=200, tol=1e-9):
    means = np.quantile(x, [0.25, 0.75]).astype(float)
    stds = np.full(n_components, x.std(ddof=1))
    weights = np.full(n_components, 1 / n_components)
    loglik_history = []

    for _ in range(max_iter):
        weighted_density = np.column_stack([
            weights[k] * stats.norm.pdf(x, means[k], stds[k])
            for k in range(n_components)
        ])
        denominator = weighted_density.sum(axis=1, keepdims=True)
        responsibilities = weighted_density / denominator

        effective_n = responsibilities.sum(axis=0)
        weights = effective_n / len(x)
        means = (responsibilities * x[:, None]).sum(axis=0) / effective_n
        variances = (
            responsibilities * (x[:, None] - means) ** 2
        ).sum(axis=0) / effective_n
        stds = np.sqrt(np.maximum(variances, 1e-12))

        loglik = np.log(denominator[:, 0]).sum()
        loglik_history.append(loglik)
        if len(loglik_history) > 1 and abs(loglik_history[-1] - loglik_history[-2]) < tol:
            break

    order = np.argsort(means)
    return weights[order], means[order], stds[order], responsibilities[:, order], loglik_history

em_weights, em_means, em_stds, responsibilities, loglik_history = gaussian_mixture_em(mixed_values)
em_result = pd.DataFrame({
    "component": ["Low-center state", "High-center state"],
    "weight": em_weights,
    "mean_mm": em_means,
    "std_mm": em_stds,
})
display(em_result.round(6))
print(f"iterations: {len(loglik_history)}")

x_grid = np.linspace(mixed_values.min() - 0.004, mixed_values.max() + 0.004, 500)
mixture_density = sum(
    em_weights[k] * stats.norm.pdf(x_grid, em_means[k], em_stds[k]) for k in range(2)
)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].hist(mixed_values, bins=35, density=True, alpha=0.45, label="Observed data")
axes[0].plot(x_grid, mixture_density, color="black", label="Fitted mixture")
for k in range(2):
    axes[0].plot(
        x_grid,
        em_weights[k] * stats.norm.pdf(x_grid, em_means[k], em_stds[k]),
        "--",
        label=f"Component {k + 1}",
    )
axes[0].set_title("EM Fit to Mixed Process Data")
axes[0].set_xlabel("Diameter (mm)")
axes[0].set_ylabel("Density")
axes[0].grid(True, alpha=0.3)
axes[0].legend()

axes[1].plot(np.arange(1, len(loglik_history) + 1), loglik_history, "o-")
axes[1].set_title("Observed Log-Likelihood by Iteration")
axes[1].set_xlabel("Iteration")
axes[1].set_ylabel("Log-likelihood")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

	component	weight	mean_mm	std_mm
0	Low-center state	0.640713	9.989861	0.004939
1	High-center state	0.359287	10.013738	0.007088

iterations: 48

png

結果の読み取り

EM により、ラベルなしデータから中心・ばらつき・構成比の異なる2状態が推定され、反復に伴って観測対数尤度が収束します。高中心側の状態が規格上限リスクを持つなら、推定責任度の高い時刻を設備ログや材料履歴と突合し、原因候補を絞れます。ただし EM は局所解、初期値、成分数、ラベル入れ替わりに影響されます。複数初期値での再計算、BIC などによる成分数比較、既知ラベルでの検証を行い、欠落している設備IDを推定で恒久的に代替しないことが重要です。

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

No.081〜No.090を通じて、推定は「データから一つの数値を出す作業」ではなく、対象・抽出・情報・精度を設計する仕事だと分かります。

母集団と抽出枠を明示し、設備・シフト・期間の漏れを防ぐ
標準誤差と許容損失から測定点数を決め、平均値だけを独り歩きさせない
不偏性、一致性、効率性は異なる性質であり、目的に応じて推定量を選ぶ
十分統計量で日常監視を簡潔にしつつ、原因分析用の粒度は保持する
尤度と情報量により、推定方法・データ量・測定系の改善余地を共通尺度で議論する
混合分布が見えたら、推定結果を設備・材料・作業履歴へ戻して仮説を検証する

統計理論は現場判断を難しくするためではなく、「どの条件ならその数値を信頼できるか」を明確にするために使います。

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

推定対象の定義：対象製品、設備、期間、将来への一般化範囲を分析仕様書に記載する
抽出設計：無作為化、層別、抽出頻度、除外条件、欠測時の扱いを標準化する
測定系の保証：校正、分解能、繰返し性・再現性を確認し、工程変動と測定誤差を分ける
推定結果の表現：点推定だけでなく、標準誤差、信頼区間、サンプルサイズ、モデル仮定を併記する
モデル診断：分布形、独立性、外れ値、工程変更、混合状態を定期的に確認する
判断ルールとの接続：追加検査、条件補正、保留、出荷の責任者と発動基準を定める
再現性と監査：seed、コード版、データ期間、パラメータ、収束状況、承認履歴を保存する
段階導入：過去データ検証、並行運用、限定工程での試行を経て、誤判定コストを確認する

PoC の評価指標は推定誤差だけでなく、検査時間、追加調査件数、不適合流出、過剰調整、判断リードタイムまで含めます。

まとめ

No.081〜No.090では、標本と母集団から始め、標本平均・標本分散、不偏性・一致性、十分統計量、最尤法、フィッシャー情報量、クラメール・ラオ下界、EMアルゴリズムまでを、精密部品の抜取検査という共通課題で確認しました。

重要なのは、推定値の桁数を増やすことではありません。誰を母集団とし、どう抽出し、どの仮定で推定し、どの程度の誤差を含み、どの判断に使うかを一貫させることです。基礎理論を押さえることで、測定数追加、測定系改善、モデル高度化のどれに投資すべきかを、費用とリスクに基づいて説明できます。

法人向けのご相談

数理工房では、製造業における抜取検査設計、工程能力評価、品質データ分析、異常検知、混合工程のモデル化について、課題整理、PoC、運用設計、企業研修まで支援しています。既存の検査点数や管理基準に根拠を持たせたい場合、測定データと現場判断を接続したい場合も、利用可能なデータと意思決定から着手範囲を整理します。

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

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

抜取検査から工程推定へ：限られたデータで製造条件を判断する

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.081：標本と母集団 — 抜取検査の対象を定義する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.082：標本平均 — 測定点数と平均値の安定性を結び付ける

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.083：標本分散 — 工程ばらつきを自由度とともに推定する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.084：不偏推定量 — 長期的な推定の偏りを点検する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.085：一致推定量 — データ増加が真値への接近につながるか確認する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.086：十分統計量 — 不適合数に推定情報を集約する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.087：最尤推定法 — 観測寸法を最も説明する工程条件を求める

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.088：フィッシャー情報量 — 測定数と精度から推定力を評価する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.089：クラメール・ラオ下界 — 不偏推定の精度限界を知る

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.090：EMアルゴリズム — ラベルのない混合工程を反復推定する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

まとめ

法人向けのご相談

Related Posts

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

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

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

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

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