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

製造業の回帰分析をPythonで実践｜線形回帰・GLM・GAM 10本ノック

2020-05-15

工程条件から品質・停止・工数を予測する：製造業の回帰モデル10本ノック

製造現場では、温度、加工速度、振動、材料状態、保全経過日数などが、エネルギー原単位、不適合、微停止、復旧工数に同時に影響します。本記事では、架空の精密部品工場の生産実績を題材に、線形回帰、Ridge、Lasso、ElasticNet、ロジスティック回帰、ポアソン回帰、負の二項回帰、GLM、GAM、分位点回帰を実装します。

目的はモデル名を覚えることではなく、「どの条件を動かすか」「不適合リスクをどこで警戒するか」「平均より悪いケースに何人時を見込むか」を、目的変数の性質に合うモデルで判断することです。対象は確率・統計Python実装100本ノックの No.061〜No.070 です。

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

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

架空の精密部品工場では、各生産バッチについて工程温度、加工速度、設備振動、材料含水率、保全後経過日数を記録しています。管理者は翌月の条件標準と要員計画を決めるため、次の問いに答える必要があります。

電力原単位を押し上げる条件は何か
多数の似たセンサー指標から、説明に必要な変数をどう絞るか
不適合の「有無」と微停止の「回数」を同じ回帰で扱ってよいか
微停止のばらつきが平均より大きいとき、発生リスクをどう表すか
平均復旧時間だけでなく、長引くケースをどう見積もるか

回帰モデルは、結果と要因の関連を定量化し、条件変更時の予測を行う道具です。ただし、目的変数が連続量、0/1、件数、正の歪んだ値のどれかによって、適切な分布とリンク関数は異なります。

現場でよくある状況

単回帰の係数を因果効果と解釈し、交絡や操業範囲を確認していない
相関の強いセンサーを同時投入し、係数の符号が会議のたびに変わる
不適合率を通常の線形回帰で予測し、0%未満や100%超の値が出る
停止件数を正規分布で扱い、ゼロ以上の整数という性質を無視する
平均予測の精度だけでモデルを選び、繁忙時の上振れを見落とす

モデル選択では、予測精度だけでなく、データの発生機構、係数の安定性、説明目的か予測目的か、意思決定時の損失が左右対称かを明確にします。

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

回帰の一般形を、平均構造として

g\{E(Y\mid X)\}=\beta_0+\beta_1x_1+\cdots+\beta_px_p

と書きます。 $g$ はリンク関数です。連続量の線形回帰では恒等リンク、二値データではロジット、件数では対数リンクが典型です。同じ説明変数でも、 $Y$ の尺度に合わせて分布とリンクを変える必要があります。

また、予測で見つかるのは原則として関連であり、因果ではありません。温度係数が正でも「温度を上げれば必ず悪化する」とは限りません。品種、設備、作業者、時間順序、条件変更のルールを考慮し、介入効果を知りたい場合は実験計画や因果推論を別途設計します。

今回扱うノックの全体像

No.	テーマ	製造業での問い
061	線形回帰	工程条件から電力原単位を説明できるか
062	Ridge回帰	相関するセンサーがあっても係数を安定化できるか
063	Lasso回帰	影響の小さい変数を候補から絞れるか
064	ElasticNet	変数選択と相関変数の保持を両立できるか
065	ロジスティック回帰	バッチ不適合確率を0〜1で推定できるか
066	ポアソン回帰	稼働時間を考慮して微停止件数を比較できるか
067	負の二項回帰	件数の過分散を織り込めるか
068	GLM	歪んだ復旧工数を分布とリンクで表現できるか
069	GAM	温度などの非線形な影響を滑らかに捉えられるか
070	Quantile Regression	復旧時間の中央値と上側リスクを分けて予測できるか

No.061〜064は連続量と正則化、No.065〜068は目的変数に応じた確率分布、No.069〜070は平均・直線だけでは表せない構造を扱います。

Python 環境の準備

NumPyで架空データを生成し、pandasで表を作ります。scikit-learnは前処理・正則化・予測評価、statsmodelsは統計モデル、matplotlibは可視化に使います。japanize_matplotlib は日本語表示のためだけに使い、seabornや外部データは使いません。乱数seedを固定し、同じ結果を再現できるようにします。

import platform
import numpy as np
import pandas as pd
import scipy
import sklearn
from sklearn.compose import TransformedTargetRegressor
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet, LogisticRegression
from sklearn.metrics import mean_squared_error, roc_auc_score, log_loss
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, SplineTransformer
import statsmodels.api as sm
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib

SEED = 20260711
rng = np.random.default_rng(SEED)
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"scikit-learn: {sklearn.__version__}")
print(f"statsmodels : {sm.__version__}")
print(f"Matplotlib  : {matplotlib.__version__}")
print(f"random seed : {SEED}")

Python      : 3.13.1
NumPy       : 2.5.1
pandas      : 3.0.3
SciPy       : 1.18.0
scikit-learn: 1.9.0
statsmodels : 0.14.6
Matplotlib  : 3.11.0
random seed : 20260711

架空データの作成

420バッチ分の生産実績を生成します。cooling_flow は工程温度と強く相関する冷却流量、sensor_temp_2 は温度の冗長な別センサーです。これらは正則化が必要になる状況を再現します。

微停止件数は、バッチ固有の見えない難しさによってばらつくようガンマ・ポアソン混合で生成します。そのため平均より分散が大きく、後でポアソン回帰と負の二項回帰の違いを確認できます。モデル学習では生成式や真の係数を使いません。

n = 420
temperature = rng.normal(72, 4.2, n)
speed = rng.normal(118, 13, n)
vibration = np.clip(rng.gamma(3.2, 0.32, n), 0.15, None)
moisture = rng.normal(0.82, 0.16, n)
maintenance_days = rng.integers(1, 121, n)
operating_hours = rng.uniform(5.5, 10.0, n)
cooling_flow = 34 - 0.32 * temperature + rng.normal(0, 0.9, n)
sensor_temp_2 = temperature + rng.normal(0, 0.7, n)
ambient_humidity = rng.normal(56, 9, n)

energy = (31 + 0.34 * speed + 0.75 * (temperature - 72)
          + 2.8 * vibration + rng.normal(0, 3.2, n))
defect_logit = (-3.1 + 0.15 * (temperature - 72) + 0.62 * (vibration - 1)
                + 0.018 * maintenance_days + 0.55 * (moisture - 0.82))
defect = rng.binomial(1, 1 / (1 + np.exp(-defect_logit)))
stop_rate = np.exp(-2.0 + 0.42 * vibration + 0.009 * maintenance_days)
latent_rate = rng.gamma(shape=1.5, scale=stop_rate / 1.5)
micro_stops = rng.poisson(operating_hours * latent_rate)
recovery_hours = rng.gamma(
    shape=2.5,
    scale=np.exp(-1.25 + 0.012 * maintenance_days + 0.27 * vibration) / 2.5,
)
cycle_time = (41 + 0.15 * speed + 4.6 * vibration + 0.020 * maintenance_days
              + rng.normal(0, 1.5 + 1.3 * vibration, n))
quality_loss = (1.8 + 0.12 * (temperature - 72) ** 2
                + 0.015 * (speed - 118) ** 2 + 1.7 * vibration + rng.normal(0, 1.2, n))

df = pd.DataFrame({
    "temperature": temperature, "sensor_temp_2": sensor_temp_2, "speed": speed,
    "vibration": vibration, "moisture": moisture, "maintenance_days": maintenance_days,
    "operating_hours": operating_hours, "cooling_flow": cooling_flow,
    "ambient_humidity": ambient_humidity, "energy": energy, "defect": defect,
    "micro_stops": micro_stops, "recovery_hours": recovery_hours,
    "cycle_time": cycle_time, "quality_loss": quality_loss,
})

display(df.head().round(3))
display(pd.DataFrame({
    "指標": ["バッチ数", "不適合率", "微停止の平均", "微停止の分散", "平均復旧工数"],
    "値": [len(df), df.defect.mean(), df.micro_stops.mean(),
           df.micro_stops.var(), df.recovery_hours.mean()],
}).round(3))

	temperature	sensor_temp_2	speed	vibration	moisture	maintenance_days	operating_hours	cooling_flow	ambient_humidity	energy	defect	micro_stops	recovery_hours	cycle_time	quality_loss
0	66.740	66.167	126.618	1.062	0.989	10	7.704	11.849	61.413	71.472	0	2	0.442	60.080	7.476
1	70.616	70.193	135.361	0.377	1.015	89	9.550	11.421	51.120	75.155	0	3	0.609	64.640	7.523
2	73.686	73.601	107.679	1.341	0.570	43	6.338	10.085	76.765	73.063	0	0	0.488	65.133	6.173
3	69.063	68.940	124.009	0.150	0.690	120	8.677	12.773	41.372	73.693	1	2	3.237	62.716	4.116
4	66.060	65.541	108.298	1.641	1.132	68	6.988	12.005	44.895	65.630	0	2	1.241	70.649	10.696

	指標	値
0	バッチ数	420.000
1	不適合率	0.143
2	微停止の平均	3.133
3	微停止の分散	10.946
4	平均復旧工数	0.811

No.061：線形回帰 — 工程条件から電力原単位を説明する

実務での意味

電力原単位と工程条件の関係を数値化すると、設定変更による省エネ余地の仮説を作れます。線形回帰は説明しやすく、基準モデルとして有用です。

分析・モデル化の考え方

線形回帰は

y_i=\beta_0+\sum_{j=1}^{p}\beta_jx_{ij}+\varepsilon_i

を置き、残差平方和 $\sum_i(y_i-\hat y_i)^2$ を最小化します。係数 $\beta_j$ は「他の変数を一定としたときの $x_j$ 1単位増加に対する平均変化」です。学習用と評価用を分け、未知データへのRMSEで確認します。

Pythonで確認する

reg_features = ["temperature", "speed", "vibration", "maintenance_days"]
X_train, X_test, y_train, y_test = train_test_split(
    df[reg_features], df["energy"], test_size=0.25, random_state=SEED
)
ols = LinearRegression().fit(X_train, y_train)
ols_pred = ols.predict(X_test)

ols_result = pd.DataFrame({"説明変数": reg_features, "係数": ols.coef_})
display(ols_result.round(3))
print(f"切片      : {ols.intercept_:.3f}")
print(f"評価RMSE  : {mean_squared_error(y_test, ols_pred)**0.5:.3f}")

plt.scatter(ols_pred, y_test, alpha=0.65)
limits = [min(ols_pred.min(), y_test.min()), max(ols_pred.max(), y_test.max())]
plt.plot(limits, limits, "--", color="tab:red", label="予測＝実測")
plt.title("線形回帰による電力原単位の予測")
plt.xlabel("予測電力原単位 [kWh/バッチ]")
plt.ylabel("実測電力原単位 [kWh/バッチ]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	説明変数	係数
0	temperature	0.695
1	speed	0.353
2	vibration	2.979
3	maintenance_days	0.005

切片      : -21.146
評価RMSE  : 2.856

png

結果の読み取り

評価RMSEは、未知バッチに対する典型的な予測誤差を元の単位で示します。係数の大きさは単位に依存するため、温度と保全日数の数値をそのまま重要度比較には使えません。また、予測と実測の関係が直線から系統的に外れるなら、非線形項や交互作用、未観測要因を検討します。係数は操業範囲内の条件付き関連であり、設定変更の因果効果を保証しません。

No.062：Ridge回帰 — 相関するセンサーの係数を安定化する

実務での意味

同じ現象を測るセンサーや、制御上連動する設定値を同時に使うと、多重共線性により係数が不安定になります。予測値が似ていても、係数の符号や大きさが標本の違いで変われば、改善施策の説明に使いにくくなります。

分析・モデル化の考え方

Ridge回帰は標準化後の係数にL2罰則を加え、

\min_{\beta}\left\{\sum_i(y_i-\beta_0-x_i^T\beta)^2+\alpha\sum_j\beta_j^2\right\}

を解きます。 $\alpha$ が大きいほど係数を0方向へ縮小しますが、通常は厳密な0にはしません。縮小によるバイアスと分散低下の交換で、予測と係数の安定化を狙います。

Pythonで確認する

corr_features = ["temperature", "sensor_temp_2", "cooling_flow", "speed", "vibration"]
Xr_train, Xr_test, yr_train, yr_test = train_test_split(
    df[corr_features], df["energy"], test_size=0.25, random_state=SEED
)
scaler = StandardScaler().fit(Xr_train)
Xr_train_z, Xr_test_z = scaler.transform(Xr_train), scaler.transform(Xr_test)

ols_z = LinearRegression().fit(Xr_train_z, yr_train)
ridge = Ridge(alpha=12).fit(Xr_train_z, yr_train)
ridge_comparison = pd.DataFrame({
    "説明変数": corr_features, "OLS標準化係数": ols_z.coef_, "Ridge標準化係数": ridge.coef_,
})
display(ridge_comparison.round(3))
print(f"OLS評価RMSE   : {mean_squared_error(yr_test, ols_z.predict(Xr_test_z))**0.5:.3f}")
print(f"Ridge評価RMSE : {mean_squared_error(yr_test, ridge.predict(Xr_test_z))**0.5:.3f}")

x = np.arange(len(corr_features))
plt.bar(x - 0.18, ols_z.coef_, width=0.36, label="OLS")
plt.bar(x + 0.18, ridge.coef_, width=0.36, label="Ridge")
plt.title("OLSとRidgeの標準化係数")
plt.xlabel("説明変数")
plt.ylabel("標準化係数")
plt.xticks(x, corr_features, rotation=25, ha="right")
plt.grid(True, axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	説明変数	OLS標準化係数	Ridge標準化係数
0	temperature	0.981	1.252
1	sensor_temp_2	1.937	1.523
2	cooling_flow	0.024	-0.079
3	speed	4.368	4.201
4	vibration	1.663	1.597

OLS評価RMSE   : 2.880
Ridge評価RMSE : 2.878

png

結果の読み取り

温度、別温度センサー、冷却流量は相関しており、OLSでは効果を変数間で不安定に分配しやすくなります。Ridgeは標準化係数を縮小し、類似情報を持つ変数の極端な係数を抑えます。RMSEが大きく悪化せず係数が安定するなら、監視モデルとして有力です。ただし、縮小係数をそのまま物理法則や因果効果と解釈しません。 $\alpha$ は本番では交差検証で決めます。

No.063：Lasso回帰 — 改善候補となる変数を絞る

実務での意味

多数のセンサーを収集しても、全てを常時保守するコストは小さくありません。Lassoは不要な係数を0にできるため、予測に必要な候補の絞り込みや簡潔なモデル作成に使えます。

分析・モデル化の考え方

LassoはL1罰則を使い、

\min_{\beta}\left\{\sum_i(y_i-\beta_0-x_i^T\beta)^2+\alpha\sum_j|\beta_j|\right\}

を解きます。変数の尺度で罰則の効き方が変わるので標準化が必須です。相関の強い変数群では、そのうち一つだけを恣意的に残すことがあるため、選択結果を物理的な唯一原因とみなしません。

Pythonで確認する

wide_features = corr_features + ["moisture", "maintenance_days", "ambient_humidity"]
Xl_train, Xl_test, yl_train, yl_test = train_test_split(
    df[wide_features], df["energy"], test_size=0.25, random_state=SEED
)
lasso_pipe = make_pipeline(StandardScaler(), Lasso(alpha=0.45, max_iter=20000))
lasso_pipe.fit(Xl_train, yl_train)
lasso_coef = lasso_pipe.named_steps["lasso"].coef_
lasso_table = pd.DataFrame({
    "説明変数": wide_features, "標準化係数": lasso_coef,
    "選択": np.where(np.abs(lasso_coef) > 1e-8, "採用候補", "0（除外候補）"),
})
display(lasso_table.round(3))
print(f"評価RMSE: {mean_squared_error(yl_test, lasso_pipe.predict(Xl_test))**0.5:.3f}")

plt.barh(wide_features, lasso_coef)
plt.title("Lassoによる変数選択")
plt.xlabel("標準化係数")
plt.ylabel("説明変数")
plt.grid(True, axis="x", alpha=0.3)
plt.tight_layout()
plt.show()

	説明変数	標準化係数	選択
0	temperature	0.779	採用候補
1	sensor_temp_2	1.639	採用候補
2	cooling_flow	-0.000	0（除外候補）
3	speed	3.886	採用候補
4	vibration	1.160	採用候補
5	moisture	0.000	0（除外候補）
6	maintenance_days	0.000	0（除外候補）
7	ambient_humidity	0.000	0（除外候補）

評価RMSE: 3.013

png

結果の読み取り

係数が0になった変数は、今回の罰則強度と他変数の存在下で追加的な予測情報が小さい候補です。「工程に影響しない」と断定する結果ではありません。選択は標本や $\alpha$ に依存するため、再標本化での選択頻度、保全費用、測定の信頼性を併せて、センサー廃止のような不可逆な判断を避けます。

No.064：ElasticNet — 選択と相関変数の保持を両立する

実務での意味

製造データには、温度センサー群や振動の周波数帯のような相関変数のまとまりがあります。Lassoの簡潔さを保ちつつ、関連する変数群をまとめて扱いたいときにElasticNetが候補になります。

分析・モデル化の考え方

ElasticNetはL1とL2を組み合わせ、

\min_{\beta}\left\{\frac{1}{2n}\|y-X\beta\|_2^2+\alpha\left(\rho\|\beta\|_1+\frac{1-\rho}{2}\|\beta\|_2^2\right)\right\}

を解きます。 $\rho=1$ はLasso、 $\rho=0$ に近づくほどRidge的です。ここでは同じ学習・評価分割でLassoと比較します。

Pythonで確認する

enet_pipe = make_pipeline(
    StandardScaler(), ElasticNet(alpha=0.35, l1_ratio=0.45, max_iter=20000)
)
enet_pipe.fit(Xl_train, yl_train)
enet_coef = enet_pipe.named_steps["elasticnet"].coef_
regularized_comparison = pd.DataFrame({
    "説明変数": wide_features, "Lasso": lasso_coef, "ElasticNet": enet_coef,
})
display(regularized_comparison.round(3))

model_rmse = pd.DataFrame({
    "モデル": ["Lasso", "ElasticNet"],
    "評価RMSE": [
        mean_squared_error(yl_test, lasso_pipe.predict(Xl_test))**0.5,
        mean_squared_error(yl_test, enet_pipe.predict(Xl_test))**0.5,
    ],
})
display(model_rmse.round(3))

x = np.arange(len(wide_features))
plt.plot(x, lasso_coef, "o-", label="Lasso")
plt.plot(x, enet_coef, "s-", label="ElasticNet")
plt.title("正則化モデルの標準化係数比較")
plt.xlabel("説明変数")
plt.ylabel("標準化係数")
plt.xticks(x, wide_features, rotation=25, ha="right")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	説明変数	Lasso	ElasticNet
0	temperature	0.779	1.117
1	sensor_temp_2	1.639	1.171
2	cooling_flow	-0.000	-0.234
3	speed	3.886	3.501
4	vibration	1.160	1.214
5	moisture	0.000	0.049
6	maintenance_days	0.000	0.050
7	ambient_humidity	0.000	0.009

	モデル	評価RMSE
0	Lasso	3.013
1	ElasticNet	3.064

png

結果の読み取り

ElasticNetはL2成分により、相関する温度系変数の情報を複数係数へ分けて保持しやすくなります。Lassoより変数数は増え得ますが、相関群を一つだけ残す不安定さを緩和できます。実務ではRMSE、係数の再現性、計測コスト、説明のしやすさを評価軸にし、 $\alpha$ と $\rho$ を交差検証で選びます。

No.065：ロジスティック回帰 — バッチ不適合確率を予測する

実務での意味

目的変数が合格・不合格の二値なら、必要なのは連続値そのものではなく不適合確率です。確率が閾値を超えたバッチを追加検査へ回すなど、見逃し費用と検査費用を結び付けられます。

分析・モデル化の考え方

ロジスティック回帰は

\log\frac{p_i}{1-p_i}=\beta_0+x_i^T\beta, \qquad p_i=P(Y_i=1\mid x_i)

と置きます。 $e^{\beta_j}$ は他の変数を一定とした1単位増加あたりのオッズ比です。ROC-AUCは順位付け、log lossは確率の確からしさを評価します。運用閾値は0.5固定ではなく、見逃しと過剰検査の費用から決めます。

Pythonで確認する

clf_features = ["temperature", "vibration", "moisture", "maintenance_days"]
Xc_train, Xc_test, yc_train, yc_test = train_test_split(
    df[clf_features], df["defect"], test_size=0.25, stratify=df["defect"], random_state=SEED
)
logit_pipe = make_pipeline(StandardScaler(), LogisticRegression(C=10, max_iter=5000))
logit_pipe.fit(Xc_train, yc_train)
prob = logit_pipe.predict_proba(Xc_test)[:, 1]
logit_coef = logit_pipe.named_steps["logisticregression"].coef_[0]
display(pd.DataFrame({
    "説明変数": clf_features, "標準化係数": logit_coef, "標準偏差1個分のオッズ比": np.exp(logit_coef),
}).round(3))
print(f"不適合率 : {yc_test.mean():.3f}")
print(f"ROC-AUC  : {roc_auc_score(yc_test, prob):.3f}")
print(f"Log loss : {log_loss(yc_test, prob):.3f}")

bins = pd.qcut(prob, q=5, duplicates="drop")
calibration = pd.DataFrame({"予測確率": prob, "実績": yc_test.to_numpy(), "bin": bins}).groupby(
    "bin", observed=True
).agg(平均予測確率=("予測確率", "mean"), 実績不適合率=("実績", "mean"))
display(calibration.round(3))
plt.plot(calibration["平均予測確率"], calibration["実績不適合率"], "o-", label="5群の実績")
plt.plot([0, 1], [0, 1], "--", color="tab:red", label="理想")
plt.title("不適合確率の簡易キャリブレーション確認")
plt.xlabel("平均予測不適合確率")
plt.ylabel("実績不適合率")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	説明変数	標準化係数	標準偏差1個分のオッズ比
0	temperature	0.501	1.651
1	vibration	0.230	1.259
2	moisture	0.478	1.613
3	maintenance_days	0.940	2.559

不適合率 : 0.143
ROC-AUC  : 0.749
Log loss : 0.373

	平均予測確率	実績不適合率
bin
(0.00366, 0.0371]	0.025	0.048
(0.0371, 0.0778]	0.054	0.048
(0.0778, 0.153]	0.106	0.048
(0.153, 0.239]	0.200	0.190
(0.239, 0.646]	0.382	0.381

png

結果の読み取り

オッズ比が1より大きい変数は、標準偏差1個分の増加と不適合オッズ上昇が関連します。ROC-AUCが高くても、予測確率と実績率が一致するとは限らないため、5群の較正も確認します。追加検査へ回す閾値は、見逃し損失、検査能力、対象比率を使って決め、品種別・期間別にも性能を監視します。

No.066：ポアソン回帰 — 稼働時間を考慮して微停止件数を比較する

実務での意味

8時間稼働したバッチと6時間のバッチでは、同じ2件でも発生率が違います。ポアソン回帰に稼働時間のoffsetを入れると、観測機会をそろえた停止率として工程条件を比較できます。

分析・モデル化の考え方

$Y_i\sim\mathrm{Poisson}(\mu_i)$ とし、

\log\mu_i=\log t_i+\beta_0+x_i^T\beta

と置きます。 $t_i$ は稼働時間で、係数1に固定した $\log t_i$ がoffsetです。 $e^{\beta_j}$ は他条件一定での発生率比です。ポアソン分布は条件付き平均と分散が等しいと仮定します。

Pythonで確認する

count_features = ["vibration", "maintenance_days"]
X_count = sm.add_constant(df[count_features])
poisson_model = sm.GLM(
    df["micro_stops"], X_count, family=sm.families.Poisson(),
    offset=np.log(df["operating_hours"]),
).fit()
pearson_ratio = np.sum(poisson_model.resid_pearson**2) / poisson_model.df_resid
poisson_table = pd.DataFrame({
    "係数": poisson_model.params,
    "発生率比 exp(係数)": np.exp(poisson_model.params),
    "p値": poisson_model.pvalues,
})
display(poisson_table.round(4))
print(f"Pearsonカイ二乗 / 自由度: {pearson_ratio:.3f}")

pred_count = poisson_model.predict(X_count, offset=np.log(df["operating_hours"]))
plt.scatter(pred_count, df["micro_stops"], alpha=0.45)
plt.title("ポアソン回帰：予測微停止件数と実績")
plt.xlabel("予測微停止件数 [件/バッチ]")
plt.ylabel("実績微停止件数 [件/バッチ]")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

	係数	発生率比 exp(係数)
const	-1.9616	0.1406
vibration	0.4876	1.6284
maintenance_days	0.0080	1.0080

Pearsonカイ二乗 / 自由度: 2.832

png

結果の読み取り

exp(係数) は1単位増加あたりの微停止発生率比です。offsetにより、稼働時間の長さを説明変数として推定せず、曝露量として補正しています。一方、Pearsonカイ二乗を自由度で割った値が1を大きく上回るなら過分散の兆候で、ポアソンの標準誤差が過小になる可能性があります。予測点の散らばりと併せ、次の負の二項回帰を検討します。

No.067：負の二項回帰 — 微停止件数の過分散を織り込む

実務での意味

微停止は、材料ロットや設備状態など未観測要因が重なる日に集中します。平均件数だけでは表せない日間差を無視すると、「異常」と判定する範囲が狭くなり、過剰なアラートを生むことがあります。

分析・モデル化の考え方

負の二項回帰も対数リンクを使いますが、分散を

\mathrm{Var}(Y_i\mid X_i)=\mu_i+\alpha\mu_i^2

とし、 $\alpha>0$ で過分散を許します。ポアソンは $\alpha\to0$ の極限です。AICは同じデータに当てたモデル間の適合度と複雑さの比較指標で、小さい方を候補にします。

Pythonで確認する

nb_model = sm.NegativeBinomial(
    df["micro_stops"], X_count, offset=np.log(df["operating_hours"])
).fit(disp=False)
nb_params = nb_model.params
nb_table = pd.DataFrame({
    "推定値": nb_params,
    "exp(推定値)": np.exp(nb_params),
    "p値": nb_model.pvalues,
})
display(nb_table.round(4))
comparison_count = pd.DataFrame({
    "モデル": ["Poisson GLM", "Negative Binomial"],
    "AIC": [poisson_model.aic, nb_model.aic],
    "対数尤度": [poisson_model.llf, nb_model.llf],
})
display(comparison_count.round(2))

plt.bar(comparison_count["モデル"], comparison_count["AIC"], color=["tab:blue", "tab:orange"])
plt.title("件数モデルのAIC比較")
plt.xlabel("モデル")
plt.ylabel("AIC（小さいほど候補）")
plt.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	推定値	exp(推定値)
const	-1.9694	0.1395
vibration	0.4821	1.6195
maintenance_days	0.0083	1.0083
alpha	0.5726	1.7729

	モデル	AIC	対数尤度
0	Poisson GLM	2141.00	-1067.50
1	Negative Binomial	1858.13	-925.06

png

結果の読み取り

alpha が正で、負の二項回帰のAICがポアソンより小さければ、未観測の異質性を含む過分散を表現した価値があります。係数の exp は発生率比ですが、alpha は説明変数の効果ではありません。ゼロが極端に多い、時間依存や設備内相関がある場合は、ゼロ過剰モデル、混合モデル、時系列モデルも候補です。

No.068：GLM — 歪んだ復旧工数をGamma回帰で表す

実務での意味

復旧工数は0より大きく、少数の長時間作業によって右に歪みます。正規誤差の線形回帰では負の予測が出る可能性があるため、正の連続量に合うGamma分布と対数リンクを使います。

分析・モデル化の考え方

一般化線形モデル（GLM）は、確率分布とリンク関数を組み合わせます。Gamma回帰では

Y_i\sim\mathrm{Gamma},\qquad \log E(Y_i\mid X_i)=\beta_0+x_i^T\beta

とし、予測平均を必ず正に保ちます。係数を指数変換した $e^{\beta_j}$ は平均工数の倍率で、 $(e^{\beta_j}-1)\times100\%$ が概ね相対変化率です。

Pythonで確認する

gamma_features = ["vibration", "maintenance_days"]
X_gamma = sm.add_constant(df[gamma_features])
gamma_model = sm.GLM(
    df["recovery_hours"], X_gamma,
    family=sm.families.Gamma(link=sm.families.links.Log()),
).fit()
gamma_table = pd.DataFrame({
    "係数": gamma_model.params,
    "平均工数の倍率": np.exp(gamma_model.params),
    "p値": gamma_model.pvalues,
})
display(gamma_table.round(4))

gamma_pred = gamma_model.predict(X_gamma)
plt.scatter(gamma_pred, df["recovery_hours"], alpha=0.45)
plt.title("Gamma GLM：予測復旧工数と実績")
plt.xlabel("予測平均復旧工数 [時間]")
plt.ylabel("実績復旧工数 [時間]")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

	係数	平均工数の倍率
const	-1.2942	0.2741
vibration	0.2387	1.2697
maintenance_days	0.0123	1.0124

png

結果の読み取り

Gamma GLMの予測値は正に保たれ、係数は加算時間ではなく倍率として解釈します。例えば保全後経過日数の倍率を30乗すれば、30日増加時の平均工数倍率を計算できます。分布を変えても外れ値が無害になるわけではないため、残差、影響点、品種別の適合を確認します。GLMは分布に合わせた平均モデルであり、個々の復旧時間を完全に当てるものではありません。

No.069：GAM — 工程条件の非線形な影響を捉える

実務での意味

品質損失は「温度が高いほど悪い」とは限らず、低すぎても高すぎても悪化することがあります。直線だけでは安全な操業帯を表せないため、滑らかな曲線で各変数の影響を学習します。

分析・モデル化の考え方

一般化加法モデル（GAM）は

g\{E(Y\mid X)\}=\beta_0+f_1(x_1)+f_2(x_2)+\cdots+f_p(x_p)

とし、 $f_j$ をスプライン基底の重み付き和で表します。ここでは各列をBスプラインへ展開し、Ridgeで曲線の過度な揺れを抑える加法モデルを作ります。交互作用を入れないため、説明しやすい一方、温度×速度の複合効果は表現しません。

Pythonで確認する

gam_features = ["temperature", "speed", "vibration"]
Xg_train, Xg_test, yg_train, yg_test = train_test_split(
    df[gam_features], df["quality_loss"], test_size=0.25, random_state=SEED
)
linear_quality = LinearRegression().fit(Xg_train, yg_train)
gam = make_pipeline(
    SplineTransformer(n_knots=6, degree=3, include_bias=False), Ridge(alpha=2.0)
).fit(Xg_train, yg_train)
print(f"線形モデル 評価RMSE: {mean_squared_error(yg_test, linear_quality.predict(Xg_test))**0.5:.3f}")
print(f"GAM相当   評価RMSE: {mean_squared_error(yg_test, gam.predict(Xg_test))**0.5:.3f}")

temp_grid = np.linspace(df.temperature.quantile(0.01), df.temperature.quantile(0.99), 160)
profile = pd.DataFrame({
    "temperature": temp_grid,
    "speed": df.speed.median(),
    "vibration": df.vibration.median(),
})
plt.scatter(df.temperature, df.quality_loss, alpha=0.18, label="実績")
plt.plot(temp_grid, gam.predict(profile), color="tab:red", lw=2.5, label="GAMの部分プロファイル")
plt.title("GAMで捉えた温度と品質損失の非線形関係")
plt.xlabel("工程温度 [℃]")
plt.ylabel("品質損失指数")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

線形モデル 評価RMSE: 4.754
GAM相当   評価RMSE: 1.769

png

結果の読み取り

直線モデルよりGAM相当モデルの評価RMSEが小さく、部分プロファイルがU字なら、温度に中間的な良好領域があることを示唆します。ただし曲線端はデータが少なく不安定です。図は速度と振動を中央値に固定した条件付き予測で、単純な散布図平均ではありません。設定標準へ反映する前に、許容範囲内での確認実験、交互作用、信頼区間を検討します。

No.070：Quantile Regression — 長引くバッチの上側リスクを見積もる

実務での意味

平均サイクル時間だけで要員や後工程の能力を決めると、遅延が集中する条件を見落とします。分位点回帰は中央値や90%分位点を直接モデル化し、「10回に1回程度はこれを超える」水準の計画に使えます。

分析・モデル化の考え方

$\tau$ 分位点回帰は、非対称なチェック損失

\rho_\tau(u)=u\{\tau-\mathbb{1}(u<0)\}

の合計を最小化します。 $\tau=0.5$ は条件付き中央値、 $\tau=0.9$ は条件付き90%分位点です。誤差分散が条件によって変わる場合、平均回帰では見えない上側リスクの変化を係数として表せます。

Pythonで確認する

quant_features = ["speed", "vibration", "maintenance_days"]
X_quant = sm.add_constant(df[quant_features])
q50 = sm.QuantReg(df["cycle_time"], X_quant).fit(q=0.5)
q90 = sm.QuantReg(df["cycle_time"], X_quant).fit(q=0.9)
quant_table = pd.DataFrame({"中央値 q=0.5": q50.params, "上側 q=0.9": q90.params})
display(quant_table.round(3))

vib_grid = np.linspace(df.vibration.quantile(0.01), df.vibration.quantile(0.99), 160)
quant_profile = pd.DataFrame({
    "const": 1.0, "speed": df.speed.median(), "vibration": vib_grid,
    "maintenance_days": df.maintenance_days.median(),
})
plt.scatter(df.vibration, df.cycle_time, alpha=0.2, label="実績")
plt.plot(vib_grid, q50.predict(quant_profile), lw=2.2, label="条件付き中央値")
plt.plot(vib_grid, q90.predict(quant_profile), lw=2.2, label="条件付き90%分位点")
plt.title("振動水準別のサイクル時間分位点")
plt.xlabel("設備振動 [mm/s]")
plt.ylabel("サイクル時間 [分] ")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

	中央値 q=0.5	上側 q=0.9
const	41.485	47.253
speed	0.149	0.113
vibration	4.090	6.336
maintenance_days	0.015	0.019

png

結果の読み取り

90%分位点の線は、各条件でサイクル時間の約90%が下側に入る水準です。振動係数が中央値より90%分位点で大きければ、振動は典型値だけでなく長時間側の尾を広げることを示唆します。納期や要員計画では上側分位点、標準作業の代表値には中央値というように、損失に応じて分位点を選びます。これは将来値の同時予測区間とは異なる点に注意します。

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

目的変数の尺度からモデルを選ぶ：連続量、二値、件数、正の歪んだ量では、妥当な分布とリンクが異なります。
相関変数には安定性の観点を加える：OLSの係数だけで施策を決めず、Ridge、Lasso、ElasticNetで予測・簡潔さ・再現性を比較します。
曝露量と過分散を無視しない：停止件数は稼働時間で補正し、平均を超えるばらつきがあれば負の二項回帰を検討します。
平均だけをKPIにしない：GAMで非線形な安全領域を、分位点回帰で上側リスクを捉えると、条件標準と能力計画を分けて設計できます。
予測と因果を区別する：高精度なモデルでも、条件を操作した効果は保証しません。改善実施には現場知識と確認実験が必要です。

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

1. 予測時点と意思決定を定義する

いつ、誰が、何を変更するための予測かを先に決めます。予測時点より後に判明する情報を説明変数へ入れると、検証精度だけが高いデータリークになります。

2. データ生成過程を管理する

単位、欠損、センサー校正、設備・品種・作業者ID、稼働時間、保全履歴を揃えます。同一設備の反復観測や時系列分割を無視したランダム分割は、性能を楽観的に見せることがあります。

3. モデル選択と閾値を業務損失で評価する

RMSE、AUC、AICだけでなく、不適合の見逃し費用、追加検査費、停止アラート対応費、要員不足損失へ翻訳します。学習・検証・最終評価を分離し、ハイパーパラメータは検証データまたは交差検証で決めます。

4. 運用後の監視と責任分界を決める

入力分布、欠損率、較正、残差、設備別性能を継続監視します。モデルが判断を支援する範囲、現場が覆せる条件、再学習・停止基準、変更履歴を明文化します。

5. 改善効果は実験で確認する

回帰係数から改善仮説を作り、安全・品質上許される範囲でDOEや段階導入を行います。予測モデルの関連と、設定変更による因果効果を分けて検証します。

まとめ

No.061〜No.070では、同じ製造データでも目的変数と意思決定に応じてモデルを変える必要性を確認しました。

線形回帰を基準に、Ridge・Lasso・ElasticNetで相関と変数数を制御した
ロジスティック回帰で不適合確率、ポアソン／負の二項回帰で稼働時間あたり微停止をモデル化した
Gamma GLMで正の歪んだ工数、GAMで非線形関係、分位点回帰で上側リスクを捉えた

実務で重要なのは、最も複雑なモデルを選ぶことではありません。データの発生機構と利用目的を明確にし、未知データで検証し、予測を費用・品質・納期の意思決定へ接続することです。

法人向けのご相談

数理工房では、製造業向けに以下のご支援が可能です。

品質・設備・生産実績を用いた予測モデルの設計と検証
不適合、停止件数、復旧工数に応じたGLM・機械学習の選定
PoCから現場運用までのKPI、閾値、監視方法の設計
回帰分析、統計モデリング、Python活用に関する企業研修

「手元のデータで何を予測できるか整理したい」「精度は出たが現場の判断につながらない」といった段階からでも、業務課題とデータ生成過程を確認しながら進めます。

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

製造業の回帰分析をPythonで実践｜線形回帰・GLM・GAM 10本ノック

工程条件から品質・停止・工数を予測する：製造業の回帰モデル10本ノック

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.061：線形回帰 — 工程条件から電力原単位を説明する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.062：Ridge回帰 — 相関するセンサーの係数を安定化する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.063：Lasso回帰 — 改善候補となる変数を絞る

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.064：ElasticNet — 選択と相関変数の保持を両立する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.065：ロジスティック回帰 — バッチ不適合確率を予測する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.066：ポアソン回帰 — 稼働時間を考慮して微停止件数を比較する

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.067：負の二項回帰 — 微停止件数の過分散を織り込む

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.068：GLM — 歪んだ復旧工数をGamma回帰で表す

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.069：GAM — 工程条件の非線形な影響を捉える

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.070：Quantile Regression — 長引くバッチの上側リスクを見積もる

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

1. 予測時点と意思決定を定義する

2. データ生成過程を管理する

3. モデル選択と閾値を業務損失で評価する

4. 運用後の監視と責任分界を決める

5. 改善効果は実験で確認する

まとめ

法人向けのご相談

Related Posts

不確実性を設備・在庫・生産能力の判断へつなぐ：製造業シミュレーション10本ノック

需要変動を読み解き、生産計画へつなぐ：製造業の時系列分析10本ノック

不確かさを更新して改善案を選ぶ：製造業のベイズ統計10本ノック

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

観測データから工程の真の姿を推定する：製造業の統計的推定10本ノック