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

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

2020-04-14

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

概要

製造現場では、温度、圧力、寸法、粗さ、電力などが同時に変動します。個別KPIの平均だけを追うと、組合せで現れる品質リスクや、条件を追加した後の予測変化を見落とします。本記事では、架空の精密成形工場における900ロットのデータを使い、同時分布・周辺分布・条件付き分布・独立性・共分散行列・相関行列・多変量正規分布・条件付き正規分布・ガウス過程・コピュラを、工程監視と意思決定へ結び付けます。

対象は「確率・統計理論100本ノック」の No.061〜No.070 です。各手法を単なる計算で終わらせず、「どの条件の組合せを監視するか」「追加情報で予測がどう変わるか」「同時リスクにどれだけ余力を持つか」という実務上の問いへつなげます。

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

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

題材は、製品A・B・Cを共通設備で生産する精密成形工場です。管理者は次の判断を求められています。

寸法と表面粗さが同時に悪化するロットを、単独KPIより早く捉えられるか
温度帯が分かったとき、品質状態の確率をどう更新するか
2つの異常が偶然重なったのか、依存関係があるのか
複数センサーを行列で管理し、冗長な指標や連動方向を発見できるか
温度を観測した後の電力原単位を、点ではなく分布で予測できるか
時間とともに滑らかに変化する設備ドリフトを、不確実性付きで補間できるか
停止時間と納期損失の周辺分布を変えず、同時発生リスクだけを表現できるか

多変量確率の役割は「列を増やすこと」ではありません。複数KPIが作る確率構造を、監視ルール、保全計画、品質保証、BCPの判断へ翻訳することです。

現場でよくある状況

寸法と粗さを別々の管理図で監視し、両方がやや悪いロットを見逃す
全体の不良率を報告するが、高温時・夜勤時など条件別の発生確率を示していない
相関が0に近いことを「独立」と読み替える
単位の大きい変数が支配する共分散行列を、そのまま関係の強さの比較に使う
多変量正規分布を、分布形状を確認せず全データへ機械的に当てはめる
補間曲線だけを見て、未観測区間の不確実性を無視する
停止時間と損失額を独立に乱数生成し、同時テールリスクを過小評価する

どれも計算そのものより、前提と用途の取り違えが問題です。本記事では、表・式・図を往復しながら、モデルが答えられる問いと答えられない問いを明確にします。

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

変数が2つになるだけで、確率は「それぞれの分布」だけでは決まりません。 $X$ と $Y$ の周辺分布が同じでも、同時分布 $p(x,y)$ は複数あり得ます。依存構造が違えば、同時異常や合計損失の確率も変わります。

また、条件が追加されると評価対象が変わります。全体確率 $P(Y)$ と、温度帯が分かった後の $P(Y\mid X)$ は別物です。後者は現場で観測できる情報を使った予測ですが、条件を細かくしすぎると該当ロットが減り、推定が不安定になります。

さらに、正規分布・ガウス過程・コピュラはいずれも「ガウス」という語に関係しますが、役割が異なります。多変量正規分布は同時ベクトル、条件付き正規分布は観測後の更新、ガウス過程は関数全体、ガウスコピュラは周辺分布と依存構造の分離を扱います。

今回扱うノックの全体像

No.	テーマ	製造業での問い
061	同時分布	寸法状態と粗さ状態の組合せは何％か
062	周辺分布	同時表から各KPI単独の分布を取り出せるか
063	条件付き分布	温度帯を知ると粗さ状態の確率はどう変わるか
064	独立性	高温と粗さ異常は独立とみなせるか
065	共分散行列	複数KPIのばらつきと共変動を一括表現できるか
066	相関行列	単位を除いて関係の強さを比較できるか
067	多変量正規分布	複数KPIの通常領域を楕円で捉えられるか
068	条件付き正規分布	温度観測後の電力分布を更新できるか
069	ガウス過程	設備ドリフトを不確実性付きで補間できるか
070	コピュラ	非正規な停止時間と損失の同時リスクを表せるか

No.061〜No.064で確率表の基本、No.065〜No.068でベクトルの分布、No.069〜No.070で関数と依存構造のモデルへ進みます。

Python 環境の準備

NumPyで乱数・行列計算、pandasで集計、matplotlibで可視化、SciPyで正規分布の分位点を扱います。外部データやseabornには依存しません。再実行して同じ結果になるよう、乱数seedを固定します。

import platform
import sys

import japanize_matplotlib  # noqa: F401  日本語フォント設定
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display
from scipy.stats import chi2, norm

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

pd.set_option("display.max_columns", 30)
pd.set_option("display.float_format", lambda x: f"{x:,.4f}")
plt.rcParams["figure.figsize"] = (8, 4.8)
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"matplotlib: {matplotlib.__version__}")
print(f"乱数 seed: {SEED}")

Python: 3.11.9
OS: Darwin 25.3.0
NumPy: 1.26.4
pandas: 2.2.2
matplotlib: 3.9.2
乱数 seed: 20260711

架空データの作成

900ロットについて、製品、勤務帯、金型温度、射出圧力、寸法偏差、表面粗さ、電力原単位を生成します。共通する潜在変動を複数KPIへ反映し、現場で起こる連動を模擬します。

高温側では寸法偏差・粗さ・電力が上がりやすい
圧力上昇は寸法偏差と粗さを抑える方向に働く
製品と勤務帯で平均水準が少し異なる
規格判定は、寸法偏差 $\pm 0.030$ mm、表面粗さ $1.05\ \mu$ m 以下とする

因果効果を推定するデータではなく、多変量確率の挙動を確認するための架空データです。

n = 900
product = rng.choice(["製品A", "製品B", "製品C"], n, p=[0.45, 0.35, 0.20])
shift = rng.choice(["日勤", "夜勤"], n, p=[0.68, 0.32])
night = (shift == "夜勤").astype(float)
product_temp = pd.Series(product).map({"製品A": 0.0, "製品B": 1.3, "製品C": -0.8}).to_numpy()
product_pressure = pd.Series(product).map({"製品A": 0.0, "製品B": 2.2, "製品C": -1.5}).to_numpy()
product_dimension = pd.Series(product).map({"製品A": 0.000, "製品B": 0.004, "製品C": -0.003}).to_numpy()
product_roughness = pd.Series(product).map({"製品A": 0.00, "製品B": 0.07, "製品C": 0.12}).to_numpy()
product_energy = pd.Series(product).map({"製品A": 0.00, "製品B": 0.10, "製品C": 0.18}).to_numpy()

z1, z2, z3, z4, z5 = rng.normal(size=(5, n))
temp_std = z1
pressure_std = 0.25 * z1 + np.sqrt(1 - 0.25**2) * z2
dimension_std = 0.55 * z1 - 0.30 * z2 + np.sqrt(1 - 0.55**2 - 0.30**2) * z3
roughness_std = 0.40 * z1 - 0.25 * z2 + 0.55 * z3 + np.sqrt(1 - 0.40**2 - 0.25**2 - 0.55**2) * z4
energy_std = 0.65 * z1 + 0.20 * z2 + np.sqrt(1 - 0.65**2 - 0.20**2) * z5

mold_temp = 180 + product_temp + 1.4 * night + 4.2 * temp_std
injection_pressure = 92 + product_pressure + 1.0 * night + 5.5 * pressure_std
dimension_deviation = product_dimension + 0.002 * night + 0.018 * dimension_std
surface_roughness = product_roughness + 0.035 * night + 0.76 + 0.20 * roughness_std
energy_kwh = product_energy + 0.05 * night + 1.85 + 0.24 * energy_std

df = pd.DataFrame({
    "ロットID": [f"L{i:04d}" for i in range(1, n + 1)],
    "製品": product,
    "勤務帯": shift,
    "金型温度_C": mold_temp,
    "射出圧力_MPa": injection_pressure,
    "寸法偏差_mm": dimension_deviation,
    "表面粗さ_um": surface_roughness,
    "電力原単位_kWh": energy_kwh,
})
df["寸法状態"] = pd.cut(
    df["寸法偏差_mm"], [-np.inf, -0.030, 0.030, np.inf],
    labels=["下限外", "規格内", "上限外"]
)
df["粗さ状態"] = np.where(df["表面粗さ_um"] <= 1.05, "規格内", "規格外")
df["温度帯"] = pd.cut(
    df["金型温度_C"], [-np.inf, 178, 184, np.inf],
    labels=["低温", "標準", "高温"]
)
df["総合判定"] = np.where(
    (df["寸法状態"] == "規格内") & (df["粗さ状態"] == "規格内"), "合格", "要確認"
)

display(df.head(8).round(4))
display(
    df.groupby("製品", observed=True).agg(
        ロット数=("ロットID", "size"),
        平均温度_C=("金型温度_C", "mean"),
        平均寸法偏差_mm=("寸法偏差_mm", "mean"),
        平均粗さ_um=("表面粗さ_um", "mean"),
        要確認率=("総合判定", lambda s: (s == "要確認").mean()),
    ).round(4)
)
print(f"欠損数: {int(df.isna().sum().sum())} / 総セル数: {df.size:,}")

	ロットID	製品	勤務帯	金型温度_C	射出圧力_MPa	寸法偏差_mm	表面粗さ_um	電力原単位_kWh	寸法状態	粗さ状態	温度帯	総合判定
0	L0001	製品A	日勤	188.8566	93.2255	0.0261	1.0195	2.0346	規格内	規格内	高温	合格
1	L0002	製品C	日勤	175.9133	88.3551	-0.0191	0.8189	1.9327	規格内	規格内	低温	合格
2	L0003	製品B	夜勤	183.1683	95.6314	0.0003	0.6809	1.9350	規格内	規格内	標準	合格
3	L0004	製品B	日勤	179.5093	95.3241	-0.0115	0.5540	1.5872	規格内	規格内	標準	合格
4	L0005	製品C	日勤	174.6242	88.3919	-0.0006	0.6642	1.7433	規格内	規格内	低温	合格
5	L0006	製品A	夜勤	184.7163	90.1155	-0.0029	0.7349	2.2325	規格内	規格内	高温	合格
6	L0007	製品A	日勤	180.6574	92.1942	-0.0207	0.6606	2.2577	規格内	規格内	標準	合格
7	L0008	製品C	日勤	182.1593	87.6065	-0.0066	0.9197	2.0499	規格内	規格内	標準	合格

	ロット数	平均温度_C	平均寸法偏差_mm	平均粗さ_um	要確認率
製品
製品A	415	180.2865	0.0010	0.7668	0.1446
製品B	321	181.7583	0.0037	0.8398	0.1838
製品C	164	179.2849	-0.0015	0.8972	0.2683

欠損数: 0 / 総セル数: 10,800

No.061：同時分布

実務での意味

寸法と粗さを別々に見るだけでは、「寸法上限外かつ粗さ規格外」のような組合せリスクは分かりません。同時分布は、複数の品質状態が同時に起こる確率を表し、複合判定や選別能力の設計に使えます。

分析・モデル化の考え方

離散変数 $X,Y$ の同時確率質量関数は

p_{X,Y}(x,y)=P(X=x,Y=y), \qquad \sum_x\sum_y p_{X,Y}(x,y)=1

です。データでは、各組合せのロット数を全900ロットで割った経験同時確率を使います。0件の組合せも含め、確率表の合計が1になることを確認します。

Pythonで確認する

dimension_order = ["下限外", "規格内", "上限外"]
roughness_order = ["規格内", "規格外"]
joint_count = pd.crosstab(df["寸法状態"], df["粗さ状態"]).reindex(
    index=dimension_order, columns=roughness_order, fill_value=0
)
joint_prob = joint_count / len(df)
display(joint_count.rename_axis("寸法状態 / ロット数"))
display(joint_prob.rename_axis("寸法状態 / 同時確率").round(4))
print(f"同時確率の合計: {joint_prob.to_numpy().sum():.6f}")

fig, ax = plt.subplots(figsize=(7, 4.5))
im = ax.imshow(joint_prob.to_numpy(), cmap="Blues", vmin=0)
ax.set_xticks(range(len(roughness_order)), roughness_order)
ax.set_yticks(range(len(dimension_order)), dimension_order)
for i in range(len(dimension_order)):
    for j in range(len(roughness_order)):
        ax.text(j, i, f"{joint_prob.iloc[i, j]:.1%}", ha="center", va="center")
fig.colorbar(im, ax=ax, label="同時確率")
ax.set_title("寸法状態と粗さ状態の同時分布")
ax.set_xlabel("粗さ状態")
ax.set_ylabel("寸法状態")
ax.grid(False)
plt.tight_layout()
plt.show()

粗さ状態	規格内	規格外
寸法状態 / ロット数
下限外	35	0
規格内	737	78
上限外	16	34

粗さ状態	規格内	規格外
寸法状態 / 同時確率
下限外	0.0389	0.0000
規格内	0.8189	0.0867
上限外	0.0178	0.0378

同時確率の合計: 1.000000

png

結果の読み取り

確率表の1セルが、2つの状態の同時発生率です。規格内×規格内が最多でも、複合的な要確認セルを合計すると選別対象の規模が分かります。実務では、この表を製品・設備・期間で比較し、件数の少ないセルは率だけでなく実数も示します。規格判定の境界が変われば同時分布も変わるため、判定定義と版管理が必要です。

No.062：周辺分布

実務での意味

同時分布から相手の状態を合計して、寸法だけ、粗さだけの分布を取り出したものが周辺分布です。全体品質KPIと複合判定表が整合しているかを検算できます。

分析・モデル化の考え方

$Y$ を合計して得る $X$ の周辺分布は

p_X(x)=\sum_y p_{X,Y}(x,y)

であり、同様に $p_Y(y)=\sum_x p_{X,Y}(x,y)$ です。「周辺」は同時確率表の行端・列端に合計が置かれることに由来します。周辺分布だけからは、2変数がどう組み合わさるかは復元できません。

Pythonで確認する

margin_dimension = joint_prob.sum(axis=1).rename("同時表からの寸法周辺確率")
margin_roughness = joint_prob.sum(axis=0).rename("同時表からの粗さ周辺確率")
direct_dimension = df["寸法状態"].value_counts(normalize=True).reindex(dimension_order).rename("直接集計")
direct_roughness = df["粗さ状態"].value_counts(normalize=True).reindex(roughness_order).rename("直接集計")

display(pd.concat([margin_dimension, direct_dimension], axis=1).round(4))
display(pd.concat([margin_roughness, direct_roughness], axis=1).round(4))

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].bar(margin_dimension.index, margin_dimension.values, color="#4C78A8")
axes[0].set_title("寸法状態の周辺分布")
axes[0].set_xlabel("寸法状態")
axes[0].set_ylabel("確率")
axes[0].grid(axis="y", alpha=0.3)
axes[1].bar(margin_roughness.index, margin_roughness.values, color="#F58518")
axes[1].set_title("粗さ状態の周辺分布")
axes[1].set_xlabel("粗さ状態")
axes[1].set_ylabel("確率")
axes[1].grid(axis="y", alpha=0.3)
fig.suptitle("同時分布から取り出した2つの周辺分布")
plt.tight_layout()
plt.show()

	同時表からの寸法周辺確率	直接集計
寸法状態
下限外	0.0389	0.0389
規格内	0.9056	0.9056
上限外	0.0556	0.0556

	同時表からの粗さ周辺確率	直接集計
粗さ状態
規格内	0.8756	0.8756
規格外	0.1244	0.1244

png

結果の読み取り

同時確率表の行和・列和は、元データを単独集計した確率と一致します。ただし、たとえば粗さ規格外率という周辺KPIだけでは、そのロットが寸法規格内だったかは分かりません。経営報告には周辺KPIが簡潔ですが、選別・原因解析・複合保証には同時分布を残す、という使い分けが必要です。

No.063：条件付き分布

実務での意味

金型温度帯が観測された後に、粗さ状態の確率を更新します。全体不良率ではなく「高温なら何％か」が分かるため、温度アラーム後の抜取強化や条件調整の判断材料になります。

分析・モデル化の考え方

$P(X=x)>0$ のとき、条件付き確率は

P(Y=y\mid X=x)=\frac{P(X=x,Y=y)}{P(X=x)}

です。各温度帯の行内で合計が1になるよう正規化します。条件付き分布の違いは関連を示しますが、温度操作の因果効果を直接示すものではありません。製品構成や勤務帯が交絡している可能性があります。

Pythonで確認する

temp_order = ["低温", "標準", "高温"]
conditional_count = pd.crosstab(df["温度帯"], df["粗さ状態"]).reindex(
    index=temp_order, columns=roughness_order, fill_value=0
)
conditional_prob = conditional_count.div(conditional_count.sum(axis=1), axis=0)
display(conditional_count.rename_axis("温度帯 / ロット数"))
display(conditional_prob.rename_axis("温度帯 / 条件付き確率").round(4))
print("各温度帯での確率合計:")
display(conditional_prob.sum(axis=1).rename("合計").to_frame())

conditional_prob.plot(kind="bar", stacked=True, color=["#54A24B", "#E45756"], figsize=(8, 4.5))
plt.title("温度帯を条件とした粗さ状態の分布")
plt.xlabel("金型温度帯")
plt.ylabel("条件付き確率")
plt.xticks(rotation=0)
plt.grid(axis="y", alpha=0.3)
plt.legend(title="粗さ状態", bbox_to_anchor=(1.02, 1), loc="upper left")
plt.tight_layout()
plt.show()

粗さ状態	規格内	規格外
温度帯 / ロット数
低温	253	7
標準	382	52
高温	153	53

粗さ状態	規格内	規格外
温度帯 / 条件付き確率
低温	0.9731	0.0269
標準	0.8802	0.1198
高温	0.7427	0.2573

各温度帯での確率合計:

	合計
温度帯
低温	1.0000
標準	1.0000
高温	1.0000

png

結果の読み取り

温度帯ごとに粗さ規格外の比率が異なり、温度情報が品質確率の更新に使えることが分かります。アラーム設計では、規格外率だけでなく各帯のロット数、誤報コスト、見逃しコストを併記します。高温帯を狭くしすぎると該当件数が減るため、しきい値は物理知識と必要検出力に基づいて決めます。

No.064：独立性

実務での意味

高温と粗さ規格外が独立なら、同時発生率は各発生率の積で見積もれます。独立でなければ、単独リスクを掛け合わせるだけのFMEAやシミュレーションは同時異常を過小・過大評価します。

分析・モデル化の考え方

$X,Y$ が独立であるとは、すべての組合せで

P(X=x,Y=y)=P(X=x)P(Y=y)

が成り立つことです。今回は「高温」と「粗さ規格外」という2事象について、観測同時確率と独立仮定下の積を比較します。有限標本では完全一致しないため、差の大きさ、業務影響、期間安定性を見ます。相関0は一般には独立を保証しません。

Pythonで確認する

is_high_temp = df["温度帯"].eq("高温")
is_rough_ng = df["粗さ状態"].eq("規格外")
p_high = is_high_temp.mean()
p_ng = is_rough_ng.mean()
p_both = (is_high_temp & is_rough_ng).mean()
p_product = p_high * p_ng

independence_table = pd.DataFrame({
    "指標": ["P(高温)", "P(粗さ規格外)", "観測 P(高温∩粗さ規格外)", "独立仮定 P(高温)P(粗さ規格外)"],
    "確率": [p_high, p_ng, p_both, p_product],
})
display(independence_table.assign(確率_pct=lambda x: x["確率"].map(lambda v: f"{v:.2%}")))
print(f"同時確率比（観測 / 独立仮定）: {p_both / p_product:.2f} 倍")
print(f"追加同時発生ロット数の目安: {(p_both - p_product) * len(df):.1f} ロット / {len(df)}ロット")

plt.bar(["観測同時確率", "独立仮定の積"], [p_both, p_product], color=["#E45756", "#4C78A8"])
plt.title("高温と粗さ規格外：観測値と独立仮定の比較")
plt.xlabel("評価方法")
plt.ylabel("同時発生確率")
plt.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	指標	確率	確率_pct
0	P(高温)	0.2289	22.89%
1	P(粗さ規格外)	0.1244	12.44%
2	観測 P(高温∩粗さ規格外)	0.0589	5.89%
3	独立仮定 P(高温)P(粗さ規格外)	0.0285	2.85%

同時確率比（観測 / 独立仮定）: 2.07 倍
追加同時発生ロット数の目安: 27.4 ロット / 900ロット

png

結果の読み取り

観測同時確率が独立仮定の積を上回るなら、2事象を独立に生成するリスクモデルは同時発生を過小評価します。ただし、これは全ての値の組合せに対する独立性検証ではありません。実務導入ではクロス集計の期待度数、製品・勤務帯による層別、複数期間での再現性を確認し、必要に応じてカイ二乗検定や回帰モデルを補助的に使います。

No.065：共分散行列

実務での意味

センサーが増えると、2変数ずつの共分散を個別管理するのは困難です。共分散行列は、各KPIのばらつきとKPI間の共変動を1つの対称行列へまとめ、線形合成KPIのリスク計算や多変量監視の基礎になります。

分析・モデル化の考え方

確率ベクトル $\mathbf{X}=(X_1,\ldots,X_p)^\top$ の共分散行列は

\boldsymbol{\Sigma}=E[(\mathbf{X}-\boldsymbol{\mu})(\mathbf{X}-\boldsymbol{\mu})^\top]

です。対角成分は分散、非対角成分は共分散で、 $\Sigma_{ij}=\Sigma_{ji}$ です。任意の重み $\mathbf{a}$ に対して $\mathrm{Var}(\mathbf{a}^\top\mathbf{X})=\mathbf{a}^\top\boldsymbol{\Sigma}\mathbf{a}$ となります。

Pythonで確認する

sensor_cols = ["金型温度_C", "射出圧力_MPa", "寸法偏差_mm", "表面粗さ_um", "電力原単位_kWh"]
cov_matrix = df[sensor_cols].cov()
display(cov_matrix.round(6))
print(f"対称性の最大誤差: {np.abs(cov_matrix.to_numpy() - cov_matrix.to_numpy().T).max():.2e}")
print(f"最小固有値: {np.linalg.eigvalsh(cov_matrix.to_numpy()).min():.8f}（0以上なら半正定値）")

# 温度と圧力を標準化し、同じ重みで合成した工程負荷指数の分散を確認
two = df[["金型温度_C", "射出圧力_MPa"]]
two_z = (two - two.mean()) / two.std(ddof=1)
sigma_two = two_z.cov().to_numpy()
a = np.array([0.6, 0.4])
load_index = two_z.to_numpy() @ a
var_matrix = a @ sigma_two @ a
var_direct = load_index.var(ddof=1)
display(pd.DataFrame({"計算方法": ["a'Σa", "合成後に直接計算"], "分散": [var_matrix, var_direct]}).round(6))

fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(cov_matrix.to_numpy(), cmap="RdBu_r")
ax.set_xticks(range(len(sensor_cols)), sensor_cols, rotation=30, ha="right")
ax.set_yticks(range(len(sensor_cols)), sensor_cols)
fig.colorbar(im, ax=ax, label="共分散（単位に依存）")
ax.set_title("工程KPIの共分散行列")
ax.set_xlabel("KPI")
ax.set_ylabel("KPI")
ax.grid(False)
plt.tight_layout()
plt.show()

	金型温度_C	射出圧力_MPa	寸法偏差_mm	表面粗さ_um	電力原単位_kWh
金型温度_C	19.8177	6.5382	0.0462	0.3915	0.6696
射出圧力_MPa	6.5382	30.4943	-0.0136	-0.1519	0.4433
寸法偏差_mm	0.0462	-0.0136	0.0003	0.0027	0.0013
表面粗さ_um	0.3915	-0.1519	0.0027	0.0421	0.0151
電力原単位_kWh	0.6696	0.4433	0.0013	0.0151	0.0606

対称性の最大誤差: 0.00e+00
最小固有値: 0.00012796（0以上なら半正定値）

	計算方法	分散
0	a'Σa	0.6477
1	合成後に直接計算	0.6477

png

結果の読み取り

共分散行列は対称で、対角には各指標の分散が並びます。また、行列式で求めた合成指数の分散と直接計算が一致します。一方、mmとMPaでは単位・尺度が大きく異なるため、色の大きさだけで関係の強さを比較できません。原単位のまま損失分散を計算する用途には共分散、関係の比較には次の相関行列を使います。

No.066：相関行列

実務での意味

相関行列は各KPIを標準化し、線形な関係の方向と強さを $-1$ から $1$ で比較します。似た情報を持つセンサー、品質と連動する候補条件、多重共線性の可能性を探索する入口です。

分析・モデル化の考え方

標準偏差を並べた対角行列を $\mathbf{D}$ とすると、相関行列は

\mathbf{R}=\mathbf{D}^{-1}\boldsymbol{\Sigma}\mathbf{D}^{-1}

です。対角成分は1で、単位換算の影響を受けません。ただしPearson相関は線形関係の指標であり、非線形関係、外れ値、群混在には注意が必要です。相関は因果を意味しません。

Pythonで確認する

std = df[sensor_cols].std(ddof=1).to_numpy()
d_inv = np.diag(1 / std)
corr_from_cov = d_inv @ cov_matrix.to_numpy() @ d_inv
corr_direct = df[sensor_cols].corr().to_numpy()
correlation = pd.DataFrame(corr_direct, index=sensor_cols, columns=sensor_cols)
display(correlation.round(3))
print(f"共分散行列からの計算との最大差: {np.abs(corr_from_cov - corr_direct).max():.2e}")

fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(corr_direct, vmin=-1, vmax=1, cmap="RdBu_r")
ax.set_xticks(range(len(sensor_cols)), sensor_cols, rotation=30, ha="right")
ax.set_yticks(range(len(sensor_cols)), sensor_cols)
for i in range(len(sensor_cols)):
    for j in range(len(sensor_cols)):
        ax.text(j, i, f"{corr_direct[i, j]:.2f}", ha="center", va="center",
                color="white" if abs(corr_direct[i, j]) > 0.55 else "black")
fig.colorbar(im, ax=ax, label="Pearson相関係数")
ax.set_title("工程KPIの相関行列")
ax.set_xlabel("KPI")
ax.set_ylabel("KPI")
ax.grid(False)
plt.tight_layout()
plt.show()

	金型温度_C	射出圧力_MPa	寸法偏差_mm	表面粗さ_um	電力原単位_kWh
金型温度_C	1.0000	0.2660	0.5650	0.4290	0.6110
射出圧力_MPa	0.2660	1.0000	-0.1340	-0.1340	0.3260
寸法偏差_mm	0.5650	-0.1340	1.0000	0.7150	0.2850
表面粗さ_um	0.4290	-0.1340	0.7150	1.0000	0.2980
電力原単位_kWh	0.6110	0.3260	0.2850	0.2980	1.0000

共分散行列からの計算との最大差: 1.94e-15

png

結果の読み取り

寸法偏差・粗さ・温度などの正負の連動を、単位によらず比較できます。強い相関を持つセンサーは冗長化や原因候補の手掛かりになりますが、片方を削除してよいとは限りません。制御応答速度、校正方法、故障モードも考慮します。製品別相関と全体相関が逆転する場合もあるため、全体行列だけで設備条件を変更しません。

No.067：多変量正規分布

実務での意味

個々のKPIが管理限界内でも、組合せとして通常領域から外れることがあります。多変量正規分布は、平均ベクトルと共分散行列で楕円状の通常領域を表し、複数センサーの同時監視や異常候補抽出に利用できます。

分析・モデル化の考え方

$p$ 次元多変量正規分布の密度は

f(\mathbf{x})=\frac{1}{(2\pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}} \exp\left[-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right]

です。指数部の二次形式はマハラノビス距離の二乗 $D_M^2$ です。正規仮定の下で $D_M^2$ は概ね自由度 $p$ のカイ二乗分布に従うため、確率楕円を作れます。ここでは混合を避け、製品A・日勤に限定します。

Pythonで確認する

subset = df.query("製品 == '製品A' and 勤務帯 == '日勤'")[["寸法偏差_mm", "表面粗さ_um"]].copy()
x2 = subset.to_numpy()
mu2 = x2.mean(axis=0)
sigma2 = np.cov(x2, rowvar=False, ddof=1)
centered = x2 - mu2
d2 = np.einsum("ni,ij,nj->n", centered, np.linalg.inv(sigma2), centered)
threshold95 = chi2.ppf(0.95, df=2)
outside = d2 > threshold95

eigvals, eigvecs = np.linalg.eigh(sigma2)
angles = np.linspace(0, 2 * np.pi, 240)
unit_circle = np.vstack([np.cos(angles), np.sin(angles)])
ellipse = mu2[:, None] + eigvecs @ np.diag(np.sqrt(eigvals * threshold95)) @ unit_circle

print(f"対象ロット数: {len(subset)}")
print(f"95%楕円外: {outside.sum()}ロット（{outside.mean():.2%}）")
display(pd.DataFrame({"平均": mu2}, index=subset.columns).round(5))

plt.scatter(x2[~outside, 0], x2[~outside, 1], s=24, alpha=0.45, color="#4C78A8", label="95%楕円内")
plt.scatter(x2[outside, 0], x2[outside, 1], s=34, alpha=0.85, color="#E45756", label="95%楕円外")
plt.plot(ellipse[0], ellipse[1], color="#222222", linewidth=2, label="多変量正規の95%楕円")
plt.title("寸法偏差と表面粗さの多変量通常領域（製品A・日勤）")
plt.xlabel("寸法偏差（mm）")
plt.ylabel("表面粗さ（μm）")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

対象ロット数: 287
95%楕円外: 17ロット（5.92%）

	平均
寸法偏差_mm	0.0003
表面粗さ_um	0.7521

png

結果の読み取り

楕円の傾きが2変数の相関方向、幅がばらつきを表します。各軸では極端でなくても、相関方向から外れた点はマハラノビス距離が大きくなります。ただし楕円外は直ちに不良ではなく、調査優先候補です。実運用前にはQ-Qプロット等で分布形状を確認し、製品切替、非線形境界、多峰性、時系列自己相関が強い場合は別モデルを検討します。

No.068：条件付き正規分布

実務での意味

当日の金型温度が分かった後、電力原単位の予測平均と予測幅を更新します。基準電力を1本の固定値にせず、観測可能な条件に応じたレンジとして持つことで、過剰な警報と見逃しを減らせます。

分析・モデル化の考え方

$(X,Y)$ が2変量正規分布に従うとき、 $X=x$ の下での $Y$ は正規分布となり、

E[Y\mid X=x]=\mu_Y+\frac{\sigma_{YX}}{\sigma_X^2}(x-\mu_X)

\mathrm{Var}(Y\mid X=x)=\sigma_Y^2-\frac{\sigma_{YX}^2}{\sigma_X^2}

です。相関があるほど、 $X$ の観測により $Y$ の不確実性が小さくなります。ここでも単一条件に近づけるため製品A・日勤を使います。

Pythonで確認する

cond_df = df.query("製品 == '製品A' and 勤務帯 == '日勤'")[["金型温度_C", "電力原単位_kWh"]]
x = cond_df["金型温度_C"].to_numpy()
y = cond_df["電力原単位_kWh"].to_numpy()
mu_x, mu_y = x.mean(), y.mean()
cov_xy = np.cov(x, y, ddof=1)
var_x, var_y = cov_xy[0, 0], cov_xy[1, 1]
cov_yx = cov_xy[1, 0]
beta = cov_yx / var_x
conditional_sd = np.sqrt(var_y - cov_yx**2 / var_x)

x_grid = np.linspace(x.min(), x.max(), 200)
conditional_mean = mu_y + beta * (x_grid - mu_x)
lower = conditional_mean - 1.96 * conditional_sd
upper = conditional_mean + 1.96 * conditional_sd
target_temp = 186.0
target_mean = mu_y + beta * (target_temp - mu_x)

comparison = pd.DataFrame({
    "予測": ["温度を知らない場合", f"温度={target_temp:.1f}℃を知った場合"],
    "平均_kWh": [mu_y, target_mean],
    "標準偏差_kWh": [np.sqrt(var_y), conditional_sd],
})
display(comparison.round(4))

plt.scatter(x, y, s=22, alpha=0.35, color="#4C78A8", label="各ロット")
plt.plot(x_grid, conditional_mean, color="#E45756", linewidth=2, label="条件付き平均")
plt.fill_between(x_grid, lower, upper, color="#E45756", alpha=0.18, label="条件付き95%範囲")
plt.axvline(target_temp, color="#222222", linestyle="--", label=f"評価温度 {target_temp:.0f}℃")
plt.title("金型温度を観測した後の電力原単位分布")
plt.xlabel("金型温度（℃）")
plt.ylabel("電力原単位（kWh/ロット）")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	予測	平均_kWh	標準偏差_kWh
0	温度を知らない場合	1.8466	0.2370
1	温度=186.0℃を知った場合	2.0652	0.1754

png

結果の読み取り

温度が平均より高いと条件付き平均電力も上がり、温度を知った後の標準偏差は無条件の標準偏差より小さくなります。これにより温度補正済みの監視帯を作れます。ただし、範囲はモデル上の条件付き分布であり平均推定の信頼区間ではありません。外挿、製品混在、分散の温度依存、設備劣化があると式の前提が崩れます。

No.069：ガウス過程

実務での意味

設備のゼロ点ドリフトは時間に対して滑らかに変化しますが、校正測定は連続ではありません。ガウス過程は、少数の測定点からドリフト曲線を補間し、データが少ない区間ほど広い不確実性を示せます。校正時期や追加測定点の判断に向きます。

分析・モデル化の考え方

ガウス過程は関数 $f(t)$ に対する確率分布で、

f(t)\sim\mathcal{GP}(m(t),k(t,t'))

と書きます。ここでは平均関数を0、カーネルをRBF

k(t,t')=\sigma_f^2\exp\left(-\frac{(t-t')^2}{2\ell^2}\right)

とし、観測ノイズ分散 $\sigma_n^2$ を加えます。予測平均と分散はカーネル行列の条件付き正規分布から得られます。 $\ell$ は変化の滑らかさを表します。

Pythonで確認する

gp_rng = np.random.default_rng(SEED + 69)
t_obs = np.array([0, 4, 9, 14, 20, 27, 34, 42, 51, 60], dtype=float)
true_drift = lambda t: 0.004 * np.sin(t / 10) + 0.00012 * t
y_obs = true_drift(t_obs) + gp_rng.normal(0, 0.0012, len(t_obs))
t_pred = np.linspace(0, 70, 281)

length_scale = 13.0
signal_sd = 0.006
noise_sd = 0.0012

def rbf_kernel(a, b, length_scale, signal_sd):
    sq_distance = (np.asarray(a)[:, None] - np.asarray(b)[None, :]) ** 2
    return signal_sd**2 * np.exp(-0.5 * sq_distance / length_scale**2)

k_oo = rbf_kernel(t_obs, t_obs, length_scale, signal_sd) + noise_sd**2 * np.eye(len(t_obs))
k_po = rbf_kernel(t_pred, t_obs, length_scale, signal_sd)
k_pp = rbf_kernel(t_pred, t_pred, length_scale, signal_sd)
alpha = np.linalg.solve(k_oo, y_obs)
gp_mean = k_po @ alpha
v = np.linalg.solve(k_oo, k_po.T)
gp_cov = k_pp - k_po @ v
gp_sd = np.sqrt(np.clip(np.diag(gp_cov), 0, None))

gp_summary = pd.DataFrame({
    "地点": ["最終観測日", "10日先"],
    "日": [60, 70],
    "予測平均_mm": [np.interp(60, t_pred, gp_mean), np.interp(70, t_pred, gp_mean)],
    "予測標準偏差_mm": [np.interp(60, t_pred, gp_sd), np.interp(70, t_pred, gp_sd)],
})
display(gp_summary.round(5))

plt.scatter(t_obs, y_obs, color="#222222", s=42, zorder=3, label="校正測定")
plt.plot(t_pred, gp_mean, color="#4C78A8", linewidth=2, label="GP予測平均")
plt.fill_between(t_pred, gp_mean - 1.96 * gp_sd, gp_mean + 1.96 * gp_sd,
                 color="#4C78A8", alpha=0.2, label="潜在関数の95%区間")
plt.axvline(t_obs.max(), color="#E45756", linestyle="--", label="最終観測日")
plt.title("ガウス過程による設備ゼロ点ドリフトの補間・予測")
plt.xlabel("稼働日")
plt.ylabel("ゼロ点ドリフト（mm）")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	地点	日	予測平均_mm	予測標準偏差_mm
0	最終観測日	60	0.0041	0.0011
1	10日先	70	0.0029	0.0035

png

結果の読み取り

観測点の近くでは区間が狭く、最終観測日より先では徐々に広がります。点予測だけでなく不確実性が広がる位置を追加校正の候補にできます。今回は説明のためハイパーパラメータを固定しました。実運用では、周辺尤度による推定、設備停止・交換による変化点、季節性、予測時点より後のデータを学習に混ぜない時系列検証が必要です。

No.070：コピュラ

実務での意味

設備停止時間と納期損失額は、ともに右裾が長く、同じ日に大きくなりやすい指標です。コピュラは各指標の周辺分布と依存構造を分離し、「分布形は実績に合わせたまま、同時悪化の強さを変える」ストレステストに利用できます。

分析・モデル化の考え方

Sklarの定理により、連続な周辺分布 $F_X,F_Y$ を持つ同時分布は

F_{X,Y}(x,y)=C(F_X(x),F_Y(y))

と表せます。 $C$ がコピュラです。今回は相関 $\rho$ の2変量正規乱数 $(Z_1,Z_2)$ を $U_i=\Phi(Z_i)$ で一様分布へ写し、右裾の長い周辺分布へ変換するガウスコピュラを使います。比較対象は同じ周辺分布を持つ独立モデルです。

Pythonで確認する

copula_rng = np.random.default_rng(SEED + 70)
n_scenarios = 50_000
rho = 0.72
z_dep = copula_rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]], size=n_scenarios)
z_ind = copula_rng.normal(size=(n_scenarios, 2))

# Phi(Z)を周辺分布の逆関数へ通すのと同値な対数正規変換
downtime_dep = np.exp(1.15 + 0.65 * z_dep[:, 0])
loss_dep = np.exp(4.80 + 0.90 * z_dep[:, 1])  # 万円
downtime_ind = np.exp(1.15 + 0.65 * z_ind[:, 0])
loss_ind = np.exp(4.80 + 0.90 * z_ind[:, 1])

downtime_q90 = np.quantile(downtime_dep, 0.90)
loss_q90 = np.quantile(loss_dep, 0.90)
joint_tail_dep = np.mean((downtime_dep > downtime_q90) & (loss_dep > loss_q90))
joint_tail_ind = np.mean((downtime_ind > downtime_q90) & (loss_ind > loss_q90))

copula_summary = pd.DataFrame({
    "モデル": ["ガウスコピュラ（依存あり）", "独立"],
    "停止時間中央値_h": [np.median(downtime_dep), np.median(downtime_ind)],
    "損失中央値_万円": [np.median(loss_dep), np.median(loss_ind)],
    "両方が各90%点超の確率": [joint_tail_dep, joint_tail_ind],
    "50,000回中の同時超過数": [joint_tail_dep * n_scenarios, joint_tail_ind * n_scenarios],
})
display(copula_summary.round(4))
print(f"90%点: 停止時間 {downtime_q90:.2f}時間、損失 {loss_q90:.1f}万円")
print(f"同時テール確率比: {joint_tail_dep / joint_tail_ind:.2f}倍")

sample_idx = copula_rng.choice(n_scenarios, 1800, replace=False)
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharex=True, sharey=True)
for ax, dx, ly, title in [
    (axes[0], downtime_dep, loss_dep, "依存あり（ガウスコピュラ）"),
    (axes[1], downtime_ind, loss_ind, "独立（同じ周辺分布）"),
]:
    ax.scatter(dx[sample_idx], ly[sample_idx], s=10, alpha=0.25, color="#4C78A8")
    ax.axvline(downtime_q90, color="#E45756", linestyle="--")
    ax.axhline(loss_q90, color="#E45756", linestyle="--")
    ax.set_title(title)
    ax.set_xlabel("停止時間（h）")
    ax.set_ylabel("納期損失（万円）")
    ax.set_xlim(0, np.quantile(downtime_dep, 0.995))
    ax.set_ylim(0, np.quantile(loss_dep, 0.995))
    ax.grid(alpha=0.3)
fig.suptitle("周辺分布が同じでも依存構造で同時テールが変わる")
plt.tight_layout()
plt.show()

	モデル	停止時間中央値_h	損失中央値_万円	両方が各90%点超の確率	50,000回中の同時超過数
0	ガウスコピュラ（依存あり）	3.1421	121.1719	0.0487	2,435.0000
1	独立	3.1645	122.0452	0.0092	460.0000

90%点: 停止時間 7.27時間、損失 384.6万円
同時テール確率比: 5.29倍

png

結果の読み取り

2モデルの停止時間・損失額の中央値はほぼ同じですが、両方が90%点を超える確率は依存ありモデルで大きくなります。周辺分布だけを個別に合わせても、同時リスクは再現できないことが分かります。なお、ガウスコピュラは極端な裾での依存を弱く見積もる場合があります。BCP用途ではtコピュラ等との比較、極端事象の定義、パラメータ推定誤差、ストレスシナリオを検討します。

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

単独KPIと組合せKPIを役割分担する
周辺分布は経営報告に簡潔ですが、複合品質・同時停止には同時分布が必要です。
観測情報で確率を更新する
温度帯や製品が分かった後は、全体確率ではなく条件付き分布を使います。条件別サンプル数も併記します。
独立仮定を便宜的に置いたままにしない
リスクの積み上げやモンテカルロ計画では、依存を無視すると同時テールを誤ります。
共分散と相関を使い分ける
原単位で合成損失の分散を計算するなら共分散、単位を越えて関係を比較するなら相関です。
モデルの通常領域と品質規格を区別する
多変量正規の楕円外は統計的に珍しい点であり、不良判定そのものではありません。
予測値には不確実性を付ける
条件付き正規分布やガウス過程は、平均だけでなく予測幅を意思決定へ渡せます。
周辺分布と依存構造を別々に検証する
コピュラでは各KPIの分布適合と同時発生構造の適合を分けて評価します。

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

1. 分析単位と時刻をそろえる

ロット、個体、設備秒次など粒度を統一し、材料投入、条件設定、検査、停止、出荷を追跡可能なキーで結びます。センサー時計のずれと集計窓も管理します。

2. KPI・規格・欠測の定義を版管理する

測定器、単位、規格上下限、丸め、再測定、欠測補完、設備交換を記録します。規格変更前後の確率を同じ母集団として混ぜません。

3. 層別と時系列構造を確認する

製品、設備、金型、材料、勤務帯で分布が異なる可能性があります。ランダム分割だけでなく、将来期間を使った検証とドリフト監視を行います。

4. 前提診断と代替モデルを用意する

散布図、外れ値、正規性、多峰性、非線形、自己相関を確認します。前提が合わなければロバスト共分散、混合分布、木モデル、時系列モデル、別コピュラを比較します。

5. 判断ルールを損失で設計する

統計的に珍しいことと、事業上重大なことは同じではありません。見逃し、誤報、停止、追加検査、保証費の損失を定義し、しきい値と対応フローを合意します。

6. 小規模な並行運用から始める

既存判定をすぐ置換せず、特定設備で並行表示し、誰がいつ確認し何をするかを標準作業へ落とします。再学習・再校正・監査ログの責任者も決めます。

まとめ

No.061〜No.070では、2変数の確率表から、行列、正規分布、関数分布、コピュラまでを一貫した製造業の題材で確認しました。

同時分布は組合せ、周辺分布は単独、条件付き分布は情報取得後の確率を表す
独立なら同時確率は積になるが、相関0だけでは一般に独立とはいえない
共分散行列は原単位での共変動、相関行列は標準化した線形関係を表す
多変量正規分布は楕円状の通常領域、条件付き正規分布は観測後の更新に使える
ガウス過程は関数の予測と不確実性、コピュラは周辺分布と依存構造の分離を扱う

実務で重要なのは、精巧なモデル名ではなく、分析単位、条件、前提、不確実性、対応行動を一続きに設計することです。まずは重要な2〜3個のKPIから同時分布を可視化し、現場の因果仮説と照合するところから始められます。

法人向けのご相談

数理工房では、製造業向けのデータ分析・統計モデリング・異常検知・リスクシミュレーションについて、課題整理からPoC、実装、現場研修までご支援しています。

品質KPI・設備センサーの多変量監視設計
条件付き予測や時系列ドリフト検知のPoC
停止・納期・在庫の同時リスクシミュレーション
自社データを用いたPython・統計研修
分析結果を現場標準と意思決定プロセスへ組み込む支援

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

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

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

概要

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.061：同時分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.062：周辺分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.063：条件付き分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.064：独立性

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.065：共分散行列

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.066：相関行列

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.067：多変量正規分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.068：条件付き正規分布

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.069：ガウス過程

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.070：コピュラ

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

1. 分析単位と時刻をそろえる

2. KPI・規格・欠測の定義を版管理する

3. 層別と時系列構造を確認する

4. 前提診断と代替モデルを用意する

5. 判断ルールを損失で設計する

6. 小規模な並行運用から始める

まとめ

法人向けのご相談

Related Posts

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

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

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

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

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