100本ノック / 確率統計 / 確率・統計マーケティング応用100本ノック
製造業の機械学習をPythonで実践｜予測・因果推論・最適化10本ノック
予測から現場アクションへ：製造業の機械学習10本ノック
本記事では、架空の産業機器メーカーを題材に、回帰、分類、クラスタリング、次元削減、異常検知、時系列予測、因果推論、ベイズ最適化、強化学習、グラフ機械学習を、製造業の意思決定へつなげます。
目的はモデル精度を競うことではありません。工数見積、品質判定、顧客対応、設備監視、需要計画、改善施策、加工条件、保全方針、サプライヤーリスクについて、どのデータを使い、どの指標で評価し、誰の判断をどう支援するかを具体化します。データはすべて Python で生成し、外部データには依存しません。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。 掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
対象企業は、受注仕様に応じてポンプ部品を加工・組立し、国内外へ販売しています。現場には、加工時間の見積誤差、不良の見逃し、設備停止、需要変動、改善効果の過大評価、サプライヤー障害の波及といった課題があります。
機械学習は、これらを自動的に解決する魔法ではありません。予測値を「残業を承認する」「検査を追加する」「保全を実施する」といった意思決定へ変換し、実績との差を監視して初めて価値を生みます。本記事では、各モデルの出力と、その出力を使う業務を対にして説明します。
現場でよくある状況
- 標準工数が製品群別の平均値しかなく、高難度品で納期回答を外す
- 不良予測モデルを作ったが、見逃しと過剰検査の費用を比較していない
- センサー項目が多く、どの変化を設備担当者へ示すべきか分からない
- 改善前後の平均を比べただけで、案件難度や設備差を改善効果と取り違える
- 予測・最適化結果が既存の作業標準、承認権限、例外処理につながっていない
なぜこの問題は判断が難しいのか
製造データは時系列・設備・品種・作業者・取引先が相互に依存します。ランダム分割では未来情報が学習側へ混ざり、平均精度だけでは重大不良の見逃しが隠れます。また、相関が高くても、条件変更による効果とは限りません。
したがって、目的変数、予測時点で利用可能な説明変数、評価単位、誤判定費用、適用範囲を先に定義します。モデルの選択は、その後の手段です。
今回扱うノックの全体像
| No. | テーマ | 製造業での主な問い | 主な評価・出力 |
|---|---|---|---|
| 061 | 回帰 | 加工時間を何分と見積もるか | MAE、残差 |
| 062 | 分類 | 不良リスク品を追加検査するか | 再現率、混同行列 |
| 063 | クラスタリング | 顧客をどの支援方針に分けるか | セグメント特性 |
| 064 | 次元削減 | 多数のセンサー変化をどう要約するか | 寄与率、主成分得点 |
| 065 | 異常検知 | 正解ラベルが少ない設備異常をどう拾うか | 異常スコア、上位候補 |
| 066 | 時系列予測 | 来月の受注量と必要能力はいくつか | 時系列MAE、予測推移 |
| 067 | 因果推論 | 新治具は本当に加工時間を短縮したか | IPWによるATE |
| 068 | ベイズ最適化 | 少ない試作で加工条件をどう探すか | 獲得関数、最良条件 |
| 069 | 強化学習 | 劣化状態ごとに保全するか継続するか | 方策、累積報酬 |
| 070 | グラフ機械学習 | 取引網から供給リスクをどう捉えるか | 近傍特徴、リスク確率 |
Python 環境の準備
NumPy と pandas はデータ生成・集計、scikit-learn は前処理と機械学習、SciPy は獲得関数、NetworkX は取引ネットワーク、matplotlib は可視化に利用します。default_rng(42) を使い、再実行しても同じ架空データになるようにします。
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib
import networkx as nx
import sklearn
from IPython.display import display
from scipy.stats import norm
from sklearn.cluster import KMeans
from sklearn.compose import ColumnTransformer
from sklearn.decomposition import PCA
from sklearn.ensemble import IsolationForest
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, Matern, WhiteKernel
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.metrics import (ConfusionMatrixDisplay, RocCurveDisplay, accuracy_score,
classification_report, mean_absolute_error, r2_score,
roc_auc_score)
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler
rng = np.random.default_rng(42)
print("numpy :", np.__version__)
print("pandas :", pd.__version__)
print("matplotlib :", matplotlib.__version__)
print("sklearn :", sklearn.__version__)
print("networkx :", nx.__version__)numpy : 2.5.1
pandas : 3.0.3
matplotlib : 3.11.0
sklearn : 1.9.0
networkx : 3.6.1架空データの作成
直近420件の加工実績を想定します。受注時点で分かる材質、製品群、ロット数、複雑度、設備、段取り回数に加え、加工後に分かる振動・温度・寸法偏差・不良フラグを持たせます。event_date は時系列評価に使えるよう日付順にします。
重要なのは「予測時点」です。たとえば受注時の工数見積には、加工後に測る振動や寸法偏差を使えません。利用可能性を無視すると、検証精度だけが高いデータリークになります。
n = 420
event_date = pd.date_range("2025-01-01", periods=n, freq="D")
material = rng.choice(["Steel", "Aluminum", "Resin"], n, p=[0.50, 0.30, 0.20])
product = rng.choice(["Pump-A", "Pump-B", "Pump-C"], n, p=[0.45, 0.35, 0.20])
machine = rng.choice(["MC-1", "MC-2", "MC-3"], n, p=[0.40, 0.35, 0.25])
lot_size = rng.integers(20, 151, n)
complexity = np.clip(rng.normal(5.0, 1.8, n), 1, 10)
setup_count = rng.integers(1, 5, n)
mat_effect = pd.Series(material).map({"Steel": 18, "Aluminum": 8, "Resin": 3}).to_numpy()
prod_effect = pd.Series(product).map({"Pump-A": 0, "Pump-B": 12, "Pump-C": 25}).to_numpy()
machine_effect = pd.Series(machine).map({"MC-1": 2, "MC-2": -4, "MC-3": 6}).to_numpy()
cycle_time = (35 + 0.48 * lot_size + 7.2 * complexity + 5.5 * setup_count
+ mat_effect + prod_effect + machine_effect + rng.normal(0, 8, n))
vibration = np.clip(1.6 + 0.12 * complexity + 0.003 * lot_size
+ (machine == "MC-3") * 0.35 + rng.normal(0, 0.22, n), 0.5, None)
temperature = 53 + 1.7 * complexity + (machine == "MC-1") * 3 + rng.normal(0, 3.2, n)
dim_error = rng.normal(0, 0.018 + 0.003 * complexity, n) + (vibration - 2.4) * 0.012
logit = -5.0 + 0.75 * complexity + 0.85 * (vibration - 2.2) + 18 * np.abs(dim_error)
defect_prob = 1 / (1 + np.exp(-logit))
defect = rng.binomial(1, defect_prob)
production = pd.DataFrame({
"event_date": event_date, "material": material, "product": product,
"machine": machine, "lot_size": lot_size, "complexity": complexity.round(2),
"setup_count": setup_count, "cycle_time_min": cycle_time.round(1),
"vibration_mm_s": vibration.round(3), "temperature_c": temperature.round(2),
"dimension_error_mm": dim_error.round(4), "defect": defect,
})
display(production.head())
print(f"行数: {len(production):,} / 不良率: {production['defect'].mean():.1%}")| event_date | material | product | machine | lot_size | complexity | setup_count | cycle_time_min | vibration_mm_s | temperature_c | dimension_error_mm | defect | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2025-01-01 | Aluminum | Pump-B | MC-1 | 77 | 4.51 | 1 | 127.1 | 2.540 | 60.91 | 0.0375 | 0 |
| 1 | 2025-01-02 | Steel | Pump-A | MC-3 | 25 | 3.33 | 1 | 91.7 | 2.363 | 61.24 | 0.0212 | 0 |
| 2 | 2025-01-03 | Resin | Pump-A | MC-2 | 47 | 3.44 | 2 | 80.7 | 2.553 | 57.05 | -0.0436 | 0 |
| 3 | 2025-01-04 | Aluminum | Pump-A | MC-2 | 72 | 2.06 | 3 | 118.6 | 1.869 | 51.99 | -0.0139 | 0 |
| 4 | 2025-01-05 | Steel | Pump-B | MC-1 | 28 | 5.43 | 4 | 143.7 | 2.383 | 68.49 | 0.0055 | 0 |
行数: 420 / 不良率: 36.9%No.061：回帰 — 受注仕様から加工時間を見積もる
実務での意味
加工時間の予測は、見積原価、納期回答、設備負荷、残業判断の共通入力です。平均標準時間だけでなく、材質・品種・設備・難度を反映すると、高難度案件の過小見積を減らせます。
分析・モデル化の考え方
連続値 を予測する線形回帰を使います。
カテゴリ変数は One-hot 化し、時系列の前半80%を学習、後半20%を評価します。誤差は現場で解釈しやすい平均絶対誤差 と、説明率を表す で確認します。
Pythonで確認する
features_061 = ["material", "product", "machine", "lot_size", "complexity", "setup_count"]
cat_061 = ["material", "product", "machine"]
num_061 = ["lot_size", "complexity", "setup_count"]
split = int(len(production) * 0.8)
X_train, X_test = production.loc[:split-1, features_061], production.loc[split:, features_061]
y_train, y_test = production.loc[:split-1, "cycle_time_min"], production.loc[split:, "cycle_time_min"]
prep_061 = ColumnTransformer([
("category", OneHotEncoder(handle_unknown="ignore"), cat_061),
("numeric", StandardScaler(), num_061),
])
model_061 = make_pipeline(prep_061, LinearRegression()).fit(X_train, y_train)
pred_061 = model_061.predict(X_test)
metrics_061 = pd.DataFrame({
"指標": ["MAE（分）", "R²"],
"値": [mean_absolute_error(y_test, pred_061), r2_score(y_test, pred_061)],
})
display(metrics_061.round(3))
fig, ax = plt.subplots(figsize=(7, 4))
ax.scatter(pred_061, y_test - pred_061, alpha=0.65)
ax.axhline(0, color="red", linestyle="--")
ax.set_title("加工時間予測の残差")
ax.set_xlabel("予測加工時間（分）")
ax.set_ylabel("実績 - 予測（分）")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| 指標 | 値 | |
|---|---|---|
| 0 | MAE（分） | 6.898 |
| 1 | R² | 0.902 |
結果の読み取り
MAE は、見積値が平均して何分ずれるかを示します。残差がゼロ周辺に散らばれば大きな系統誤差はありませんが、長時間側で残差が広がるなら、難度区分の細分化や予測区間が必要です。実務では MAE をそのまま納期バッファにせず、製品群別の過小見積率と、過小見積が残業・納期へ与える費用を確認します。
No.062：分類 — 不良リスクに応じて追加検査を割り当てる
実務での意味
全数検査は品質を守りやすい一方、検査能力とリードタイムを消費します。分類モデルは、不良確率が高いロットを重点検査へ回し、見逃しと検査負荷を両立する判断材料になります。
分析・モデル化の考え方
ロジスティック回帰で不良確率を推定します。
不良は少数クラスなので class_weight="balanced" を使います。ROC-AUCだけでなく混同行列を確認し、特に不良を正常と判定する偽陰性を管理します。今回は加工後の品質ゲートを想定し、振動・温度・寸法偏差を入力に含めます。
Pythonで確認する
features_062 = ["complexity", "vibration_mm_s", "temperature_c", "dimension_error_mm"]
X_062 = production[features_062]
y_062 = production["defect"]
X_train, X_test, y_train, y_test = train_test_split(
X_062, y_062, test_size=0.25, random_state=42, stratify=y_062
)
model_062 = make_pipeline(
StandardScaler(), LogisticRegression(class_weight="balanced", random_state=42)
).fit(X_train, y_train)
prob_062 = model_062.predict_proba(X_test)[:, 1]
pred_062 = (prob_062 >= 0.45).astype(int)
print(f"ROC-AUC: {roc_auc_score(y_test, prob_062):.3f}")
print(classification_report(y_test, pred_062, target_names=["良品", "不良"], digits=3))
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
ConfusionMatrixDisplay.from_predictions(y_test, pred_062, display_labels=["良品", "不良"], ax=axes[0], colorbar=False)
axes[0].set_title("追加検査ルールの混同行列")
axes[0].set_xlabel("予測クラス")
axes[0].set_ylabel("実績クラス")
axes[0].grid(False)
RocCurveDisplay.from_predictions(y_test, prob_062, ax=axes[1])
axes[1].set_title("不良分類のROC曲線")
axes[1].set_xlabel("偽陽性率")
axes[1].set_ylabel("真陽性率")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()ROC-AUC: 0.811
precision recall f1-score support
良品 0.837 0.621 0.713 66
不良 0.554 0.795 0.653 39
accuracy 0.686 105
macro avg 0.695 0.708 0.683 105
weighted avg 0.732 0.686 0.691 105
結果の読み取り
混同行列の右下は検出できた不良、左下は見逃した不良です（横軸が予測、縦軸が実績）。しきい値を下げれば不良再現率は上がりやすい一方、良品の追加検査も増えます。実装時は、見逃し費用 × 偽陰性件数 + 追加検査費用 × 検査件数 をしきい値ごとに比較し、品種別にも性能を監視します。
No.063：クラスタリング — 法人顧客を支援ニーズで分ける
実務での意味
産業財では、売上規模だけで顧客対応を分けると、少額でも高頻度の技術問い合わせがある顧客や、納期変動の大きい顧客を見落とします。クラスタリングは、似た購買・サービス特性を持つ顧客群を探索する手段です。
分析・モデル化の考え方
K-means は、各点と所属クラスタ中心の二乗距離の総和を最小化します。
単位の大きい年間売上だけが距離を支配しないよう標準化します。クラスタ番号自体に優劣はなく、集計表を見て業務上の名前を後付けします。
Pythonで確認する
n_customer = 75
base_segment = rng.choice(3, n_customer, p=[0.40, 0.35, 0.25])
annual_sales = np.exp(rng.normal(np.choose(base_segment, [4.0, 5.0, 4.5]), 0.35))
order_frequency = np.clip(rng.normal(np.choose(base_segment, [10, 28, 17]), 3.5), 3, None)
urgent_ratio = np.clip(rng.normal(np.choose(base_segment, [0.10, 0.18, 0.42]), 0.06), 0, 0.8)
support_hours = np.clip(rng.normal(np.choose(base_segment, [8, 16, 34]), 4), 1, None)
customers = pd.DataFrame({
"customer_id": [f"C{i:03d}" for i in range(1, n_customer + 1)],
"annual_sales_million_yen": annual_sales,
"order_frequency": order_frequency,
"urgent_order_ratio": urgent_ratio,
"support_hours": support_hours,
})
cluster_features = customers.columns[1:]
X_scaled = StandardScaler().fit_transform(customers[cluster_features])
customers["cluster"] = KMeans(n_clusters=3, random_state=42, n_init=20).fit_predict(X_scaled)
profile_063 = customers.groupby("cluster")[cluster_features].mean()
profile_063["customers"] = customers.groupby("cluster").size()
display(profile_063.round(2))
fig, ax = plt.subplots(figsize=(7, 4))
for cluster, group in customers.groupby("cluster"):
ax.scatter(group["annual_sales_million_yen"], group["support_hours"], label=f"Cluster {cluster}", alpha=0.75)
ax.set_title("法人顧客の購買規模と技術支援負荷")
ax.set_xlabel("年間売上（百万円）")
ax.set_ylabel("年間技術支援時間")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| annual_sales_million_yen | order_frequency | urgent_order_ratio | support_hours | customers | |
|---|---|---|---|---|---|
| cluster | |||||
| 0 | 93.73 | 17.02 | 0.38 | 34.58 | 25 |
| 1 | 66.72 | 9.13 | 0.10 | 8.12 | 28 |
| 2 | 152.65 | 28.49 | 0.19 | 15.65 | 22 |
結果の読み取り
プロファイル表から、「高頻度・大口」「緊急対応・支援負荷大」「標準対応」のように意味付けできます。散布図で重なりがあるのは自然であり、クラスタを絶対的な顧客分類にしてはいけません。営業・サービス担当者が例外顧客を確認し、半年ごとに移動を追う用途が現実的です。
No.064：次元削減 — 多変量センサーの変化を二つの軸で見る
実務での意味
設備には温度、振動、電流、音圧など多数の信号があります。個別グラフを並べるだけでは同時変化を捉えにくいため、主成分分析（PCA）で共通変動を少数の軸へ要約します。
分析・モデル化の考え方
PCA は標準化済みデータの分散が最大となる直交方向を順に求めます。第1主成分は で、 は分散を最大化する単位ベクトルです。寄与率は情報量の目安ですが、主成分に現れない局所異常もあるため、元信号を捨てる処理ではありません。
Pythonで確認する
m = 260
load = rng.uniform(0.3, 1.0, m)
wear = np.linspace(0, 1, m) + rng.normal(0, 0.08, m)
sensors = pd.DataFrame({
"temperature": 48 + 22 * load + 8 * wear + rng.normal(0, 1.5, m),
"vibration": 1.0 + 1.2 * load + 1.5 * wear + rng.normal(0, 0.15, m),
"current": 10 + 9 * load + rng.normal(0, 0.8, m),
"sound": 62 + 8 * load + 5 * wear + rng.normal(0, 1.2, m),
"oil_particles": 18 + 30 * wear + rng.normal(0, 3, m),
"pressure": 4.0 + 2.8 * load - 0.5 * wear + rng.normal(0, 0.25, m),
})
pca_064 = PCA(n_components=2).fit(StandardScaler().fit_transform(sensors))
score_064 = pca_064.transform(StandardScaler().fit_transform(sensors))
loading_064 = pd.DataFrame(pca_064.components_.T, index=sensors.columns, columns=["PC1", "PC2"])
display(loading_064.round(3))
print("累積寄与率:", np.cumsum(pca_064.explained_variance_ratio_).round(3))
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].bar(["PC1", "PC2"], pca_064.explained_variance_ratio_)
axes[0].set_title("主成分の寄与率")
axes[0].set_xlabel("主成分")
axes[0].set_ylabel("寄与率")
axes[0].grid(True, axis="y", alpha=0.3)
sc = axes[1].scatter(score_064[:, 0], score_064[:, 1], c=wear, cmap="viridis", alpha=0.75)
axes[1].set_title("設備状態の主成分得点")
axes[1].set_xlabel("第1主成分")
axes[1].set_ylabel("第2主成分")
axes[1].grid(True, alpha=0.3)
fig.colorbar(sc, ax=axes[1], label="摩耗度（架空）")
plt.tight_layout()
plt.show()| PC1 | PC2 | |
|---|---|---|
| temperature | 0.490 | 0.068 |
| vibration | 0.448 | -0.327 |
| current | 0.389 | 0.431 |
| sound | 0.466 | -0.101 |
| oil_particles | 0.323 | -0.575 |
| pressure | 0.296 | 0.601 |
累積寄与率: [0.633 0.902]
結果の読み取り
寄与率は二軸で全変動のどの程度を表せるか、負荷量表は各軸をどのセンサーが構成するかを示します。得点が摩耗に沿って移動するなら、点検トリガーの補助指標にできます。ただし主成分の符号は任意であり、「PC1が高ければ必ず危険」と固定せず、保全記録と対応付けて解釈します。
No.065：異常検知 — 故障ラベルが少ない設備を監視する
実務での意味
重大故障は件数が少なく、教師あり学習に十分な正解ラベルを集めにくいのが普通です。異常検知は、通常と異なる観測を点検候補として順位付けし、状態基準保全を支援します。
分析・モデル化の考え方
Isolation Forest は、ランダムな特徴と分割点で観測を孤立させ、少ない分割で孤立する点ほど異常とみなします。contamination は期待する異常比率であり、真の故障率そのものではありません。スコアは故障確率ではなく、通常状態からの相対的な外れ度です。
Pythonで確認する
monitor = sensors.copy()
monitor["timestamp"] = pd.date_range("2026-01-01", periods=len(monitor), freq="h")
anomaly_idx = np.array([68, 151, 218, 242])
monitor.loc[anomaly_idx, "vibration"] += [2.2, 2.5, 2.8, 2.4]
monitor.loc[anomaly_idx, "temperature"] += [12, 15, 11, 14]
features_065 = list(sensors.columns)
model_065 = IsolationForest(n_estimators=250, contamination=0.03, random_state=42)
monitor["anomaly_label"] = model_065.fit_predict(monitor[features_065])
monitor["anomaly_score"] = -model_065.score_samples(monitor[features_065])
display(monitor.nlargest(8, "anomaly_score")[["timestamp", "vibration", "temperature", "anomaly_score"]].round(3))
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(monitor["timestamp"], monitor["anomaly_score"], label="異常スコア")
flag = monitor["anomaly_label"] == -1
ax.scatter(monitor.loc[flag, "timestamp"], monitor.loc[flag, "anomaly_score"], color="red", label="点検候補", zorder=3)
ax.set_title("設備センサーの異常スコア推移")
ax.set_xlabel("時刻")
ax.set_ylabel("異常スコア（大きいほど異常）")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| timestamp | vibration | temperature | anomaly_score | |
|---|---|---|---|---|
| 242 | 2026-01-11 02:00:00 | 5.392 | 78.802 | 0.639 |
| 68 | 2026-01-03 20:00:00 | 4.714 | 81.011 | 0.633 |
| 218 | 2026-01-10 02:00:00 | 5.692 | 79.421 | 0.629 |
| 151 | 2026-01-07 07:00:00 | 5.178 | 78.464 | 0.615 |
| 16 | 2026-01-01 16:00:00 | 1.670 | 52.597 | 0.604 |
| 10 | 2026-01-01 10:00:00 | 1.299 | 53.733 | 0.589 |
| 1 | 2026-01-01 01:00:00 | 2.039 | 67.991 | 0.572 |
| 230 | 2026-01-10 14:00:00 | 3.797 | 74.965 | 0.567 |
結果の読み取り
赤点は停止命令ではなく、保全員が波形・加工条件・直前作業を確認する点検候補です。上位表に同じ時間帯が連続するなら一つの事象としてまとめます。誤報履歴と実故障履歴を蓄積し、設備・運転モード別にしきい値を更新する運用が必要です。
No.066：時系列予測 — 季節性を含む補修部品需要を見通す
実務での意味
補修部品の需要予測は、生産枠、材料発注、安全在庫の基礎です。ランダム分割では未来の観測が学習へ混ざるため、過去から未来へ進む順序を守って評価します。
分析・モデル化の考え方
ここではトレンドと週次・月次の周期を Fourier 特徴で表し、線形回帰します。
最後の30日をテスト期間とし、その後30日を将来予測します。実務では販促、休業日、大口案件など、予測時点で既知の説明変数も追加します。
Pythonで確認する
days = 240
t = np.arange(days)
demand = (82 + 0.09 * t + 13 * np.sin(2 * np.pi * t / 7)
+ 8 * np.sin(2 * np.pi * t / 30) + rng.normal(0, 6, days))
demand = np.clip(np.rint(demand), 0, None)
demand_df = pd.DataFrame({"date": pd.date_range("2025-09-01", periods=days, freq="D"), "demand": demand})
def calendar_features(index):
index = np.asarray(index)
return np.column_stack([
index,
np.sin(2 * np.pi * index / 7), np.cos(2 * np.pi * index / 7),
np.sin(2 * np.pi * index / 30), np.cos(2 * np.pi * index / 30),
])
train_end = days - 30
model_066 = LinearRegression().fit(calendar_features(t[:train_end]), demand[:train_end])
test_pred = model_066.predict(calendar_features(t[train_end:]))
future_t = np.arange(days, days + 30)
future_pred = model_066.predict(calendar_features(future_t))
print(f"直近30日の時系列MAE: {mean_absolute_error(demand[train_end:], test_pred):.2f} 個/日")
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(demand_df["date"].iloc[-90:], demand[-90:], label="実績")
ax.plot(demand_df["date"].iloc[train_end:], test_pred, label="テスト予測", linestyle="--")
future_date = pd.date_range(demand_df["date"].iloc[-1] + pd.Timedelta(days=1), periods=30, freq="D")
ax.plot(future_date, future_pred, label="将来30日予測", color="red")
ax.set_title("補修部品の日次需要予測")
ax.set_xlabel("日付")
ax.set_ylabel("需要（個/日）")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()直近30日の時系列MAE: 4.30 個/日
結果の読み取り
時系列MAEは直近30日を未知期間として測った誤差です。将来線は平均的な予測であり、上振れリスクを表しません。調達・在庫には残差分布から予測区間を作り、欠品費用と在庫費用に応じて上側分位を使います。構造変化が起きた場合に備え、毎月「予測時点・予測値・実績値」を保存します。
No.067：因果推論 — 新治具の純粋な時間短縮効果を推定する
実務での意味
改善治具が難しい案件へ優先配備されると、導入群の加工時間は高く見えます。導入前後や導入・未導入の単純平均差は、案件構成の違い（交絡）を施策効果と取り違える可能性があります。
分析・モデル化の考え方
傾向スコア を推定し、逆確率重み付け（IPW）で比較可能な疑似母集団を作ります。
必要な仮定は、測定した交絡因子で条件付ければ割付が独立であること、両群が存在すること、処置定義が一貫することです。観測されない熟練度などは補正できません。
Pythonで確認する
n_067 = 700
complexity_067 = rng.uniform(1, 10, n_067)
lot_067 = rng.integers(30, 150, n_067)
assign_logit = -2.0 + 0.40 * complexity_067 + 0.004 * lot_067
prop_true = 1 / (1 + np.exp(-assign_logit))
treatment = rng.binomial(1, prop_true)
true_effect = -11.0
outcome = 45 + 8.0 * complexity_067 + 0.30 * lot_067 + true_effect * treatment + rng.normal(0, 8, n_067)
causal = pd.DataFrame({"complexity": complexity_067, "lot_size": lot_067, "new_jig": treatment, "cycle_time": outcome})
naive = causal.loc[causal.new_jig == 1, "cycle_time"].mean() - causal.loc[causal.new_jig == 0, "cycle_time"].mean()
ps_model = make_pipeline(StandardScaler(), LogisticRegression()).fit(causal[["complexity", "lot_size"]], treatment)
ps = np.clip(ps_model.predict_proba(causal[["complexity", "lot_size"]])[:, 1], 0.05, 0.95)
treated_mean = np.sum(treatment * outcome / ps) / np.sum(treatment / ps)
control_mean = np.sum((1-treatment) * outcome / (1-ps)) / np.sum((1-treatment) / (1-ps))
ipw_ate = treated_mean - control_mean
result_067 = pd.DataFrame({
"推定方法": ["単純平均差", "IPW", "データ生成上の真値"],
"時間短縮効果（分）": [naive, ipw_ate, true_effect],
})
display(result_067.round(2))
fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(ps[treatment == 0], bins=15, alpha=0.6, label="未導入")
ax.hist(ps[treatment == 1], bins=15, alpha=0.6, label="導入")
ax.set_title("新治具導入の傾向スコア分布")
ax.set_xlabel("導入確率（傾向スコア）")
ax.set_ylabel("案件数")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 推定方法 | 時間短縮効果（分） | |
|---|---|---|
| 0 | 単純平均差 | 8.69 |
| 1 | IPW | -12.05 |
| 2 | データ生成上の真値 | -11.00 |
結果の読み取り
難案件ほど導入されるため、単純比較は真の短縮効果を過小評価しやすく、IPWで真値に近づきます。傾向スコア分布の重なりが乏しい領域は、データから比較できません。可能なら段階導入や無作為化比較を設計し、安全・品質に悪影響がないことも同時に評価します。
No.068：ベイズ最適化 — 少ない試作で加工条件を探索する
実務での意味
切削速度などの条件探索は、試作ごとに材料・設備時間・検査費用がかかります。ベイズ最適化は、少数の評価結果から未知条件の性能を推定し、「良さそう」かつ「不確実」な条件を次に試します。
分析・モデル化の考え方
ガウス過程で目的関数の平均 と標準偏差 を表し、Expected Improvement（EI）を最大化する候補を選びます。
ここでは表面粗さとサイクル時間を合わせた損失を最小化します。実務では工具寿命、品質規格、設備上限を制約として明示します。
Pythonで確認する
def machining_loss(speed, feed):
base = 1.2 + ((speed - 185) / 38) ** 2 + ((feed - 0.20) / 0.055) ** 2
interaction = 0.35 * np.sin(speed / 18) * np.cos(feed * 30)
return base + interaction
grid_speed = np.linspace(120, 240, 45)
grid_feed = np.linspace(0.10, 0.30, 45)
S, F = np.meshgrid(grid_speed, grid_feed)
grid_X = np.column_stack([S.ravel(), F.ravel()])
scale = np.array([120.0, 0.20])
X_obs = np.array([[130, 0.12], [155, 0.26], [205, 0.14], [230, 0.28]], dtype=float)
y_obs = np.array([machining_loss(*x) for x in X_obs]) + rng.normal(0, 0.04, len(X_obs))
for _ in range(10):
kernel = ConstantKernel(1.0) * Matern(length_scale=[0.3, 0.3], nu=2.5) + WhiteKernel(0.002)
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True, random_state=42, n_restarts_optimizer=1)
gp.fit(X_obs / scale, y_obs)
mu, sigma = gp.predict(grid_X / scale, return_std=True)
improvement = y_obs.min() - mu
z = improvement / np.maximum(sigma, 1e-9)
ei = improvement * norm.cdf(z) + sigma * norm.pdf(z)
next_x = grid_X[np.argmax(ei)]
X_obs = np.vstack([X_obs, next_x])
y_obs = np.append(y_obs, machining_loss(*next_x) + rng.normal(0, 0.04))
best_idx = np.argmin(y_obs)
best_068 = pd.DataFrame({"切削速度": [X_obs[best_idx, 0]], "送り量": [X_obs[best_idx, 1]], "観測損失": [y_obs[best_idx]]})
display(best_068.round(3))
fig, ax = plt.subplots(figsize=(8, 5))
contour = ax.contourf(S, F, machining_loss(S, F), levels=18, cmap="viridis")
ax.scatter(X_obs[:4, 0], X_obs[:4, 1], color="white", edgecolor="black", label="初期条件")
ax.scatter(X_obs[4:, 0], X_obs[4:, 1], color="red", s=28, label="逐次探索")
ax.scatter(X_obs[best_idx, 0], X_obs[best_idx, 1], marker="*", s=180, color="gold", edgecolor="black", label="観測最良")
ax.set_title("ベイズ最適化による加工条件探索")
ax.set_xlabel("切削速度（m/min）")
ax.set_ylabel("送り量（mm/rev）")
ax.grid(True, alpha=0.25)
ax.legend()
fig.colorbar(contour, ax=ax, label="加工損失（小さいほど良い）")
plt.tight_layout()
plt.show()| 切削速度 | 送り量 | 観測損失 | |
|---|---|---|---|
| 0 | 188.182 | 0.2 | 0.882 |
結果の読み取り
赤点が逐次選ばれた試作条件で、星印が観測上の最良条件です。探索点が最良付近へ集まりつつ、不確実な領域も試していることが分かります。この結果をそのまま標準条件にせず、再現試験、品質規格、工具摩耗、安全限界を確認し、探索した範囲の外へ外挿しないことが重要です。
No.069：強化学習 — 劣化状態に応じた保全方針を学ぶ
実務での意味
予防保全を早く行いすぎると部品寿命を使い切れず、遅らせすぎると故障損失が増えます。強化学習は、現在の劣化状態と行動、その後の状態・費用を繰り返し評価し、長期費用を考慮した方策を学びます。
分析・モデル化の考え方
状態を健全度0〜4、行動を「運転継続」「保全」とします。Q-learning は次式で行動価値を更新します。
報酬設計が結論を決めるため、保全費、故障費、生産利益を事業数値で定義する必要があります。実設備へ直接探索させず、まずシミュレータと過去データで検証します。
Pythonで確認する
n_states, n_actions = 5, 2
Q = np.zeros((n_states, n_actions))
alpha, gamma, epsilon = 0.12, 0.96, 0.15
episode_returns = []
def step_069(state, action, random_gen):
if action == 1: # 保全
return 0, -18.0
if state == 4: # 故障状態で継続
return 0, -85.0
deteriorate_prob = 0.18 + 0.12 * state
next_state = min(4, state + int(random_gen.random() < deteriorate_prob))
reward = 8.0 - 1.5 * state
if next_state == 4:
reward -= 55.0
return next_state, reward
for episode in range(5000):
state, total = 0, 0.0
for _ in range(80):
if rng.random() < epsilon:
action = rng.integers(n_actions)
else:
action = int(np.argmax(Q[state]))
next_state, reward = step_069(state, action, rng)
Q[state, action] += alpha * (reward + gamma * Q[next_state].max() - Q[state, action])
state, total = next_state, total + reward
episode_returns.append(total)
policy = np.argmax(Q, axis=1)
policy_069 = pd.DataFrame({
"劣化状態": np.arange(n_states),
"継続Q値": Q[:, 0], "保全Q値": Q[:, 1],
"推奨行動": np.where(policy == 0, "運転継続", "保全"),
})
display(policy_069.round(2))
rolling = pd.Series(episode_returns).rolling(200).mean()
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(rolling, color="tab:blue")
ax.set_title("保全方策学習の平均累積報酬")
ax.set_xlabel("学習エピソード")
ax.set_ylabel("直近200回の平均累積報酬")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| 劣化状態 | 継続Q値 | 保全Q値 | 推奨行動 | |
|---|---|---|---|---|
| 0 | 0 | 125.80 | 105.74 | 運転継続 |
| 1 | 1 | 112.67 | 104.70 | 運転継続 |
| 2 | 2 | 106.98 | 105.37 | 運転継続 |
| 3 | 3 | 64.38 | 105.16 | 保全 |
| 4 | 4 | 37.63 | 105.25 | 保全 |
結果の読み取り
Q値表は状態ごとの長期価値と推奨行動を示し、劣化が進んだ段階で保全へ切り替わります。切替点は、設定した故障費・保全費・遷移確率に依存します。現場導入では、既存ルールとシミュレーション方策をオフライン比較し、安全制約で許可された行動だけを提示する意思決定支援から始めます。
No.070：グラフ機械学習 — サプライヤー網の近傍リスクを集約する
実務での意味
供給停止リスクは各社の財務・品質だけでなく、同じ物流拠点や二次サプライヤーへの依存を通じて波及します。表形式では失われる「誰とつながっているか」を特徴量に加えるのがグラフ機械学習の基本です。
分析・モデル化の考え方
各ノード の特徴 に、隣接ノードの平均特徴を集約します。
これはグラフニューラルネットワークの message passing を単純化した一段です。ロジスティック回帰でリスクを推定し、個社特徴のみの場合と比較します。取引がないこととデータ欠損を区別する必要があります。
Pythonで確認する
n_nodes = 90
rng_070 = np.random.default_rng(42)
G = nx.stochastic_block_model([30, 30, 30], [[0.20, 0.02, 0.01], [0.02, 0.18, 0.03], [0.01, 0.03, 0.22]], seed=42)
community = np.repeat(np.arange(3), 30)
quality_issue = np.clip(rng_070.beta(2, 8, n_nodes) + 0.12 * (community == 2), 0, 1)
delivery_delay = np.clip(rng_070.beta(2, 7, n_nodes) + 0.12 * (community == 1), 0, 1)
financial_stress = np.clip(rng_070.beta(2, 6, n_nodes), 0, 1)
base_risk = 1.4 * quality_issue + 1.2 * delivery_delay + 1.0 * financial_stress
neighbor_base = np.array([np.mean([base_risk[u] for u in G.neighbors(v)]) for v in G.nodes])
risk_prob = 1 / (1 + np.exp(-(-2.0 + 0.3 * base_risk + 2.2 * neighbor_base)))
risk_label = rng_070.binomial(1, risk_prob)
own_X = np.column_stack([quality_issue, delivery_delay, financial_stress])
neighbor_X = np.array([own_X[list(G.neighbors(v))].mean(axis=0) for v in G.nodes])
graph_X = np.column_stack([own_X, neighbor_X])
train_idx, test_idx = train_test_split(np.arange(n_nodes), test_size=0.35, random_state=42, stratify=risk_label)
own_model = make_pipeline(StandardScaler(), LogisticRegression(class_weight="balanced", random_state=42)).fit(own_X[train_idx], risk_label[train_idx])
graph_model = make_pipeline(StandardScaler(), LogisticRegression(class_weight="balanced", random_state=42)).fit(graph_X[train_idx], risk_label[train_idx])
own_prob = own_model.predict_proba(own_X[test_idx])[:, 1]
graph_prob = graph_model.predict_proba(graph_X[test_idx])[:, 1]
comparison_070 = pd.DataFrame({
"モデル": ["個社特徴のみ", "個社 + 近傍集約"],
"テストROC-AUC": [roc_auc_score(risk_label[test_idx], own_prob), roc_auc_score(risk_label[test_idx], graph_prob)],
})
display(comparison_070.round(3))
all_graph_prob = graph_model.predict_proba(graph_X)[:, 1]
pos = nx.spring_layout(G, seed=42)
fig, ax = plt.subplots(figsize=(8, 6))
nodes = nx.draw_networkx_nodes(G, pos, node_color=all_graph_prob, cmap="Reds", node_size=110, ax=ax)
nx.draw_networkx_edges(G, pos, alpha=0.22, ax=ax)
ax.set_title("サプライヤー取引網と推定供給リスク")
ax.set_xlabel("ネットワーク配置 X（表示用）")
ax.set_ylabel("ネットワーク配置 Y（表示用）")
ax.grid(True, alpha=0.15)
fig.colorbar(nodes, ax=ax, label="推定リスク確率")
plt.tight_layout()
plt.show()| モデル | テストROC-AUC | |
|---|---|---|
| 0 | 個社特徴のみ | 0.496 |
| 1 | 個社 + 近傍集約 | 0.609 |
結果の読み取り
近傍集約を加えたモデルが常に勝つとは限らず、テストAUCの差はネットワーク情報の追加価値を検証する材料です。赤いノードが集中する領域では、個社監査だけでなく共通物流・二次取引先・代替調達を確認します。予測確率だけで発注停止せず、説明可能なリスク要因と調達担当者の確認を組み合わせます。
対象ノックを通して見える実務上の示唆
- 予測対象より意思決定を先に決める：工数予測なら納期回答、不良分類なら追加検査というように、利用者・タイミング・選択肢を定義します。
- 分割方法と評価指標を業務に合わせる：時系列は未来をテストにし、不均衡分類は再現率と誤判定費用を見ます。平均精度だけでは不十分です。
- 相関・予測・因果を区別する：予測精度が高くても、条件を変えた効果は分かりません。改善効果には実験設計か交絡調整が必要です。
- 点ではなく関係と時間を扱う：時系列、強化学習、グラフは、独立した一行の表では表しにくい依存関係をモデルへ入れます。
- 人が介入できる出力にする：確率、異常スコア、推奨方策には、しきい値、根拠、例外処理、最終権限を設定します。
実務導入する場合に必要なこと
1. 業務定義とベースライン
モデル導入前に、現行ルール、判断頻度、担当者、誤判断費用を記録します。高度なモデルだけでなく、「品種別平均」「前週値」「現在の点検基準」と比較し、追加価値を測ります。
2. データ品質と時点整合
設備・品種・ロット・取引先IDを統一し、欠損理由とセンサー校正履歴を残します。特徴量が予測時点で利用可能か、後から修正された実績値が混ざっていないかを検査します。
3. 検証、安全、責任分界
シャドー運用で推奨と実判断を比較し、品種・設備・期間別に性能を確認します。安全・品質に関わる判断は制約ルールを優先し、自動化範囲、承認者、停止・手動復帰手順を明文化します。
4. 継続監視
入力分布、予測精度、介入率、事業KPIを監視し、学習データと本番の差を検知します。再学習の周期だけでなく、モデルを据え置く・停止する条件も決めます。
まとめ
No.061〜No.070では、加工時間と不良の予測、顧客・設備状態の構造把握、需要と異常の時間的監視、改善効果の因果推定、試作・保全の逐次意思決定、取引ネットワークのリスク評価を扱いました。
製造業で機械学習を定着させる鍵は、複雑なモデルそのものではなく、予測時点を守ったデータ、費用に沿った評価、現場の制約、人の最終判断、実績フィードバックです。まず限定工程でベースラインを作り、意思決定の質が改善したかを測りながら適用範囲を広げます。
法人向けのご相談
数理工房では、製造業における機械学習・統計解析のテーマ選定から、PoC、データ基盤、業務実装、精度・効果の継続監視、内製化支援までご相談を承ります。
- 品質予測・設備異常検知・需要予測の設計と検証
- 改善施策の実験設計・因果効果評価
- 加工条件探索・保全方策・サプライチェーンリスク分析
- 現場で運用できるKPI、説明画面、承認フローの設計
📩 お問い合わせ: surikobo.co.jp/contact まずはお気軽にご相談ください。