100本ノック / 確率統計 / 確率・統計Python100本ノック

改善効果を「偶然」で終わらせない:製造業の仮説検定10本ノック

改善効果を「偶然」で終わらせない:製造業の仮説検定10本ノック

製造現場では、設備条件の変更、治具の更新、仕入先の切替、保全後の確認など、日々さまざまな比較が行われます。しかし、標本平均や不適合件数に差が見えたというだけでは、その差が工程変更によるものか、偶然のばらつきかを区別できません。

本記事では、架空の精密部品工場を題材に、t検定からLjung–Box検定まで10種類の仮説検定をPythonで実装します。検定名を覚えることよりも、データの対応関係、分布、標本数、分散、時系列依存を確認し、効果量や信頼区間とともに意思決定へつなげることを重視します。

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

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

今回のテーマは「観測された差を、変更を採用する根拠へ変えること」です。充填量が基準からずれていないか、省エネ設定で消費電力が下がったか、ライン間の寸法に差があるか、仕入先によって強度が異なるか、といった問いを扱います。

仮説検定は、帰無仮説 H0H_0 のもとで、観測結果以上に極端なデータが得られる確率であるp値を計算します。p値が小さいことは「H0H_0 とデータの整合性が低い」ことを示しますが、改善効果の大きさや投資価値を直接示すものではありません。本記事では有意水準を原則 alpha=0.05alpha=0.05 としながら、差、効果量、信頼区間、前提条件も併記します。

現場でよくある状況

  • 改善前後の平均値だけを比べ、「下がったから効果あり」と結論づける
  • 同じ設備を変更前後で測ったデータを、独立した2群として扱う
  • 分散が大きく異なるライン比較に、等分散を仮定した検定を使う
  • 発生件数が少ない不適合に、近似が不安定な検定を使う
  • 非正規・外れ値を含むデータにも平均値ベースの検定を機械的に使う
  • 時系列データを独立な標本とみなし、連続する偏りを見落とす

検定手法の選択ミスは、効果のない変更の横展開、必要な改善の見送り、品質リスクの過小評価につながります。

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

統計的な差と実務的な差は同じではありません。標本数が大きければごく小さな差でも有意になり、標本数が小さければ重要な差を検出できないことがあります。また、複数の検定を繰り返すと、偶然の有意差を拾う確率が増えます。

実務では、検定前に次の点を定義する必要があります。

  1. 比較単位:製品、ロット、設備、日、作業者のどれか
  2. 対応関係:同じ対象を前後比較したのか、独立した群か
  3. 主要評価指標:平均、中央値、比率、分布形状、自己相関のどれか
  4. 最小重要差:現場として意味がある差はいくつか
  5. 分析計画:片側・両側、有意水準、除外規則、多重比較への対応

検定は、この設計があって初めて再現可能な意思決定になります。

今回扱うノックの全体像

No.手法架空の実務課題主に確認するもの
051t検定充填量の平均が基準値と異なるか1群の平均
052対応のあるt検定同一設備の設定変更前後で電力が変わったか対応差の平均
053Welch検定ばらつきの異なる2ラインで寸法平均が異なるか独立2群の平均
054カイ二乗検定シフトと不適合分類に関連があるかカテゴリ間の独立性
055Fisher正確検定少数事象でファームウェアと停止警報に関連があるか2×2表の比率
056ANOVA3仕入先の破断荷重平均が異なるか独立3群以上の平均
057Mann–Whitney検定2治具の非正規な段取時間に差があるか独立2群の順位
058Kruskal–Wallis検定3保全班の復旧時間分布が異なるか独立3群以上の順位
059Kolmogorov–Smirnov検定サイクル時間が事前の基準分布と一致するか分布全体
060Ljung–Box検定炉温モデル残差に時系列依存が残っていないか複数ラグの自己相関

Python 環境の準備

NumPyで架空データを生成し、pandasで表を整え、SciPyで検定を実行します。グラフはmatplotlibを使用し、日本語表示にはjapanize-matplotlibを利用します。乱数生成器は一度だけ作り、seedを固定することで、実行のたびに同じ結果を再現します。

import platform

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

SEED = 202606
rng = np.random.default_rng(SEED)
ALPHA = 0.05

pd.set_option("display.precision", 4)
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.unicode_minus"] = False

print(f"Python     : {platform.python_version()}")
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
NumPy      : 2.5.1
pandas     : 3.0.3
SciPy      : 1.18.0
matplotlib : 3.11.0
乱数seed   : 202606

架空データの作成

一つの工場内で行われた工程能力確認、設備改善、品質分類、仕入先評価、保全評価を想定します。外部ファイルは使わず、連続量は正規分布または右裾の長い分布、カテゴリデータは集計表、時系列残差は自己回帰過程から作成します。

データ生成時には差を埋め込んでいますが、実務では「期待した結果が出るまでデータを取り直す」ことはできません。分析前に標本数、主要指標、除外条件を決めることが重要です。

# No.051: 充填量(g)、基準500 g
fill_weight = rng.normal(loc=500.48, scale=1.10, size=40)

# No.052: 同一設備12台の1バッチ当たり消費電力(kWh)
energy_before = rng.normal(loc=126.0, scale=5.0, size=12)
energy_after = energy_before - rng.normal(loc=4.2, scale=2.0, size=12)

# No.053: ラインA/Bの軸径偏差(μm)。Bは標本数が少なく分散が大きい
diameter_a = rng.normal(loc=0.00, scale=0.16, size=36)
diameter_b = rng.normal(loc=0.16, scale=0.38, size=22)

# No.054: 行=シフト、列=不適合分類
defect_table = pd.DataFrame(
    [[18, 7, 5], [8, 17, 15]],
    index=["昼勤", "夜勤"], columns=["外観", "寸法", "異物"]
)

# No.055: 行=旧/新ファームウェア、列=停止警報あり/なし
alarm_table = pd.DataFrame(
    [[7, 23], [1, 29]],
    index=["旧版", "新版"], columns=["停止警報あり", "停止警報なし"]
)

# No.056: 3仕入先の破断荷重(kN)
strength = {
    "仕入先A": rng.normal(50.00, 0.18, 20),
    "仕入先B": rng.normal(50.10, 0.18, 20),
    "仕入先C": rng.normal(50.32, 0.18, 20),
}

# No.057: 独立した2治具の段取時間(分、右裾が長い)
setup_old = rng.lognormal(mean=np.log(18.0), sigma=0.28, size=28)
setup_new = rng.lognormal(mean=np.log(15.5), sigma=0.25, size=26)

# No.058: 3保全班の復旧時間(分、外れ値を含み得る)
recovery = {
    "班A": rng.lognormal(np.log(38), 0.30, 18),
    "班B": rng.lognormal(np.log(31), 0.28, 18),
    "班C": rng.lognormal(np.log(47), 0.35, 18),
}

# No.059: サイクル時間(秒)。事前に固定した基準分布 N(45, 2.5^2) と比較
cycle_time = 40.5 + rng.gamma(shape=2.0, scale=2.25, size=80)

# No.060: 炉温予測モデルの残差(℃)。AR(1)として連続性を持たせる
innovations = rng.normal(0, 0.55, 120)
temp_residual = np.zeros(120)
for t in range(1, len(temp_residual)):
    temp_residual[t] = 0.62 * temp_residual[t - 1] + innovations[t]

data_overview = pd.DataFrame({
    "分析テーマ": ["充填量", "設備電力", "軸径偏差", "不適合分類", "停止警報",
                 "破断荷重", "段取時間", "復旧時間", "サイクル時間", "炉温モデル残差"],
    "観測単位/群": [len(fill_weight), len(energy_before), f"{len(diameter_a)} / {len(diameter_b)}",
                int(defect_table.to_numpy().sum()), int(alarm_table.to_numpy().sum()),
                " / ".join(str(len(v)) for v in strength.values()),
                f"{len(setup_old)} / {len(setup_new)}",
                " / ".join(str(len(v)) for v in recovery.values()), len(cycle_time), len(temp_residual)],
    "使用する検定": ["1標本t", "対応t", "Welch", "カイ二乗", "Fisher",
                "ANOVA", "Mann–Whitney", "Kruskal–Wallis", "KS", "Ljung–Box"],
})
display(data_overview)
分析テーマ 観測単位/群 使用する検定
0 充填量 40 1標本t
1 設備電力 12 対応t
2 軸径偏差 36 / 22 Welch
3 不適合分類 70 カイ二乗
4 停止警報 60 Fisher
5 破断荷重 20 / 20 / 20 ANOVA
6 段取時間 28 / 26 Mann–Whitney
7 復旧時間 18 / 18 / 18 Kruskal–Wallis
8 サイクル時間 80 KS
9 炉温モデル残差 120 Ljung–Box

No.051:t検定 — 充填量の平均が基準値と異なるか

実務での意味

充填工程の平均が表示量500 gからずれているかを確認します。過充填は材料損失、過少充填は顧客要求違反につながるため、工程中心の調整判断に関係します。ただし、規格内かどうかの保証には工程能力分析も別途必要です。

分析・モデル化の考え方

帰無仮説を H0:μ=500H_0:\mu=500、対立仮説を H1:μ500H_1:\mu\neq500 とします。母分散が未知の1標本t統計量は

t=xˉμ0s/nt=\frac{\bar{x}-\mu_0}{s/\sqrt{n}}

です。観測の独立性と、平均の標本分布がt近似できることを仮定します。p値に加え、平均差、95%信頼区間、標準化効果量Cohen’s ddを確認します。

Pythonで確認する

target_weight = 500.0
t51 = stats.ttest_1samp(fill_weight, popmean=target_weight)
mean51 = fill_weight.mean()
sd51 = fill_weight.std(ddof=1)
ci51 = stats.t.interval(0.95, len(fill_weight) - 1, loc=mean51, scale=stats.sem(fill_weight))
d51 = (mean51 - target_weight) / sd51

display(pd.DataFrame({
    "平均[g]": [mean51], "基準との差[g]": [mean51 - target_weight],
    "95%CI下限": [ci51[0]], "95%CI上限": [ci51[1]],
    "t値": [t51.statistic], "p値": [t51.pvalue], "Cohen's d": [d51],
    "判定(alpha=0.05)": ["基準と差あり" if t51.pvalue < ALPHA else "差を確認できず"]
}))

fig, ax = plt.subplots()
ax.hist(fill_weight, bins=10, color="#4C78A8", edgecolor="white", alpha=0.85)
ax.axvline(target_weight, color="#E45756", linestyle="--", label="基準 500 g")
ax.axvline(mean51, color="#2A9D8F", linewidth=2, label=f"標本平均 {mean51:.2f} g")
ax.set_title("充填量の分布と基準値")
ax.set_xlabel("充填量 [g]")
ax.set_ylabel("個数")
ax.grid(axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
平均[g] 基準との差[g] 95%CI下限 95%CI上限 t値 p値 Cohen's d 判定(alpha=0.05)
0 500.3297 0.3297 499.9755 500.6839 1.8827 0.0672 0.2977 差を確認できず

png

結果の読み取り

出力では平均が500 gを上回り、95%信頼区間とp値から工程中心のずれを評価できます。ただし、統計的に差があっても、調整コストや計量器の不確かさに比べて差が小さければ、即時調整が最善とは限りません。最小重要差を事前に決め、管理図で時間方向の安定性も確認したうえで調整します。

No.052:対応のあるt検定 — 同一設備の省エネ設定前後を比べる

実務での意味

同じ12台の設備について、設定変更前後の1バッチ当たり消費電力を比較します。設備ごとの基礎的な消費量の違いを差し引くことで、設定変更そのものの効果を評価しやすくなります。

分析・モデル化の考え方

設備 ii の差を di=xbefore,ixafter,id_i=x_{\mathrm{before},i}-x_{\mathrm{after},i} とすると、対応のあるt検定は H0:μd=0H_0:\mu_d=0 を検定する1標本t検定です。重要な仮定は、前後それぞれの正規性ではなく、did_i の独立性と近似的な正規性です。同じ設備の対応を崩してはいけません。

Pythonで確認する

diff52 = energy_before - energy_after
t52 = stats.ttest_rel(energy_before, energy_after)
ci52 = stats.t.interval(0.95, len(diff52) - 1, loc=diff52.mean(), scale=stats.sem(diff52))

display(pd.DataFrame({
    "変更前平均[kWh]": [energy_before.mean()], "変更後平均[kWh]": [energy_after.mean()],
    "平均削減量[kWh]": [diff52.mean()], "削減量95%CI下限": [ci52[0]],
    "削減量95%CI上限": [ci52[1]], "t値": [t52.statistic], "p値": [t52.pvalue],
    "判定": ["前後差あり" if t52.pvalue < ALPHA else "前後差を確認できず"]
}))

fig, ax = plt.subplots()
for i, (before, after) in enumerate(zip(energy_before, energy_after), start=1):
    ax.plot(["変更前", "変更後"], [before, after], marker="o", alpha=0.55)
ax.set_title("同一設備の設定変更前後の消費電力")
ax.set_xlabel("設定")
ax.set_ylabel("1バッチ当たり消費電力 [kWh]")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
変更前平均[kWh] 変更後平均[kWh] 平均削減量[kWh] 削減量95%CI下限 削減量95%CI上限 t値 p値 判定
0 124.6055 120.4945 4.1109 2.7853 5.4366 6.8255 2.8555e-05 前後差あり

png

結果の読み取り

多くの線が右下がりで、平均削減量の信頼区間が0をまたがなければ、省エネ効果を支持する材料になります。横展開前には、生産量、品種、外気温、立上げ時間など前後で変わった条件を確認します。削減電力量を年間バッチ数と電力単価へ換算すると、実務的な投資効果として説明できます。

No.053:Welch検定 — ばらつきの異なる2ラインの寸法平均を比べる

実務での意味

ラインAとBの軸径偏差を比較します。Bは標本数が少なく、設備状態や材料ロットの影響でばらつきも大きい想定です。ライン間補正を行う前に、平均差が偶然の範囲かを評価します。

分析・モデル化の考え方

Welchのt検定は、独立2群の平均差を、等分散を仮定せずに検定します。統計量は

t=xˉAxˉBsA2/nA+sB2/nBt=\frac{\bar{x}_A-\bar{x}_B}{\sqrt{s_A^2/n_A+s_B^2/n_B}}

で、自由度はWelch–Satterthwaite近似で求めます。設備から連続して採った製品に強い系列相関がある場合は、ロット単位へ集約するなど、独立な分析単位を設計します。

Pythonで確認する

t53 = stats.ttest_ind(diameter_a, diameter_b, equal_var=False)
va, vb = diameter_a.var(ddof=1), diameter_b.var(ddof=1)
se53 = np.sqrt(va / len(diameter_a) + vb / len(diameter_b))
df53 = (va / len(diameter_a) + vb / len(diameter_b)) ** 2 / (
    (va / len(diameter_a)) ** 2 / (len(diameter_a) - 1)
    + (vb / len(diameter_b)) ** 2 / (len(diameter_b) - 1)
)
diff53 = diameter_a.mean() - diameter_b.mean()
ci53 = diff53 + np.array([-1, 1]) * stats.t.ppf(0.975, df53) * se53

display(pd.DataFrame({
    "A平均[μm]": [diameter_a.mean()], "B平均[μm]": [diameter_b.mean()],
    "A標準偏差": [np.sqrt(va)], "B標準偏差": [np.sqrt(vb)],
    "平均差A-B": [diff53], "差の95%CI下限": [ci53[0]], "差の95%CI上限": [ci53[1]],
    "Welch自由度": [df53], "t値": [t53.statistic], "p値": [t53.pvalue],
}))

fig, ax = plt.subplots()
ax.boxplot([diameter_a, diameter_b], tick_labels=["ラインA", "ラインB"], showmeans=True)
ax.axhline(0, color="#E45756", linestyle="--", label="設計中心")
ax.set_title("ライン別の軸径偏差")
ax.set_xlabel("製造ライン")
ax.set_ylabel("軸径偏差 [μm]")
ax.grid(axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
A平均[μm] B平均[μm] A標準偏差 B標準偏差 平均差A-B 差の95%CI下限 差の95%CI上限 Welch自由度 t値 p値
0 -0.0185 0.1991 0.1734 0.3552 -0.2176 -0.3839 -0.0513 27.2124 -2.6841 0.0122

png

結果の読み取り

Welch検定は標本数と分散の違いを織り込んだ平均差の評価です。箱ひげ図からは、平均差だけでなくBのばらつきが大きいことも読み取れます。平均補正だけでは品質リスクを解消できない可能性があるため、Bではばらつき要因の層別と工程能力の確認を優先します。

No.054:カイ二乗検定 — シフトと不適合分類の関連を調べる

実務での意味

昼勤と夜勤で、外観・寸法・異物という不適合の構成が同じかを調べます。総不適合件数だけでは見えない構成比の違いを把握し、照明、段取り、清掃、検査教育などの重点を決めます。

分析・モデル化の考え方

独立性のカイ二乗検定では、帰無仮説を「シフトと不適合分類は独立」とします。期待度数は

Eij=(i行の合計)(j列の合計)総数,χ2=i,j(OijEij)2EijE_{ij}=\frac{(i\text{行の合計})(j\text{列の合計})}{\text{総数}},\qquad \chi^2=\sum_{i,j}\frac{(O_{ij}-E_{ij})^2}{E_{ij}}

です。期待度数が小さいセルが多い場合は近似が不安定になります。標準化残差 (OE)/E(O-E)/\sqrt{E} を見て、どの組合せが差へ寄与したかを確認します。

Pythonで確認する

chi54, p54, dof54, expected54 = stats.chi2_contingency(defect_table)
expected54 = pd.DataFrame(expected54, index=defect_table.index, columns=defect_table.columns)
residual54 = (defect_table - expected54) / np.sqrt(expected54)

display(pd.concat({"観測度数": defect_table, "期待度数": expected54.round(2),
                   "Pearson残差": residual54.round(2)}, axis=1))
display(pd.DataFrame({"カイ二乗値": [chi54], "自由度": [dof54], "p値": [p54],
                      "最小期待度数": [expected54.min().min()],
                      "判定": ["関連あり" if p54 < ALPHA else "関連を確認できず"]}))

rates54 = defect_table.div(defect_table.sum(axis=1), axis=0)
ax = rates54.plot(kind="bar", color=["#4C78A8", "#F2CF5B", "#E45756"])
ax.set_title("シフト別の不適合分類構成比")
ax.set_xlabel("勤務シフト")
ax.set_ylabel("シフト内構成比")
ax.grid(axis="y", alpha=0.3)
ax.legend(title="不適合分類")
plt.xticks(rotation=0)
plt.tight_layout()
plt.show()
観測度数 期待度数 Pearson残差
外観 寸法 異物 外観 寸法 異物 外観 寸法 異物
昼勤 18 7 5 11.14 10.29 8.57 2.05 -1.02 -1.22
夜勤 8 17 15 14.86 13.71 11.43 -1.78 0.89 1.06
カイ二乗値 自由度 p値 最小期待度数 判定
0 11.8256 2 0.0027 8.5714 関連あり

png

結果の読み取り

p値が小さければ、シフトと不適合分類を独立とはみなしにくいと判断します。構成比と残差から、昼勤では外観、夜勤では寸法・異物が相対的に多い構造を確認できます。ただしシフトそのものが原因とは限らず、品種構成、材料ロット、設備、検査者が交絡していないかを追跡します。

No.055:Fisher正確検定 — 少数の停止警報を2×2表で評価する

実務での意味

センサー制御の旧版・新版ファームウェアを各30設備で試し、停止警報の有無を集計します。発生件数が少ない段階でも、新版が警報率に関係するかを評価したい場面です。

分析・モデル化の考え方

Fisher正確検定は、2×2分割表の周辺度数を固定した条件付き分布から、観測表以上に極端な表の確率を正確に計算します。小標本や期待度数が小さい場合に有用です。オッズ比

OR=adbc\mathrm{OR}=\frac{ad}{bc}

も併記します。ここでは「旧版の警報オッズ / 新版の警報オッズ」なので、1より大きいほど旧版の警報が多い方向です。

Pythonで確認する

odds55, p55 = stats.fisher_exact(alarm_table.to_numpy(), alternative="two-sided")
rates55 = alarm_table["停止警報あり"] / alarm_table.sum(axis=1)

display(alarm_table)
display(pd.DataFrame({
    "旧版警報率": [rates55["旧版"]], "新版警報率": [rates55["新版"]],
    "警報率差(旧-新)": [rates55["旧版"] - rates55["新版"]],
    "オッズ比(旧/新)": [odds55], "両側p値": [p55],
    "判定": ["関連あり" if p55 < ALPHA else "関連を確認できず"]
}))

fig, ax = plt.subplots()
ax.bar(rates55.index, rates55.values, color=["#E45756", "#2A9D8F"])
ax.set_title("ファームウェア別の停止警報率")
ax.set_xlabel("ファームウェア")
ax.set_ylabel("停止警報率")
ax.set_ylim(0, max(rates55.max() * 1.3, 0.1))
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()
停止警報あり 停止警報なし
旧版 7 23
新版 1 29
旧版警報率 新版警報率 警報率差(旧-新) オッズ比(旧/新) 両側p値 判定
0 0.2333 0.0333 0.2 8.8261 0.0523 関連を確認できず

png

結果の読み取り

新版の警報率は低く、オッズ比も効果の方向を示します。一方、件数が少ないためp値は境界的になり得ます。「有意でない=同等」とは結論できません。安全に関わる指標では、効果方向、許容リスク、追加試験の必要台数、重大事象の内容を合わせて段階導入を判断します。

No.056:ANOVA — 3仕入先の破断荷重平均を比較する

実務での意味

同一仕様の部材を供給する3社について、破断荷重の平均に違いがあるかを比較します。仕入先認定や受入検査水準の見直しに関係します。

分析・モデル化の考え方

一元配置分散分析は、H0:μA=μB=μCH_0:\mu_A=\mu_B=\mu_C を、群間変動と群内変動の比

F=MSbetweenMSwithinF=\frac{MS_{\mathrm{between}}}{MS_{\mathrm{within}}}

で検定します。独立性、各群の近似的な正規性、等分散性が前提です。ANOVAが有意でも「全組合せが異なる」とは限らないため、事前対比または多重比較が必要です。ここでは全体効果量 η2=SSbetween/SStotal\eta^2=SS_{\mathrm{between}}/SS_{\mathrm{total}} とTukey HSDも確認します。

Pythonで確認する

f56 = stats.f_oneway(*strength.values())
all56 = np.concatenate(list(strength.values()))
grand56 = all56.mean()
ss_between56 = sum(len(x) * (x.mean() - grand56) ** 2 for x in strength.values())
ss_total56 = ((all56 - grand56) ** 2).sum()
eta56 = ss_between56 / ss_total56
tukey56 = stats.tukey_hsd(*strength.values())

display(pd.DataFrame({
    "n": [len(x) for x in strength.values()],
    "平均[kN]": [x.mean() for x in strength.values()],
    "標準偏差[kN]": [x.std(ddof=1) for x in strength.values()],
}, index=strength.keys()))
display(pd.DataFrame({"F値": [f56.statistic], "p値": [f56.pvalue], "eta二乗": [eta56],
                      "判定": ["少なくとも1群に差あり" if f56.pvalue < ALPHA else "差を確認できず"]}))
display(pd.DataFrame(tukey56.pvalue, index=strength.keys(), columns=strength.keys()).round(4))

fig, ax = plt.subplots()
ax.boxplot(list(strength.values()), tick_labels=list(strength.keys()), showmeans=True)
ax.set_title("仕入先別の破断荷重")
ax.set_xlabel("仕入先")
ax.set_ylabel("破断荷重 [kN]")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()
n 平均[kN] 標準偏差[kN]
仕入先A 20 49.9516 0.1885
仕入先B 20 50.1246 0.2153
仕入先C 20 50.3190 0.1516
F値 p値 eta二乗 判定
0 19.3208 3.9264e-07 0.404 少なくとも1群に差あり
仕入先A 仕入先B 仕入先C
仕入先A 1.0000 0.0135 0.0000
仕入先B 0.0135 1.0000 0.0049
仕入先C 0.0000 0.0049 1.0000

png

結果の読み取り

ANOVAのp値は全体として平均差があるか、η2\eta^2 は仕入先が総変動のどの程度を説明するかを示します。Tukey HSDのp値行列で差のある組合せを絞れます。ただし、平均が高い仕入先が常に最適とは限りません。規格下限に対する余裕、ばらつき、価格、納期、ロット間変動も含めて認定基準を設計します。

No.057:Mann–Whitney検定 — 2治具の段取時間を順位で比べる

実務での意味

旧治具と新治具を別の作業で使用し、段取時間の分布に差があるかを評価します。段取時間は0未満にならず、トラブル時に右裾が長くなりやすいため、平均値と正規分布の仮定が扱いにくい指標です。

分析・モデル化の考え方

Mann–WhitneyのU検定は、独立2群をまとめて順位化し、順位の偏りを評価するノンパラメトリック検定です。2群の分布形状が同じとみなせる場合は位置の差として解釈できますが、形状も違う場合は「中央値だけの検定」とは言えません。ここでは中央値、四分位範囲、順位二分相関も併記します。

Pythonで確認する

u57 = stats.mannwhitneyu(setup_old, setup_new, alternative="two-sided")
rbc57 = 2 * u57.statistic / (len(setup_old) * len(setup_new)) - 1

summary57 = pd.DataFrame({
    "n": [len(setup_old), len(setup_new)],
    "中央値[分]": [np.median(setup_old), np.median(setup_new)],
    "第1四分位": [np.quantile(setup_old, .25), np.quantile(setup_new, .25)],
    "第3四分位": [np.quantile(setup_old, .75), np.quantile(setup_new, .75)],
}, index=["旧治具", "新治具"])
display(summary57)
display(pd.DataFrame({"U値": [u57.statistic], "p値": [u57.pvalue],
                      "順位二分相関": [rbc57],
                      "判定": ["分布位置に差あり" if u57.pvalue < ALPHA else "差を確認できず"]}))

fig, ax = plt.subplots()
for values, label, color in [(setup_old, "旧治具", "#E45756"), (setup_new, "新治具", "#2A9D8F")]:
    x = np.sort(values)
    y = np.arange(1, len(x) + 1) / len(x)
    ax.step(x, y, where="post", label=label, color=color)
ax.set_title("治具別の段取時間の経験累積分布")
ax.set_xlabel("段取時間 [分]")
ax.set_ylabel("累積割合")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
n 中央値[分] 第1四分位 第3四分位
旧治具 28 18.6360 15.0231 22.3628
新治具 26 15.6962 14.1175 20.4698
U値 p値 順位二分相関 判定
0 437.0 0.2094 0.2005 差を確認できず

png

結果の読み取り

経験累積分布が左にあるほど短時間で段取が完了する傾向です。p値に加え、中央値差と順位二分相関から差の大きさを確認します。独立群の比較なので、作業者熟練度や品種難易度が偏ると治具効果と混ざります。可能なら同じ作業者・同程度の品種で無作為な使用順序を設計します。

No.058:Kruskal–Wallis検定 — 3保全班の復旧時間を順位で比べる

実務での意味

3つの保全班について故障復旧時間を比較し、手順標準化や教育の重点を検討します。復旧時間は故障の難易度によって外れ値が生じやすく、正規性・等分散性が疑わしい想定です。

分析・モデル化の考え方

Kruskal–Wallis検定は、独立3群以上をまとめて順位化し、群ごとの平均順位の違いを評価します。帰無仮説は各群の分布が同じことです。検定統計量 HH は大標本で自由度 k1k-1 のカイ二乗分布に近似します。全体検定が有意でも、どの班間が違うかには多重比較が別途必要です。

Pythonで確認する

h58 = stats.kruskal(*recovery.values())
n58 = sum(len(x) for x in recovery.values())
k58 = len(recovery)
epsilon58 = max(0, (h58.statistic - k58 + 1) / (n58 - k58))

display(pd.DataFrame({
    "n": [len(x) for x in recovery.values()],
    "中央値[分]": [np.median(x) for x in recovery.values()],
    "平均順位": [stats.rankdata(np.concatenate(list(recovery.values())))[
        sum(len(v) for v in list(recovery.values())[:i]):sum(len(v) for v in list(recovery.values())[:i+1])
    ].mean() for i in range(k58)],
}, index=recovery.keys()))
display(pd.DataFrame({"H値": [h58.statistic], "自由度": [k58 - 1], "p値": [h58.pvalue],
                      "epsilon二乗": [epsilon58],
                      "判定": ["少なくとも1群に差あり" if h58.pvalue < ALPHA else "差を確認できず"]}))

fig, ax = plt.subplots()
ax.boxplot(list(recovery.values()), tick_labels=list(recovery.keys()), showmeans=True)
ax.set_title("保全班別の故障復旧時間")
ax.set_xlabel("保全班")
ax.set_ylabel("復旧時間 [分]")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()
n 中央値[分] 平均順位
班A 18 39.0264 31.2222
班B 18 29.2456 17.6111
班C 18 38.9934 33.6667
H値 自由度 p値 epsilon二乗 判定
0 10.8853 2 0.0043 0.1742 少なくとも1群に差あり

png

結果の読み取り

全体p値と効果量から、班による分布差が無視できるかを判断します。中央値が長い班を直ちに「能力不足」とみなしてはいけません。担当設備、故障モード、夜間呼出し、部品待ち時間の構成が違う可能性があります。故障難易度を層別し、全体検定後の事後比較ではp値補正を行います。

No.059:Kolmogorov–Smirnov検定 — 基準分布とのずれを検出する

実務での意味

サイクル時間が、設備導入時に固定した基準分布 N(45,2.52)N(45, 2.5^2) と一致するかを確認します。平均と標準偏差だけでは見落とし得る、歪み、裾、局所的な滞留を分布全体で捉えます。

分析・モデル化の考え方

1標本Kolmogorov–Smirnov検定の統計量は、経験累積分布 Fn(x)F_n(x) と基準累積分布 F0(x)F_0(x) の最大差

D=supxFn(x)F0(x)D=\sup_x\left|F_n(x)-F_0(x)\right|

です。連続分布を仮定します。重要なのは、基準分布の平均・標準偏差を今回の標本から推定していないことです。同じデータで推定すると通常のKS検定のp値は適切でなく、Lilliefors型の補正やブートストラップが必要です。

Pythonで確認する

target_mu59, target_sd59 = 45.0, 2.5
target_cdf59 = lambda x: stats.norm.cdf(x, loc=target_mu59, scale=target_sd59)
ks59 = stats.kstest(cycle_time, target_cdf59)
x59 = np.sort(cycle_time)
ecdf59 = np.arange(1, len(x59) + 1) / len(x59)
grid59 = np.linspace(min(x59.min(), 37.5), max(x59.max(), 52.5), 300)

display(pd.DataFrame({
    "標本数": [len(cycle_time)], "標本平均[秒]": [cycle_time.mean()],
    "標本標準偏差": [cycle_time.std(ddof=1)], "基準平均": [target_mu59],
    "基準標準偏差": [target_sd59], "KS統計量D": [ks59.statistic], "p値": [ks59.pvalue],
    "判定": ["基準分布と差あり" if ks59.pvalue < ALPHA else "差を確認できず"]
}))

fig, ax = plt.subplots()
ax.step(x59, ecdf59, where="post", label="観測ECDF", color="#4C78A8")
ax.plot(grid59, stats.norm.cdf(grid59, target_mu59, target_sd59),
        label="基準 N(45, 2.5²)", color="#E45756", linestyle="--")
ax.set_title("サイクル時間の経験分布と基準分布")
ax.set_xlabel("サイクル時間 [秒]")
ax.set_ylabel("累積確率")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
標本数 標本平均[秒] 標本標準偏差 基準平均 基準標準偏差 KS統計量D p値 判定
0 80 44.8608 2.5536 45.0 2.5 0.0838 0.5986 差を確認できず

png

結果の読み取り

ECDFと基準CDFの縦方向の最大差がKS統計量です。今回のような右裾の長い分布では、平均だけでなく遅いサイクルの蓄積を検出できます。差があれば、短時間停止、材料供給待ち、作業介入などの混合を疑います。KS検定は差の場所や原因までは示さないため、分位点差やログと突き合わせます。

No.060:Ljung–Box検定 — 炉温モデル残差の自己相関を調べる

実務での意味

炉温予測モデルの残差が時間方向に連続していないかを調べます。残差に自己相関が残ると、モデルが熱慣性や制御周期を捉えておらず、異常判定の連続誤報や予測区間の過小評価につながります。

分析・モデル化の考え方

Ljung–Box検定は、ラグ1から hh までの自己相関をまとめて検定します。

Q=n(n+2)k=1hρ^k2nkQ=n(n+2)\sum_{k=1}^{h}\frac{\hat{\rho}_k^2}{n-k}

帰無仮説は「指定ラグまでの自己相関がすべて0」です。QQ は近似的にカイ二乗分布に従います。ARIMA等の推定済みモデル残差へ適用する場合、自由度は推定したAR・MAパラメータ数を差し引く必要があります。ここでは説明用の固定予測モデル残差として自由度を hh とします。

Pythonで確認する

def ljung_box(x, max_lag):
    x = np.asarray(x) - np.mean(x)
    n = len(x)
    denominator = np.dot(x, x)
    acf = np.array([np.dot(x[k:], x[:-k]) / denominator for k in range(1, max_lag + 1)])
    lags = np.arange(1, max_lag + 1)
    q_values = n * (n + 2) * np.cumsum(acf ** 2 / (n - lags))
    p_values = stats.chi2.sf(q_values, df=lags)
    return lags, acf, q_values, p_values

lags60, acf60, q60, p60 = ljung_box(temp_residual, max_lag=12)
result60 = pd.DataFrame({"ラグ": lags60, "自己相関": acf60,
                         "累積Q値": q60, "p値": p60})
display(result60.round(4))
display(pd.DataFrame({"評価ラグ": [12], "Ljung-Box Q": [q60[-1]], "p値": [p60[-1]],
                      "判定": ["自己相関あり" if p60[-1] < ALPHA else "自己相関を確認できず"]}))

fig, ax = plt.subplots()
ax.stem(lags60, acf60, basefmt=" ")
bound60 = 1.96 / np.sqrt(len(temp_residual))
ax.axhline(bound60, color="#E45756", linestyle="--", label="参考95%限界")
ax.axhline(-bound60, color="#E45756", linestyle="--")
ax.axhline(0, color="black", linewidth=0.8)
ax.set_title("炉温モデル残差の自己相関")
ax.set_xlabel("ラグ")
ax.set_ylabel("自己相関係数")
ax.set_xticks(lags60)
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
ラグ 自己相関 累積Q値 p値
0 1 0.5605 38.6467 0.0
1 2 0.2884 48.9647 0.0
2 3 0.1936 53.6523 0.0
3 4 0.1218 55.5232 0.0
4 5 0.1156 57.2235 0.0
5 6 0.0181 57.2658 0.0
6 7 -0.0316 57.3950 0.0
7 8 -0.0229 57.4635 0.0
8 9 -0.0249 57.5455 0.0
9 10 0.1029 58.9560 0.0
10 11 0.0745 59.7013 0.0
11 12 0.0141 59.7282 0.0
評価ラグ Ljung-Box Q p値 判定
0 12 59.7282 2.5300e-08 自己相関あり

png

結果の読み取り

小さいラグで正の自己相関が見られ、ラグ12までのp値が小さければ、残差を独立なノイズとはみなしにくいと判断します。炉の熱慣性、設定変更履歴、搬送周期、センサー平滑化などを特徴量または時系列モデルへ組み込みます。Ljung–Box検定は「どのモデルが正しいか」ではなく、残差に未説明の時間構造が残ることを知らせます。

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

10種類の検定は、対象とする問いによって役割が分かれます。

  • 平均を問う場合も、1群、対応あり、独立2群、3群以上で手法が変わる
  • カテゴリ比率では、期待度数が十分ならカイ二乗、小標本2×2表ならFisherが候補になる
  • 外れ値や歪みが強い連続量では、順位検定が有力だが、検定対象の解釈に注意する
  • 平均・中央値だけでなく分布全体を問うならKS検定、時系列の独立性を問うならLjung–Box検定を使う
  • p値は効果の大きさでも、帰無仮説が正しい確率でもない。差、信頼区間、効果量を併記する
  • 「有意差なし」は同等性の証明ではない。同等性を示したい場合は許容差を定めた同等性検定を設計する

最も重要なのは、データを見てから都合のよい検定を選ぶのではなく、現場の意思決定、比較単位、最小重要差を先に定義することです。

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

  1. 測定系を確認する:校正、分解能、測定者差、欠測・丸めを評価する
  2. 分析単位を揃える:製品を独立標本とみなせるか、ロット・設備単位にすべきか決める
  3. 事前計画を残す:仮説、片側・両側、有意水準、必要標本数、除外規則を記録する
  4. 実務閾値を決める:統計的有意差とは別に、採用に必要な最小効果を設定する
  5. 交絡を管理する:品種、材料、作業者、時間帯、設備状態を無作為化・層別化する
  6. 多重性へ対応する:指標・群・期間を多数比較する場合は主要仮説を絞り、補正する
  7. 継続監視へ接続する:一度の検定で終わらせず、管理図やKPIで効果の持続を追う
  8. 再現可能にする:元データ、コード、環境、実行日時、分析者、判断を版管理する

安全・法規・顧客保証に関わる変更は、p値だけで自動承認せず、FMEA、工程能力、変更管理、専門家レビューと組み合わせます。

まとめ

No.051〜No.060では、製造現場の代表的な比較課題を、t検定、対応のあるt検定、Welch検定、カイ二乗検定、Fisher正確検定、ANOVA、Mann–Whitney検定、Kruskal–Wallis検定、Kolmogorov–Smirnov検定、Ljung–Box検定で確認しました。

検定結果を意思決定へ変えるには、手法のAPIよりも、問いとデータ構造を合わせることが重要です。p値、効果量、信頼区間、グラフ、現場の許容差を一つの判断資料にまとめることで、改善効果を過大評価せず、見逃しも抑えられます。

法人向けのご相談

数理工房では、製造業における品質データ分析、実験計画、工程改善効果の検証、異常検知、需要・在庫・生産計画、社内データ人材育成をご支援しています。

「現場データに合う検定を選べない」「有意差は出たが経営判断へつながらない」「Excel分析を再現可能なPython運用へ移したい」といった段階から、課題整理、分析設計、実装、運用定着までご相談いただけます。

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