品質・生産リスクを数字で捉える：製造業のためのNumPy・確率シミュレーション10本ノック

本記事は、確率・統計をPythonで実装する100本ノックの第1回です。全100本では、NumPyの基礎、確率分布、可視化、統計量、推定、仮説検定、回帰、ベイズ統計、時系列、シミュレーションを段階的に扱います。目的は文法や公式の暗記ではなく、品質、生産能力、在庫、設備保全といった製造業の判断を、再現可能な計算へ落とし込めるようになることです。

今回はNo.001〜No.010として、架空の精密部品工場を題材に、NumPy配列からモンテカルロ法・計算速度比較までを扱います。同じデータを段階的に読み替えることで、配列計算が「集計を速くする道具」にとどまらず、不確実性を含む意思決定の土台になることを確認します。

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

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

精密部品工場では、ライン別の生産数、検査測定値、不良件数、停止時間が日々蓄積されます。ところが、会議で必要なのは単なる実績表ではありません。「規格中心からどれだけずれているか」「翌月の利益計画をどの程度の確率で達成できるか」「同じ分析を大量データにも適用できるか」という判断です。

本Notebookでは、3ライン・30日分の生産実績と、各ラインから得た寸法測定値を生成します。配列の形を明示し、要素ごとの計算と行列演算を使い分け、乱数による将来シナリオを再現可能に作るところまでを一続きで実装します。

現場でよくある状況

• ラインごとに生産量が違うのに、不良「件数」だけで良し悪しを比較している
• 製品別の規格値を各測定値へ適用する処理が、表計算のコピー＆ペーストになっている
• 乱数シミュレーションを実行するたび結果が変わり、レビュー時に数字を再現できない
• 平均値は報告されるが、分布の裾や「規格以下となる確率」が共有されない
• ループ処理が長時間化し、日次分析やシナリオ比較の運用が定着しない

これらは別々の問題に見えますが、データを適切な配列として設計し、計算の軸と乱数生成条件を明示することで、かなりの部分を共通の方法で整理できます。

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

製造データには、日、ライン、製品、測定項目など複数の軸があります。配列の軸を取り違えると、コードは動いても意味の違う集計になります。また、将来の生産量や不良数は確定値ではないため、平均シナリオだけでは下振れリスクを把握できません。

さらに、シミュレーション結果は乱数の標本です。試行回数、乱数シード、仮定した分布を記録しなければ、結果の再現も妥当性の検証もできません。本Notebookでは、各配列の shape、集計軸、単位、乱数生成器を意識的に表示します。

今回扱うノックの全体像

No.	テーマ	製造業での問い
001	NumPy配列の作成	日×ラインの実績をどう表現するか
002	ベクトル演算	生産規模の違いをならして品質を比較できるか
003	行列演算	日次実績を週次・ライン別へ一括集計できるか
004	ブロードキャスト	ライン別の基準値を全日に安全に適用できるか
005	乱数生成	未知の測定値や停止時間をどう模擬するか
006	乱数シード	分析結果を第三者が再現できるか
007	ヒストグラム描画	平均値の背後にあるばらつきや裾を読めるか
008	累積分布関数の描画	規格値以下となる割合を直接読めるか
009	モンテカルロシミュレーション	月間利益の下振れ確率を見積もれるか
010	計算速度比較	大量シナリオを運用可能な時間で計算できるか

No.001〜No.004はデータ構造と確定的な計算、No.005〜No.009は不確実性の生成・可視化・評価、No.010は継続運用に必要な計算効率を扱います。

Python環境の準備

外部データには依存せず、numpy、pandas、matplotlibを使います。日本語ラベルを表示するため japanize_matplotlib を読み込みます。乱数にはNumPyの Generator を使い、基準シードを1か所で管理します。

import sys
import timeit

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

SEED = 20260711
rng = np.random.default_rng(SEED)
np.set_printoptions(precision=3, suppress=True)
pd.options.display.float_format = "{:,.3f}".format

print(f"Python: {sys.version.split()[0]}")
print(f"NumPy: {np.__version__}")
print(f"乱数シード: {SEED}")

Python: 3.13.1
NumPy: 2.5.1
乱数シード: 20260711

架空データの作成

工場にはA・B・Cの3ラインがあり、30日間稼働したとします。日次生産数は曜日とライン能力で変動し、不良件数は生産数とライン別の基礎不良率から生成します。また、寸法規格が10.00 mmの部品について、各ライン500個ずつ測定した値を作ります。

この生成モデルは現実の因果構造を主張するものではありません。計算方法を比較するための共通シナリオであり、実務適用時には設備履歴、品種、ロット、材料、測定器などを含めて仮定を再設計します。

days = np.arange(1, 31)
lines = np.array(["A", "B", "C"])
weekday_factor = np.array([1.00 if (d - 1) % 7 < 5 else 0.82 for d in days])
line_capacity = np.array([1_050, 930, 810])

expected_output = weekday_factor[:, None] * line_capacity[None, :]
production = rng.poisson(expected_output)
base_defect_rate = np.array([0.018, 0.026, 0.034])
defects = rng.binomial(production, base_defect_rate)

dimension_mean = np.array([10.000, 10.012, 9.992])
dimension_sd = np.array([0.032, 0.041, 0.052])
measurements = rng.normal(
    loc=dimension_mean[:, None],
    scale=dimension_sd[:, None],
    size=(3, 500),
)

daily_df = pd.DataFrame(production, columns=lines, index=pd.Index(days, name="日"))
print("生産数配列:", production.shape, "不良数配列:", defects.shape,
      "寸法測定配列:", measurements.shape)
display(daily_df.head())

生産数配列: (30, 3) 不良数配列: (30, 3) 寸法測定配列: (3, 500)

	A	B	C
日
1	1019	933	845
2	1008	942	837
3	1057	941	800
4	1092	976	769
5	1070	967	831

No.001：NumPy配列の作成

実務での意味

日×ラインの数表を2次元配列にすると、「行は日、列はライン」というデータ契約を計算へ持ち込めます。製造データでは値だけでなく、軸の意味と単位を固定することが重要です。ここが曖昧だと、ライン合計と日合計を取り違える事故につながります。

分析・モデル化の考え方

NumPy配列は同じデータ型の値を規則的な形で保持します。shape=(30, 3) は30日×3ラインを意味します。axis=0 で集計すれば30日を畳み込んでライン別、axis=1 なら3ラインを畳み込んで日別の値が得られます。実務では、コードコメントやデータ辞書に各軸と単位を残します。

Pythonで確認する

array_design = pd.DataFrame({
    "確認項目": ["次元数", "形状", "要素数", "データ型", "先頭日の3ライン"],
    "値": [production.ndim, str(production.shape), production.size,
           str(production.dtype), production[0].tolist()],
})
display(array_design)
print("1日目・ラインBの生産数:", production[0, 1], "個")

	確認項目	値
0	次元数	2
1	形状	(30, 3)
2	要素数	90
3	データ型	int64
4	先頭日の3ライン	[1019, 933, 845]

1日目・ラインBの生産数: 933 個

結果の読み取り

生産数は30行×3列、合計90要素の整数配列です。production[0, 1] のように対象を明示すれば、1日目・ラインBという位置を一意に取得できます。実システムでは配列単体だけで軸ラベルを失わないよう、入力段階ではDataFrameやメタデータを併用し、数値計算部分でNumPyへ変換する設計が安全です。

No.002：ベクトル演算

実務での意味

不良件数の多いラインが、必ずしも品質の悪いラインとは限りません。生産数が多ければ不良件数も増えやすいため、ライン別不良率と品質損失額を同じ計算規則で比較します。

分析・モデル化の考え方

ライン別生産数をベクトル $\mathbf{n}$、不良数を $\mathbf{d}$ とすると、不良率ベクトルは要素ごとの除算

$$\mathbf{r}=\mathbf{d}\oslash\mathbf{n}$$

で求めます。さらに不良1個当たり損失を $\mathbf{c}$ とすれば、ライン別損失は $\mathbf{d}\odot\mathbf{c}$ です。ここで必要なのは内積ではなく、同じ位置の要素同士を計算する要素演算です。

Pythonで確認する

line_production = production.sum(axis=0)
line_defects = defects.sum(axis=0)
line_defect_rate = line_defects / line_production
loss_per_defect = np.array([4_500, 5_000, 5_800])
quality_loss = line_defects * loss_per_defect

line_kpi = pd.DataFrame({
    "生産数": line_production,
    "不良数": line_defects,
    "不良率": line_defect_rate,
    "不良1個当たり損失(円)": loss_per_defect,
    "推定品質損失(円)": quality_loss,
}, index=pd.Index(lines, name="ライン"))
display(line_kpi.style.format({"不良率": "{:.2%}", "推定品質損失(円)": "{:,.0f}"}))

	生産数	不良数	不良率	不良1個当たり損失(円)	推定品質損失(円)
ライン
A	29709	531	1.79%	4500	2,389,500
B	26751	710	2.65%	5000	3,550,000
C	23240	824	3.55%	5800	4,779,200

結果の読み取り

ベクトル演算により、3ラインへ同じKPI定義を一括適用できました。改善優先順位は不良率だけでなく、生産量と損失単価を反映した品質損失も合わせて判断します。ただし、ここでの損失単価は一定という仮定です。顧客流出、納期遅延、選別工数などを含める場合は、損失モデルを別途精緻化する必要があります。

No.003：行列演算

実務での意味

月次データから週次・ライン別集計を毎回手作業で作ると、期間の切り方や数式のコピー範囲が揺れます。集計ルールを行列として定義すると、複数ラインへ同じ期間変換を一括適用できます。

分析・モデル化の考え方

日次生産行列を $X\in\mathbb{R}^{30\times3}$、各日が第何週に属するかを表す集計行列を $W\in\{0,1\}^{5\times30}$ とします。週次生産行列は

$$Y=WX\in\mathbb{R}^{5\times3}$$

です。行列積では内側の次元30が一致し、日次の軸が集約されます。行列積 @ と要素積 * は意味が異なるため、モデル式に対応させて選びます。

Pythonで確認する

week_id = np.minimum((days - 1) // 7, 4)
week_matrix = np.eye(5, dtype=int)[week_id].T
weekly_production = week_matrix @ production

weekly_df = pd.DataFrame(
    weekly_production,
    index=pd.Index([f"第{i}週" for i in range(1, 6)], name="期間"),
    columns=lines,
)
weekly_df["全ライン"] = weekly_df.sum(axis=1)
display(weekly_df)
print("集計行列 @ 日次行列:", week_matrix.shape, "@", production.shape,
      "=", weekly_production.shape)

	A	B	C	全ライン
期間
第1週	6976	6340	5423	18739
第2週	6843	6187	5391	18421
第3週	6856	6277	5525	18658
第4週	6970	6102	5281	18353
第5週	2064	1845	1620	5529

集計行列 @ 日次行列: (5, 30) @ (30, 3) = (5, 3)

結果の読み取り

5×30の集計行列と30×3の日次行列から、5×3の週次行列が得られました。第5週は29〜30日の2日分だけなので、他週より小さい点に注意が必要です。単純な週合計を能力差として比較せず、稼働日数や計画時間で正規化して解釈します。

No.004：ブロードキャスト

実務での意味

ラインごとに計画生産数や管理基準が異なる場合、30日×3ラインの実績に対して、3個の基準値を適用したいことがあります。ブロードキャストを使えば、基準値のコピーを明示的に作らず計算できます。

分析・モデル化の考え方

形状 $(30,3)$ の実績行列から形状 $(3,)$ のライン別基準ベクトルを引くと、NumPyは末尾の次元を合わせ、基準を各行へ仮想的に展開します。計画達成率は

$$a_{ij}=\frac{x_{ij}}{t_j}$$

です。ブロードキャストは便利ですが、意図しない形状でも計算が成立する場合があります。shape と数件の手計算を検証します。

Pythonで確認する

daily_target = np.array([1_000, 900, 780])
achievement = production / daily_target
shortfall = production - daily_target

broadcast_check = pd.DataFrame({
    "ライン": lines,
    "日次目標": daily_target,
    "月平均実績": production.mean(axis=0),
    "平均達成率": achievement.mean(axis=0),
    "目標未達日数": (shortfall < 0).sum(axis=0),
}).set_index("ライン")
display(broadcast_check.style.format({"月平均実績": "{:,.1f}", "平均達成率": "{:.1%}"}))
print("実績", production.shape, "/ 目標", daily_target.shape, "-> 達成率", achievement.shape)

	日次目標	月平均実績	平均達成率	目標未達日数
ライン
A	1000	990.3	99.0%	9
B	900	891.7	99.1%	8
C	780	774.7	99.3%	12

実績 (30, 3) / 目標 (3,) -> 達成率 (30, 3)

結果の読み取り

1本のライン別目標ベクトルが全30日へ適用され、ライン別の平均達成率と未達日数を計算できました。ただし、休日にも平日と同じ目標を置けば未達日数を過大評価します。実務ではカレンダー、計画停止、品種構成を基準配列へ反映し、「比較可能な日だけを比較する」ことが必要です。

No.005：乱数生成

実務での意味

新条件での測定値や設備停止時間は、実施前には確定していません。乱数は未来を当てる道具ではなく、仮定したばらつきの下で起こり得るシナリオを多数作り、判断の感度を調べる道具です。

分析・モデル化の考え方

寸法値には正規分布、ある時間内の不良件数には二項分布など、変数の生成過程に合う分布を選びます。ここでは改善後のラインCの寸法を

$$X\sim\mathcal{N}(\mu=10.000,\sigma=0.038)$$

と仮定し、1,000個を生成します。分布選択は観測データ、工程知識、残差診断に基づき、外れ値やロット間差を無視していないか確認します。

Pythonで確認する

improved_c = rng.normal(loc=10.000, scale=0.038, size=1_000)
lower_spec, upper_spec = 9.90, 10.10
random_summary = pd.DataFrame({
    "指標": ["標本数", "平均(mm)", "標準偏差(mm)", "最小(mm)", "最大(mm)", "規格外率"],
    "値": [improved_c.size, improved_c.mean(), improved_c.std(ddof=1),
           improved_c.min(), improved_c.max(),
           ((improved_c < lower_spec) | (improved_c > upper_spec)).mean()],
})
display(random_summary)

	指標	値
0	標本数	1,000.000
1	平均(mm)	9.999
2	標準偏差(mm)	0.038
3	最小(mm)	9.890
4	最大(mm)	10.109
5	規格外率	0.007

結果の読み取り

生成標本の平均と標準偏差は設定値の近くになりますが、有限標本なので完全には一致しません。推定規格外率も試行ごとに変動します。改善案の評価では、平均を規格中心へ寄せる効果と、標準偏差を小さくする効果を分けて検討すると、設備調整とばらつき削減のどちらを優先すべきか説明しやすくなります。

No.006：乱数シード

実務での意味

経営会議、監査、モデルレビューでは、「同じ入力とコードから同じ結果が得られること」が重要です。シード固定により疑似乱数列を再現し、計算変更の影響だけを比較できます。

分析・モデル化の考え方

疑似乱数は決定的なアルゴリズムで生成されます。同じ生成器・同じシード・同じ呼び出し順なら同じ列になります。一方、シード固定はモデルの正しさを保証しません。また、1つのシードに都合のよい結論だけを依存させないよう、本番評価では複数シナリオや十分な試行回数を用います。

Pythonで確認する

rng_a = np.random.default_rng(SEED)
rng_b = np.random.default_rng(SEED)
rng_c = np.random.default_rng(SEED + 1)

sample_a = rng_a.integers(0, 100, size=8)
sample_b = rng_b.integers(0, 100, size=8)
sample_c = rng_c.integers(0, 100, size=8)
seed_check = pd.DataFrame({"同じseed_A": sample_a, "同じseed_B": sample_b,
                           "異なるseed": sample_c})
display(seed_check.T)
print("同じシードの配列は一致:", np.array_equal(sample_a, sample_b))
print("異なるシードの配列は一致:", np.array_equal(sample_a, sample_c))

	0	1	2	3	4	5	6	7
同じseed_A	84	19	17	90	7	53	61	66
同じseed_B	84	19	17	90	7	53	61	66
異なるseed	33	69	35	48	2	86	51	41

同じシードの配列は一致: True
異なるシードの配列は一致: False

結果の読み取り

同じシードから作った2列は完全に一致し、異なるシードの列は一致しません。実務ではシードに加え、ライブラリ版、分布パラメータ、試行回数、コード版を分析成果物へ記録します。並列処理を行う場合は、独立した乱数ストリームの設計も必要です。

No.007：ヒストグラム描画

実務での意味

平均寸法が規格中心に近くても、ばらつきが大きければ規格外が発生します。ヒストグラムは中心、広がり、歪み、複数の山、裾を可視化し、平均値だけでは見えない工程状態を確認する入口です。

分析・モデル化の考え方

ヒストグラムは値の範囲をビンに分け、各区間の度数を数えます。ビン幅で見え方が変わるため、比較するラインには共通のビンを使います。面積を1に正規化する density=True は分布形状の比較に便利ですが、実際の不良個数を示すものではありません。

Pythonで確認する

bins = np.linspace(9.82, 10.18, 28)
fig, ax = plt.subplots(figsize=(9, 5))
for i, line_name in enumerate(lines):
    ax.hist(measurements[i], bins=bins, alpha=0.45, density=True,
            label=f"ライン{line_name}")
ax.axvline(lower_spec, color="crimson", linestyle="--", label="規格下限")
ax.axvline(upper_spec, color="crimson", linestyle="--", label="規格上限")
ax.set_title("ライン別 寸法測定値の分布")
ax.set_xlabel("寸法 (mm)")
ax.set_ylabel("確率密度")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

hist_summary = pd.DataFrame({
    "平均(mm)": measurements.mean(axis=1),
    "標準偏差(mm)": measurements.std(axis=1, ddof=1),
    "標本規格外率": ((measurements < lower_spec) | (measurements > upper_spec)).mean(axis=1),
}, index=pd.Index(lines, name="ライン"))
display(hist_summary.style.format({"平均(mm)": "{:.4f}", "標準偏差(mm)": "{:.4f}",
                                   "標本規格外率": "{:.2%}"}))

	平均(mm)	標準偏差(mm)	標本規格外率
ライン
A	10.0001	0.0311	0.00%
B	10.0108	0.0406	1.60%
C	9.9906	0.0495	4.00%

結果の読み取り

ラインBは中心が上側へずれ、ラインCは分布の幅が相対的に広い設定です。この2つは対策が異なり、前者は中心調整、後者は材料・設備・測定系を含むばらつき要因の調査が候補になります。ヒストグラムだけで工程能力や分布適合性を断定せず、時系列、層別、管理図、測定システム解析も確認します。

No.008：累積分布関数の描画

実務での意味

「10.05 mm以下は何%か」「規格下限を下回る割合はどの程度か」のように、しきい値を使う意思決定では累積分布関数（CDF）が有効です。調達基準、選別基準、予防保全の警報値などの説明にも使えます。

分析・モデル化の考え方

累積分布関数は

$$F(x)=P(X\le x)$$

です。未知の母集団CDFに対し、標本を小さい順に並べ、各点までの割合を描く経験累積分布関数（ECDF）を使います。特定の分布形を仮定せず、観測標本に基づいてしきい値以下の割合を読めます。

Pythonで確認する

fig, ax = plt.subplots(figsize=(9, 5))
for i, line_name in enumerate(lines):
    x = np.sort(measurements[i])
    y = np.arange(1, x.size + 1) / x.size
    ax.step(x, y, where="post", label=f"ライン{line_name}")
ax.axvline(10.05, color="black", linestyle="--", label="判断しきい値 10.05 mm")
ax.set_title("ライン別 寸法測定値の経験累積分布（ECDF）")
ax.set_xlabel("寸法 (mm)")
ax.set_ylabel("累積確率 P(Xがx以下)")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

threshold = 10.05
cdf_at_threshold = (measurements <= threshold).mean(axis=1)
display(pd.DataFrame({"10.05mm以下の割合": cdf_at_threshold},
                     index=pd.Index(lines, name="ライン")).style.format("{:.1%}"))

	10.05mm以下の割合
ライン
A	95.0%
B	81.4%
C	88.0%

結果の読み取り

同じしきい値でも、ラインごとに累積割合が異なります。ECDFの立ち上がりが急なら値が狭い範囲へ集中し、緩やかならばらつきが大きいと読めます。上限超過率が欲しい場合は $1-F(x)$ を使いますが、標本点と同値を含むかどうかで境界の定義が変わるため、< と <= を規格定義に合わせます。

No.009：モンテカルロシミュレーション

実務での意味

月間利益は、生産量、不良、設備停止など複数の不確実性を同時に受けます。期待利益が計画を上回るだけでは、未達リスクは分かりません。多数の将来シナリオを生成し、下振れ確率と分位点を評価します。

分析・モデル化の考え方

各シナリオの利益を

$$\Pi = q(1-r)m - dc - h$$

とします。$q$ は生産数、$r$ は不良率、$m$ は良品1個当たり限界利益、$d$ は不良数、$c$ は不良損失、$h$ は停止損失です。ここでは説明を簡潔にするため独立な分布を仮定しますが、需要と稼働率、停止と不良率に相関がある場合は同時分布や条件付きモデルが必要です。

Pythonで確認する

mc_rng = np.random.default_rng(SEED + 9)
n_scenarios = 20_000
monthly_units = mc_rng.normal(line_production.sum(), 2_200, size=n_scenarios).clip(0)
scenario_defect_rate = mc_rng.beta(55, 1_945, size=n_scenarios)
scenario_defects = mc_rng.binomial(monthly_units.astype(int), scenario_defect_rate)
downtime_hours = mc_rng.gamma(shape=4.0, scale=5.0, size=n_scenarios)

margin_per_good = 620
loss_per_bad = 5_200
downtime_cost_per_hour = 85_000
profit = ((monthly_units - scenario_defects) * margin_per_good
          - scenario_defects * loss_per_bad
          - downtime_hours * downtime_cost_per_hour)
profit_million = profit / 1_000_000
profit_target = 33.0

mc_summary = pd.DataFrame({
    "指標": ["期待利益", "5%点", "中央値", "95%点", "33百万円未達確率"],
    "値": [profit_million.mean(), *np.quantile(profit_million, [0.05, 0.50, 0.95]),
           (profit_million < profit_target).mean()],
})
display(mc_summary)

fig, ax = plt.subplots(figsize=(9, 5))
ax.hist(profit_million, bins=45, color="steelblue", alpha=0.8)
ax.axvline(profit_target, color="crimson", linestyle="--", label="利益計画 33百万円")
ax.axvline(profit_million.mean(), color="black", linestyle=":", label="シナリオ平均")
ax.set_title("月間利益のモンテカルロシミュレーション")
ax.set_xlabel("月間利益 (百万円)")
ax.set_ylabel("シナリオ数")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	指標	値
0	期待利益	34.947
1	5%点	31.296
2	中央値	34.984
3	95%点	38.403
4	33百万円未達確率	0.181

結果の読み取り

平均だけでなく5%点と計画未達確率を示すことで、利益計画の安全余裕を議論できます。5%点は「必ず守れる下限」ではなく、仮定したモデル下で95%のシナリオが上回った値です。設備停止時間や不良率の分布が実態とずれれば結果も変わるため、過去データによる校正、悲観・標準・楽観条件の感度分析、モデル外リスクの定性評価を併記します。

No.010：計算速度比較

実務での意味

分析が正しくても、計算に数時間かかれば日次運用や条件探索に使えません。NumPyのベクトル化は、同じ処理を多数の要素へ適用するとき、Pythonの明示的なループより短く高速に書けることがあります。

分析・モデル化の考え方

ここでは100万個の測定値について、規格外率をループとベクトル演算で計算します。両者の答えが一致することを先に確認し、その後 timeit で複数回測定します。速度はCPU、ライブラリ版、配列サイズの影響を受けるため、絶対値ではなくこの実行環境での相対比較として読みます。

Pythonで確認する

speed_rng = np.random.default_rng(SEED + 10)
large_measurements = speed_rng.normal(10.0, 0.045, size=1_000_000)

def defect_rate_loop(values, lower, upper):
    count = 0
    for value in values:
        if value < lower or value > upper:
            count += 1
    return count / len(values)

def defect_rate_vectorized(values, lower, upper):
    return np.mean((values < lower) | (values > upper))

loop_result = defect_rate_loop(large_measurements, lower_spec, upper_spec)
vector_result = defect_rate_vectorized(large_measurements, lower_spec, upper_spec)
loop_time = timeit.timeit(
    lambda: defect_rate_loop(large_measurements, lower_spec, upper_spec), number=3
) / 3
vector_time = timeit.timeit(
    lambda: defect_rate_vectorized(large_measurements, lower_spec, upper_spec), number=10
) / 10

speed_comparison = pd.DataFrame({
    "方法": ["Pythonループ", "NumPyベクトル演算"],
    "規格外率": [loop_result, vector_result],
    "平均実行時間(秒)": [loop_time, vector_time],
    "ループ比": [1.0, vector_time / loop_time],
})
display(speed_comparison.style.format({"規格外率": "{:.4%}",
                                       "平均実行時間(秒)": "{:.6f}", "ループ比": "{:.3f}"}))
print("結果は一致:", np.isclose(loop_result, vector_result))
print(f"この環境ではベクトル演算が約 {loop_time / vector_time:.1f} 倍高速")

	方法	規格外率	平均実行時間(秒)	ループ比
0	Pythonループ	2.6316%	0.047111	1.000
1	NumPyベクトル演算	2.6316%	0.001195	0.025

結果は一致: True
この環境ではベクトル演算が約 39.4 倍高速

結果の読み取り

2つの方法は同じ規格外率を返し、この配列サイズではベクトル演算が大幅に高速です。性能改善では、まず答えの一致をテストし、実データに近いサイズで測定します。ベクトル化は一時配列によってメモリを消費するため、巨大データでは分割処理、集約方法、データ型、I/Oも含めて最適化します。

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

1. 配列設計は業務定義そのもの：日、ライン、品種、測定項目のどれを各軸に置くかを明示すると、集計ミスを減らせます。
2. 件数と率を使い分ける：品質状態の比較には率、経済影響の評価には件数と損失単価が必要です。
3. 基準値はデータとして管理する：ブロードキャストで一括適用できても、休日・品種・設備状態を無視した基準では正しい判断になりません。
4. 平均から分布へ視点を広げる：ヒストグラムとCDFにより、中心だけでなく裾、しきい値超過、ばらつきを判断材料にできます。
5. 不確実性を確率で伝える：モンテカルロ法では期待値に加え、分位点と計画未達確率を提示できます。
6. 再現性と速度を運用要件に含める：シード、環境、コード版を記録し、処理時間を測定して初めて継続運用へ移せます。

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

• 製品ID、ロット、ライン、設備、品種、シフト、日時、測定器を追跡できるデータモデル
• 寸法、不良率、生産数、停止時間、損失額についての単位・分母・除外条件の統一
• 欠測、重複、測定器校正、時刻ずれ、計画停止を検知するデータ品質ルール
• 分布仮定とパラメータを定期的に更新する仕組み、および実績とのバックテスト
• シード、Python環境、入力スナップショット、コード版、承認履歴を残す再現性管理
• 現場の管理限界、品質保証上の規格、経営上のリスク許容度を分けた判断基準
• 分析結果を誰が確認し、誰が設備条件や生産計画を変更するかという業務フロー

PoCでは、まず意思決定を1つに絞り、現行判断、必要なリードタイム、誤判断の損失、利用可能データを整理します。モデル精度だけでなく、改善行動につながったかを評価指標に含めることが重要です。

まとめ

No.001〜No.010では、製造データをNumPy配列として設計し、ベクトル・行列・ブロードキャストで集計し、乱数、ヒストグラム、CDF、モンテカルロ法で不確実性を評価しました。最後に、同じ規格外率計算でもベクトル化によって処理時間を短縮できることを確認しました。

大切なのは、コードが動くことだけではありません。配列の軸、数値の単位、乱数の仮定、しきい値の定義、結果が再現できる条件を説明できて初めて、分析は製造現場の意思決定に使える形になります。

法人向けのご相談

数理工房では、製造業の品質管理、工程能力評価、異常検知、需要・生産シミュレーション、数理最適化、データ分析人材育成をご支援しています。現場データの棚卸しからPoC、意思決定プロセスへの組み込み、社内研修まで、課題とデータ成熟度に合わせて設計します。

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