趣旨
小売業では今なおエクセルで在庫管理をしているのが大半だと思います。かつ在庫管理には一定の勘と経験が役立つ事があって、高齢化社会が始まっている昨今ではその勘と経験を持つ方々の退職により、ある種の"伝承"が途切れようとしています。
それに対処する術としては、流行りの生成AI関連のソフトウェアが雨後の筍の如く開発されています。しかしながらある程度高い導入コストと学習コストを必要とします。小売業には導入300万円もするようなソフトウェアを簡単に承認されません。（ざっくり3000万円～4000万円の売上を立てなければなりませんからね）
以降で説明する内容は、そんな導入コストを必要とせず学習コストのみで導入が出来る方法の一つである、統計学問のCox比例ハザードモデルを使った在庫損失確率（Stock out probability）と安全在庫の設計について紹介します。
Cox比例ハザードモデルとは
元来は医療統計で盛んに活用されてきたモデルの一つです。例えば様々な投薬や環境に応じて人が何日経過するとどれくらいの確率で死亡するのか、などに用いられてきました。これを 商品Aが何日経過するとどれくらいの確率でゼロになるのか？ という課題に応用します。
1. なにを解くモデルか（要点）
「あるイベントが“いつ”起きるか（time-to-event）」を、共変量（特徴量）を使って確率的に説明・予測する半パラメトリックモデルです。イベントまでの時間を扱い、**打ち切り（観測終了までイベントが起きない）**を自然に取り込めるのが強みです。
2. このモデルの何が面白いのか
Coxは「時間 × 打ち切り × 説明変数」がある場面なら幅広く使えます。医療以外の実務用途を、イベント・時間軸・共変量に分けて具体化すると導入が容易です。ポイントは、イベント定義（何を「起きた」とみなすか）と打ち切り（観測終了時点で未発生をどう扱うか）を最初に固定し、現場データの粒度に合わせた共変量設計を行うことです。
-
ドメイン：サプライチェーン/在庫
- イベント：ストックアウト発生
- 時間：補充/期首から欠品までの日数
- よく効く共変量例：需要水準・変動係数、リードタイム、安全在庫比率、仕入先OTIF、需要予測バイアス、プロモ有無
-
ドメイン：HR/人事
- イベント：退職
- 時間：入社から退職まで
- よく効く共変量例：部署、評価、異動、残業、マネジャー変更
3. 各種関数と解釈
ハザード関数（瞬間リスク） を次のように定義します。
$$
h(t \mid \mathbf{x}) = h_0(t),\exp\big(\mathbf{x}^\top \boldsymbol{\beta}\big)
$$
- $h_0(t)$：ベースラインハザード（時間依存、分布は仮定しない＝半パラメトリック）
- $\mathbf{x}$：共変量ベクトル、$\boldsymbol{\beta}$：係数
ハザード比（HR） 共変量 $x_k$ を1だけ増やしたときの倍率を次のように定義します。
$$
\mathrm{HR}_k = \exp(\beta_k)
$$
例えば $\exp(\beta_k)=1.20$ は「瞬間リスクが 20% 増」、$0.80$ は「20% 減」となります。
生存関数（イベント未発生確率）
在庫がゼロにならない確率は以下の数式で表現します。
$$
S(t \mid \mathbf{x}) = \exp\Big(-H_0(t),\exp\big(\mathbf{x}^\top \boldsymbol{\beta}\big)\Big),
\quad
H_0(t) = \int_0^t h_0(u),du
$$
部分尤度（推定式）
イベントが起きた観測集合を $\mathcal{E}$、時点 $t_i$ のリスク集合を $\mathcal{R}(t_i)$ とすると，
$$L(\beta)=\prod_{i\in E}\frac{\exp(x_i^T\beta)}{\sum_{j\in R(t_i)}\exp(x_j^T\beta)}$$
$$
\mathcal{E}:\ \text{イベント発生観測の集合},\quad
\mathcal{R}(t_i):\ \text{時点 } t_i \text{ のリスク集合}
$$
積から和に変換するために対数尤度に変換します。
$$\ell(\beta)=\sum_{i\in E}\Big[x_i^T\beta-\log\Big(\sum_{j\in R(t_i)}\exp(x_j^T\beta)\Big)\Big]$$
この尤度から $h_0(t)$ は消えるため，分布仮定なしに $\boldsymbol{\beta}$ を推定できます。
比例ハザード（PH）仮定
任意の2個体 A, B に対してハザード比は時間によらず一定であることが前提となっています。
$$\frac{h(t|x_A)}{h(t|x_B)}=\exp\big((x_A - x_B)^T\beta\big)$$
4. 在庫最適化へのつなげ方（安全在庫を“目標サービス水準から逆算”）
在庫領域では、「期間 L 日で欠品しない確率」$S(L\mid x)$ が目標 $S^*$（例：0.95）以上になるように安全在庫を決めます。Cox では生存関数は
$$S(L\mid x)=\exp\big(-H_0(L),\exp(x^T\beta)\big)$$
と書けます。ここで $H_0(L)$ はベースライン累積ハザード、$x^T\beta$ は線形予測子です。安全在庫の操作変数を $b$（例：buffer_ratio）とし、他の要因を固定した線形予測子を $\eta(b)$ と置くと、目標を満たす条件は
$$S^* \le \exp\big(-H_0(L),\exp(\eta(b))\big)$$
これを $\eta(b)$ について解くと
$$\eta(b) \le \log\Big(\frac{-\log S^*}{H_0(L)}\Big)$$
となります。実務では $H_0(L)$ も推定量で、$\eta(b)$ は $b$ に対して単調（$b$ を増やすと $S$ が増える）であることが多いので、$b$ の二分探索で最小の $b$ を求めるのが堅実です（後述する付属の recommend_buffer_ratio() が実装例です）
Pythonによるコーディング
ここからは、実際に模擬データを作成するところから記述をしていきます。
ライブラリインストールとユーティリティ
python バージョンは3.11.9です。
import math
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from lifelines import CoxPHFitter, KaplanMeierFitter
from statsmodels.nonparametric.smoothers_lowess import lowess
def _set_seed(seed: int | None = 42) -> None:
if seed is not None:
np.random.seed(seed)
模擬データ生成
今回の共変量は仮で作成しています。1行が商品AのSKU x 1拠点の1スパンの観測結果というデータになっています。n = 1000ならば商品Aの販売拠点が1000か所あるという事ですね。100均一の商品がイメージに近いです。
一旦Cox比例ハザードモデルの線形予測子からβ真値を設定し、各変量は各変量のドメイン知識および変数の特徴からランダム確率分布を選定しています。そこから指数分布の親戚であるワイブル分布を使って潜在的な欠品時間を計算し、打ち切り時間よりも欠品時間のほうが短ければ欠品というフラグを立てることとします。
def generate_inventory_survival_data(n: int = 1000, seed: int | None = 42) -> pd.DataFrame:
"""
在庫損失（ストックアウト）を対象にしたサバイバル用模擬データを生成。
変数:
- duration_days: 観測開始から欠品までの期間（日）。打ち切りあり
- event : 1=欠品発生, 0=打ち切り
- mean_demand : 平均出荷数量(単位/日)
- demand_cv : 需要の変動係数 [0.1, 1.0], 変動係数 = 標準偏差/平均, 学習用の簡易版
- lead_time_days: 平均リードタイム（日）, (サプライヤー調達+輸送+入庫)
- buffer_ratio : 安全在庫の相対量, リードタイム期間の平均出荷数量に対して安全在庫がどれだけ積まれているか（例: safety_stock / (mean_demand * lead_time_days))
- otif : 仕入先OTIF(On Time In Full) [0,1]
- forecast_bias: 需要予測バイアス（+で不足方向）
データ生成ロジック:
- Weibull(shape=k, scale=lambda)から潜在時間T~Weibullを発生
- 共変量によりscaleをexp(-linear/shape)で調整
- 検閲時間C~Uniform(60, 365)日とし、duration=min(T, C), event=1[T<=C]
注: 値域や係数は学習用に設定しており、実務では確認/校正が必要
"""
_set_seed(seed)
N = n
# 需要・LT・OTIF等の分布設定（例示）
mean_demand = np.random.lognormal(mean=math.log(8.0), sigma=0.5, size=N) # ~ 3〜30程度
demand_cv = np.random.beta(a=2.5, b=5, size=N)
lead_time_days = np.random.lognormal(mean=math.log(20.0), sigma=0.4, size=N).astype(int) # ~ 10〜40程度
otif = np.clip(np.random.beta(a=18, b=2, size=N), 0.7, 0.999) # 多くが0.9台、下限0.7
forecast_bias = np.random.normal(loc=0.0, scale=0.2, size=N) # ±20%程度のバイアス
# 安全在庫の相対量（例: safety_stock / (mean_demand * lead_time)）。負を許容してクリップ。
buffer_ratio = np.clip(np.random.normal(loc=0.3, scale=0.15, size=N), -0.5, 2.0)
# 線形予測子（欠品リスクを上げる/下げる方向に係数設定）
# スケールの違いを考慮して係数は控えめに（後でCoxで回収可能）。
beta_mean_demand = 0.015 # 大きい需要ほど欠品しやすい
beta_demand_cv = 0.8 # 変動が大きいほど欠品しやすい
beta_lead_time = 0.02 # LTが長いほど欠品しやすい
beta_buffer = -1.2 # バッファが大きいほど欠品しにくい
beta_otif = -1.0 # OTIF高いほど欠品しにくい
beta_bias = 0.6 # 予測が不足方向（+）ほど欠品しやすい
linear = (
beta_mean_demand * mean_demand
+ beta_demand_cv * demand_cv
+ beta_lead_time * lead_time_days
+ beta_buffer * buffer_ratio
+ beta_otif * otif
+ beta_bias * forecast_bias
)
# Weibullのshape/scale
shape = 1.3 # 時間とともにリスクがやや増加
base_scale = 180.0 # ベースのスケール（日）
adjusted_scale = base_scale * np.exp(-linear / shape)
# 潜在欠品時間Tと打ち切り時間C
T = np.random.weibull(shape, size=N) * adjusted_scale
C = np.random.uniform(60.0, 365.0, size=N)
duration = np.minimum(T, C).astype(int)
event = (T <= C).astype(int)
df = pd.DataFrame(
{
"duration_days": duration,
"event": event,
"mean_demand": mean_demand,
"demand_cv": demand_cv,
"lead_time_days": lead_time_days,
"buffer_ratio": buffer_ratio,
"otif": otif,
"forecast_bias": forecast_bias,
}
)
return df
Cox比例ハザードモデルの一番良いところは、打ち切りデータも扱えるという点です。すべての行データが全く同じ時間で測定できていることは実社会においては稀なほうです。これをいくつかの仮定を置いた状態でモデル化できます。
加えて、このデータテーブルと同じようにデータを準備していく際に、やや面倒なのがforecast_biasの部分です。今回は 観測を始める前の例えば過去12週間分のforecast_bias を表現しています。実環境だと、都度一定期間におけるforecastを出して予実管理をしますが、このデータテーブルでは一つのスカラー値でeventが1になるのかを見ています。
(Option)実務でのforecast_biasの算出の仕方
最も使いやすいのは 符号付き平均パーセント誤差（MPE） のローリング版です。ここではコードと整合性を合わせるため 正の値＝過小予測（不足見積もり） に統一します。
週次バケット例（SKU×拠点ごと）
期間 $W$（直近8〜12週など）で、“1週先”に対して当時のスナップショットで予測した値を使います。
$$\mathrm{forecast_bias}=\frac{\sum_{t\in W}\big(A_t - F_{t\mid t-1}\big)}{\sum_{t\in W}\max(A_t,\varepsilon)}$$
- $F_{t\mid t-1}$：週 $t-1$ 時点のスナップショットが示す 翌週 $t$ の予測
- $\varepsilon$：ゼロ除算回避用の小値（例：0.1）
代替例（APICS系）
「APICS流のバイアス」は、予測が実績に対してどちらに寄って外れているかを、かんたんに数で表したものです。考え方は次の通りです。
- A … 実績（その期間に実際に起きた需要・受注・出荷など）
- F … その当時に立てていた 予測（スナップショットを使うのがポイント）
まず、予測誤差を足し合わせた RSFE（Running Sum of Forecast Errors） を作ります。
RSFE が プラスに大きいほど「実績 > 予測」＝いつも控えめに見積もっている（過小予測）、
マイナスに大きいほど「実績 < 予測」＝**いつも盛って見積もっている（過大予測）**という結果を表現します。
$$\mathrm{RSFE}=\sum(A-F)$$
この RSFE を実績の合計で割ったものが Forecast Bias です。値のスケールをそろえるための正規化だと考えてください。
$$\mathrm{Forecast\ Bias}=\frac{\mathrm{RSFE}}{\sum A}$$
- 正の値 → 過小予測（足りなく見積もっている傾向）
- 負の値 → 過大予測（多めに見積もっている傾向）
さらに、どれくらい“偏りが強い”のかを見たいときは 追跡信号（Tracking Signal） を使います。
誤差の“ふつうの大きさ”＝ MAD（平均絶対誤差） を分母にして、偏りを相対化します。
$$\mathrm{MAD}=\frac{1}{n}\sum |A-F|,\qquad \mathrm{Tracking\ Signal}=\frac{\mathrm{RSFE}}{\mathrm{MAD}}$$
- 直感：TS の絶対値が大きいほど、誤差が一方向に偏って続いている合図
- 運用：しきい値（例：±3〜4 など）はチームで決め、超えたら原因を点検（新製品・販促・受給逼迫など）
作成したデータをpairplotで確認
# データ生成（既存の関数を利用）
df = generate_inventory_survival_data(n=1000, seed=42)
# イベントのラベル化（任意）
df["event_label"] = df["event"].map({1: "Stockout", 0: "Censored"})
# 色の設定。2クラス用に取得して辞書に割り当て
base_palette = sns.color_palette("muted", n_colors=2)
palette = {"Stockout": base_palette[0], "Censored": base_palette[1]}
sns.set_theme(style="whitegrid", context="notebook")
# 可視化する変数を指定（必要に応じて調整）
vars_to_plot = [
"duration_days",
"mean_demand",
"demand_cv",
"lead_time_days",
"buffer_ratio",
"otif",
"forecast_bias",
]
g = sns.pairplot(
df,
vars=vars_to_plot,
hue="event_label", # クラス別に色分け
palette=palette,
corner=False, # 両方表示
diag_kind="hist", # 対角はヒストグラム
plot_kws=dict(alpha=0.5, s=14, edgecolor="none"), # 重なり対策
diag_kws=dict(edgecolor="white", linewidth=0.5),
)
g.figure.suptitle("Inventory Survival Data — Pairplot", y=1.02, fontsize=12)
plt.show()
在庫予測を実際に行ってみる
データフレームの値は真値に対してランダム生成しているので収束しますが、実務の場合だと収束しなかったり分散が大きすぎてp値が0.05を下回らなかったりしますので、要注意です。
def perform_cox_inventory_analysis(df: pd.DataFrame, penalizer: float = 0.0) -> CoxPHFitter:
"""
在庫欠品サバイバルデータにCoxPHを適合させ、主要指標を表示して返す。
"""
required = {"duration_days", "event"}
missing = required - set(df.columns)
if missing:
raise KeyError(f"missing columns: {missing}")
print("=" * 70)
print("Cox比例ハザードモデル[在庫欠品]")
print("=" * 70)
print(f"サンプル数 : {len(df):,}")
print(f"イベント発生数 : {df['event'].sum():,} ({df['event'].mean():.1%})")
print(f"平均観察期間(日) : {df['duration_days'].mean():.1f}")
cph = CoxPHFitter(penalizer=penalizer)
covariates = [
"mean_demand",
"demand_cv",
"lead_time_days",
"buffer_ratio",
"otif",
"forecast_bias",
]
cph.fit(df[["duration_days", "event"] + covariates], duration_col="duration_days", event_col="event")
print("\n-- 係数とHR --")
print(cph.summary[["coef", "exp(coef)", "exp(coef) lower 95%", "exp(coef) upper 95%", "p"]].round(4))
print("\nConcordance Index:", f"{cph.concordance_index_:.4f}")
print("Partial AIC :", f"{cph.AIC_partial_:.2f}")
return cph
Concordance IndexがCox比例ハザードモデルの評価指標の一つとなります。行からランダムに二つ選んでそれぞれのイベント発生が観測値に対して早いのか遅いのかを0から1の範囲で表現しています。ドメインによりますが0.6であれば十分当たっているという意味付けができます。0.75を超えてくると実用に耐えうる精度の良いモデルと言われています。
偏回帰係数βを確認
def plot_hazard_ratios(cph: CoxPHFitter) -> None:
summary = cph.summary
vars_ = summary.index.to_list()
hr = summary["exp(coef)"].values
lo = summary["exp(coef) lower 95%"].values
hi = summary["exp(coef) upper 95%"].values
y = np.arange(len(vars_))
plt.figure(figsize=(7, 4.5))
plt.errorbar(hr, y, xerr=[hr - lo, hi - hr], fmt="o", capsize=4)
plt.axvline(1.0, linestyle="--")
plt.yticks(y, vars_)
plt.xlabel("Hazard Ratio")
plt.title("Hazard Ratios (95% CI)")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
Hazard Ratioはコードの通りexp(β) です。つまり1を超えるか超えないかは欠品が発生する倍率が1を超えるか超えないかと同じと解釈できます。
これにより、需要(demand_cv)が1単位増加するとおおよそ1.5 ~ 4倍ほど欠品が起きる確率が増加することを意味しています。おなじくforcast_biasは過小評価（思っていたよりも需要が多かった≒＋側）ほど欠品しやすくなっていることが数字で表現できています。
一方、otif(On time in full:顧客の希望日までに指定された場所へ、製品をどの程度適切に全量納品するか)などは、1単位増加すればするほど欠品が起きる倍率は1よりも小さくなっていることが表現できています。
サバイバル曲線の可視化
def km_by_buffer_tertile(df: pd.DataFrame) -> None:
df = df.copy()
# buffer_ratioの三分位で層別
df["buffer_bucket"] = pd.qcut(df["buffer_ratio"], q=3, labels=["Low", "Mid", "High"])
kmf = KaplanMeierFitter()
plt.figure(figsize=(7, 4.5))
for label in ["Low", "Mid", "High"]:
sub = df[df["buffer_bucket"] == label]
kmf.fit(sub["duration_days"], event_observed=sub["event"], label=f"buffer {label}")
kmf.plot_survival_function()
plt.xlabel("Days")
plt.ylabel("Survival: P(no stockout by t)")
plt.title("KM survival by buffer tertile")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
バッファはもともと連続値としていましたが、閾値を設けてカテゴリ化し、それぞれのバッファレベルに応じてサバイバル関数の曲線がどのよう変化するのかを可視化しています。
日数が経過するにつれて、その時間までに欠品なしの確率（no stockout by t)を表しています。ここが医療統計だと患者が生きている確率（死亡というイベントが発生しない確率）となります。これが応用の部分の 「在庫管理でやりたかった事」 の一つに直結しています。つまり、縦軸が0.5を下回る時、buffer ratio（安全在庫がどれだけ積まれているか）のLowとHighの差が予測範囲こみで最大約120日ほど差が出ているということを予測しています。
さらに、縦軸が0.5を下回るようになると、buffer highの予測範囲（透過させている部分）は他の2区分と重なりが少なくなっていますが、MidとLowは非常に重なりが増えてきています。このことから、0.5を下回るようになるとバッファを高めに積んでおくと、ある程度Highの設定に信頼性が高そうですですが、MidとLowに関しては重なる部分が多く、MidとLowに十分有意な差があるとはいえず、どちらの設定値にしたとしても大した違いは出ないだろう、ということもわかります。
モデルのズレを可視化
予測結果が実際とどれくらいズレているか（残差がどれくらいあるか）を可視化して結果を確認します。
def residual_plot(cph: CoxPHFitter, df: pd.DataFrame) -> None:
res = cph.compute_residuals(df, kind="martingale")
ph = cph.predict_partial_hazard(df)
plt.figure(figsize=(7, 4.5))
plt.scatter(ph.values.ravel(), res["martingale"].values.ravel(), alpha=0.5)
plt.axhline(0.0, linestyle="--")
plt.xlabel("Predicted partial hazard")
plt.ylabel("Martingale residuals")
plt.title("Residual diagnostic")
plt.tight_layout()
plt.show()
横軸は、モデルが付けた相対リスクのスコアで、スコアが高くなるほどハザードが起きやすいことを表現しています。縦軸は、ある時間tまでに実際にイベントが発生したかどうか（欠品したかどうか）を0/1で表現しています。プロット一つが1行のデータです。
なかなか馴染みがないグラフですが、要は縦軸が「1」側に多すぎると予測が観測よりも早くイベント（欠品）が発生したことを意味しています。一方「0」側に多すぎると予測が観測よりも遅くイベント（欠品）が発生したことを意味しています。言い換えると、+側にプロットが多いと「まだまだ欠品しないでしょ！余裕余裕！」と思っていたら欠品したケース、－側にプロットが多いと「すぐ欠品するはずだ！危険！」と思っていたらまだ欠品しなかったケースを意味しています。
今回のケースを見るとプラスマイナス両方に上手く分散しており均等になっているなという印象を持てますが、横軸が0に近いほど縦軸がマイナス側へ大きく「0」から離れているプロットも見受けられるので、今回の予測モデルは、ハザードが起きにくい領域ほどやや悲観的なモデルと言えそうです。曲線はLOWESS（ローカル回帰）でざっくりいうと"動く平均線"という所ですが、これが部分ハザード（横軸）が0.5付近は縦軸がマイナス側に触れていることから、可視化したプロットの分布傾向の読み取りが、この曲線と一致していることが分かります。
(Option) 共変量一つ一つに対してマルチンゲール残差を確認する
# 共変量ごとに確認する場合の可視化
def residual_vs_covariate(cph, df, col):
y = cph.compute_residuals(df, kind="martingale")["martingale"].values.ravel()
x = df[col].values
xy = lowess(y, x, frac=0.3, return_sorted=True)
plt.figure(figsize=(7, 4.5))
plt.scatter(x, y, alpha=0.35, s=12)
plt.plot(xy[:,0], xy[:,1], linewidth=2)
plt.axhline(0, ls="--")
plt.xlabel(col)
plt.ylabel("Martingale residuals")
plt.title(f"Residuals vs {col}")
plt.tight_layout()
plt.show()
リードタイムのみを変化させてそれ以外を固定したときのリードタイムに対するマルチンゲール残差の予測と観測結果の差もプロット可能です。これを見るとリードタイムが短い部分はややマイナス側（悲観的）な予測モデルといえそうです。
モデル結果を使って、個別拠点（個別行）のデータから予測する
def predict_survival_at(cph: CoxPHFitter, row: pd.Series, t_days: float) -> float:
"""
観測個体(row)のt_daysにおける生存確率 S(t) を返す。
row: Series あるいは単一行DataFrameでも可
"""
X = pd.DataFrame([row]) if isinstance(row, pd.Series) else row
sf = cph.predict_survival_function(X, times=[t_days])
# timesが行index、個体が列: 単一個体なら[0,0]
return float(sf.iloc[0, 0])
t_days = 230.0
row_value = 110
predict_survival = predict_survival_at(cph=cph,row=df.loc[row_value,:],t_days=t_days)
display(df.loc[row_value,:])
print("生存確率（在庫が無くなっていない確率）: ",predict_survival)
この関数により、各行（各拠点）に対して起点からの日数に対してどれくらいの確率でまだ在庫が欠品していないかを確率で算出できます。※回帰結果のP値の良し悪しは忘れずに算出結果を吟味します。
じゃあ実際に安全在庫のバッファはいくらにすればいいの？
def recommend_buffer_ratio(
cph: CoxPHFitter,
row: pd.Series,
horizon_days: float,
target_survival: float = 0.95,
bounds: tuple = (-0.5, 2.0),
tol: float = 1e-3,
max_iter: int = 50,
) -> dict[str, float]:
"""
指定した期間(horizon_days)における「欠品なし確率」が target_survival 以上になる
最小の buffer_ratio を二分探索で求める。
戻り値: {"buffer_ratio": 推奨値, "survival": 達成S(t)}
備考:
- buffer_ratio を増やすと生存確率が単調増加することを前提（係数が負方向だから）
- 実務では buffer_ratio の物理定義に合わせ、最終的な安全在庫量(単位)へ換算すること
safety_stock_units ≈ buffer_ratio * mean_demand * lead_time_days
"""
lo, hi = bounds
# 単調性チェック（境界）
r = row.copy()
r["buffer_ratio"] = lo
s_lo = predict_survival_at(cph, r, horizon_days)
r["buffer_ratio"] = hi
s_hi = predict_survival_at(cph, r, horizon_days)
if s_hi < target_survival:
# 上限でも達成できない
return {"buffer_ratio": hi, "survival": s_hi}
if s_lo >= target_survival:
# 下限でも達成できる
return {"buffer_ratio": lo, "survival": s_lo}
for _ in range(max_iter):
mid = (lo + hi) / 2.0
r["buffer_ratio"] = mid
s_mid = predict_survival_at(cph, r, horizon_days)
if s_mid >= target_survival:
hi = mid
else:
lo = mid
if abs(hi - lo) < tol:
break
r["buffer_ratio"] = hi
s = predict_survival_at(cph, r, horizon_days)
return {"buffer_ratio": float(hi), "survival": float(s)}
buffer_ratio_condition_dict = recommend_buffer_ratio(
cph=cph,
row=df.loc[row_value,:],
horizon_days=t_days,
)
print(buffer_ratio_condition_dict["buffer_ratio"])
print(buffer_ratio_condition_dict["survival"])
この関数を用いて、指定した起点からの期間で欠品が起きない確率が95%以上（いったんは固定値）になる時の最も小さいバッファ率を算出します。
これで、過去の経験から統計的に算出されるバッファ率を用いて人の勘と経験を裏付けるための意思決定材料の一つを提案出来ます。
注意したいは、引数に割り当てるboundsの2値です。この2つの値は業務上の制約に相当します。このboundsの意味をもう一度整理すると、
bounds,lo = 0 : バッファ無し
bounds,lo = 1 : バッファはリードタイム期間中の平均需要と同じ量を持つ
という意味となります。
ですので上記関数は会社のルール上、バッファ率は少なくとも0.5を持っておく事、というルールが付けられた場合の関数になります。
一方hiのほうは、「運用上上限,かつ,統計上の目安,かつ,出力結果の数値が安定」の3つの条件付き数値といえます。
- 運用上上限：倉庫の空きスペースなど
- 統計上の目安：リードタイム中の需要は正規分布に従うと仮定し、z-score * 需要の変動係数 / sqrt(リードタイム日数）の式から、z-score = 1.96(97.5%) or 1.64(95%)を入れてみる
- 出力結果の数値が安定：実際に数値を変えながら結果を見比べて、survivalの値が例えば0.9に近くなるように設定しなおすなどの調整を自分で決める。
という3条件によって最終的に決定します。