需要予測と新聞売り子問題による在庫最適化

はじめに

本記事は NTTドコモ Advent Calendar 2024 の12日目の記事です。

こんにちは。NTTドコモサービスイノベーション部の淺田です。普段はデータ分析・AI技術を活用した業務効率化や意思決定支援を行っています。

今回のテーマ「需要予測と新聞売り子問題による在庫最適化」は、私が業務の中で取り組んでいる課題の一つである在庫の最適化を基にした内容です。

概要

現代のビジネス環境、特に小売・製造業において在庫管理は収益性を左右する重要な課題です。過剰な在庫は廃棄コストや保管コストを増大させ、一方で在庫の不足は機会損失を招きます。この問題に対処するため、本記事では以下の流れで在庫の最適化を行うアプローチを紹介します：

需要予測
- 過去の販売数データを用いて、将来の需要量を推定します。これには、通常時系列分析や機械学習モデルを活用します。
誤差分布の取得
- 需要予測と実際の販売数の差(誤差)を統計的に分析し、その分布を把握します。誤差分布は需要の不確実性を定量化する鍵です。
新聞売り子問題の応用
- 需要の不確実性を考慮した「新聞売り子問題」のモデルを適用し、粗利を最大化するための在庫を最適化します。このアプローチでは、販売価格、原価、欠品による機会損失、余剰在庫コストを考慮します。

最初に理論的な部分を説明し、後半で実際にPythonを用いて実装を行います。

対象読者

データサイエンティストやアナリスト
- 機械学習を活用した予測モデルに関心がある方。また、需要予測と在庫最適化の応用例を学びたい方。
サプライチェーンや在庫管理の担当者
- 日々の在庫管理業務に携わり、予測精度や効率化を改善したい方。

新聞売り子問題とは

概要

新聞売り子モデル(Newsvendor model)とは、古典的な確率的在庫モデルです。

新聞売り子モデルとは、需要量がとある分布に従う時、粗利の期待値を最大化するためには「在庫はどの程度であるべきか」という問題です。
- コストを最小化する考え方もありますが、同じ結論にたどり着きます。
- ここで、需要量(demand)とは「在庫が無限にあった場合の販売数」を指し、粗利(gross profit)とは「販売金額合計から原価合計を差し引いたもの」を指します。

機械学習モデルによる将来の需要推定量を用いる場合

機械学習モデルにより需要予測を行い、「需要量が需要推定量からどの程度ぶれるか」という誤差の分布が一定であると仮定した時、「需要推定量に対して在庫をどの程度多く用意しておくべきか」という値は実は解析的に求めることができます。

以下、数式を使って求めていきますが、気になる方以外は適宜読み飛ばしてください。

まず、変数を次のように定義しておきます。

$t$ : 何日かを表す添字
$D_{t}$ : $t$ 日の需要量(demand)。確率変数。
$\hat{D_{t}}$ : $t$ 日の需要量の推定値、需要推定量
$q_{t} \in \mathbb{Z}_{\geq 0}$ : $t$ 日の初期在庫数(quantity)
$s_{t} \in \mathbb{Z}_{\geq 0}$ : $t$ 日の販売数(sales)
$p \in \mathbb{R}$ : 1個あたりの売価(price)
$c \in \mathbb{R}$ : 1個あたりの原価(cost)
$\alpha_{t} = \frac{D_{t} - \hat{D_{t}}}{\hat{D_{t}}}$ : $t$ 日目の需要量が需要推定量からどの程度ぶれるかを表す量。確率変数。
$x_{t} = \frac{q_{t} - \hat{D_{t}}}{\hat{D_{t}}} \in \mathbb{R}$ : $t$ 日の需要推定量に対して在庫をどの程度多く用意しておくべきか
$E[\cdot]$ : 期待値
$F$ : $\alpha_{t}$ の累積分布関数
$F^{-1}$ : $F$ の逆関数

また、

$f(\alpha) = F'(\alpha)$
$F(\alpha) = \int_{-\infty}^{\alpha} f(u) du$
$f(\alpha) > 0$
$F(\infty) = 1$
$F(-\infty) = 0$

としておきます。

この時、

$\displaystyle{ \begin{aligned} D_{t} = \hat{D_{t}} ( \alpha_{t} + 1 ) \end{aligned} }$

$\displaystyle{ \begin{aligned} q_{t} &= \hat{D_{t}} \cdot \frac{q_{t}}{\hat{D_{t}}} \\ &= \hat{D_{t}} \cdot \left( \frac{q_{t} - \hat{D_{t}}}{\hat{D_{t}}} + 1 \right) \\ &= \hat{D_{t}} ( x_{t} + 1 ) \end{aligned} }$

が成り立っています。

販売数について、在庫以上には売れないため在庫と需要量の内小さい方の値を取ります。

$\displaystyle{ \begin{aligned} s_{t} &= \min \{ q_{t}, D_{t} \} \\ &= \min \{ \hat{D_{t}} ( x_{t} + 1 ), \hat{D_{t}} ( \alpha_{t} + 1 ) \} \\ &= \hat{D_{t}} ( \min \{ x_{t}, \alpha_{t} \} + 1 ) \end{aligned} }$

また、需要推定量 $\hat{D_{t}}$ が与えられた時

$\displaystyle{ \begin{aligned} E[\min \{ x_{t}, \alpha_{t} \} \mid \hat{D_{t}}] &= E[\alpha_{t} \mid \hat{D_{t}}, \alpha_{t} \leq x_{t}] F(x_{t}) + x_{t} (1 - F(x_{t})) \\ \\ &= \left( \frac{1}{ F(x_{t}) } \int_{u \leq x_{t}} u f(u) du \right) F(x_{t}) + x_{t} (1 - F(x_{t})) \\ \\ &= \int_{u \leq x_{t}} u f(u) du + x_{t} (1 - F(x_{t})) \\ \\ &= \int_{u \leq x_{t}} u f(u) du + x_{t} - x_{t} \int_{u \leq x_{t}} f(u) du \\ \\ &= \int_{u \leq x_{t}} ( u - x_{t} ) f(u) du + x_{t} \end{aligned} }$

より、販売数 $s_{t}$ の期待値は次式で与えられます。

$\displaystyle{ \begin{aligned} E[s_{t} \mid \hat{D_{t}}] &= E[\min \{ q_{t}, D_{t} \} \mid \hat{D_{t}}] \\ &= E[\hat{D_{t}} ( \min \{ x_{t}, \alpha_{t} \} + 1 ) \mid \hat{D_{t}}] \\ &= \hat{D_{t}} ( E[\min \{ x_{t}, \alpha_{t} \} \mid \hat{D_{t}}] + 1 ) \\ \\ &= \hat{D_{t}} \{ E[\alpha_{t} \mid \hat{D_{t}}, \alpha_{t} \leq x_{t}] F(x_{t}) + x_{t} (1 - F(x_{t})) + 1 \} \\ \\ &= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} \left( \frac{1}{ F(x_{t}) } \int_{u \leq x_{t}} u f(u) du \right) F(x_{t}) + x_{t} (1 - F(x_{t})) + 1 \end{aligned} \end{Bmatrix} \\ \\ &= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} \int_{u \leq x_{t}} u f(u) du + x_{t} (1 - F(x_{t})) + 1 \end{aligned} \end{Bmatrix} \\ \\ &= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} \int_{u \leq x_{t}} ( u - x_{t} ) f(u) du + x_{t} + 1 \end{aligned} \end{Bmatrix} \end{aligned} }$

次に、 $t$ 日の粗利(gross profit)を次式で定義します。

$\displaystyle{ \begin{aligned} \mathrm{profit}_{t} &:= p s_{t} - c q_{t} \\ &= p \hat{D_{t}} ( \min \{ x_{t}, \alpha_{t} \} + 1 ) - c \hat{D_{t}} ( x_{t} + 1 ) \\ &= \hat{D_{t}} ( p \cdot \min \{ x_{t}, \alpha_{t} \} + p - c x_{t} - c ) \\ &= \hat{D_{t}} ( p \cdot \min \{ x_{t}, \alpha_{t} \} - c x_{t} + p - c ) \end{aligned} }$

前述の通り需要が在庫を超えると販売数は一定になるため、需要と粗利の関係性は次のグラフの通りです。

この粗利の式は

$\displaystyle{ \begin{aligned} &\mathrm{profit}_{t} \\ &= (p - c) D_{t} - (p - c)(D_{t} - s_{t}) - c (q_{t} - s_{t}) = \text{差益額} \cdot \text{需要量}_{t} - \text{機会損失額}_{t} - \text{廃棄額}_{t} \\ &= (p - c) s_{t} - c (q_{t} - s_{t}) = \text{差益額} \cdot \text{販売量}_{t} - \text{廃棄額}_{t} \\ &= p s_{t} - c q_{t} = \text{売上合計}_{t} - \text{原価合計}_{t} \end{aligned} }$

と解釈することも可能です。

よってここまでの内容から、需要推定量 $\hat{D_{t}}$ が与えられた時の粗利(gross profit)の期待値は次のように計算することができます。

$\displaystyle{ \begin{aligned} E[\mathrm{profit}_{t} \mid \hat{D_{t}}] &= E[p s_{t} - c q_{t} \mid \hat{D_{t}}] \\ &= E[\hat{D_{t}} ( p \cdot \min \{ x_{t}, \alpha_{t} \} - c x_{t} + p - c ) \mid \hat{D_{t}}] \\ &= \hat{D_{t}} E[p \cdot \min \{ x_{t}, \alpha_{t} \} - c x_{t} + p - c \mid \hat{D_{t}}] \\ &= \hat{D_{t}} ( E[p \cdot \min \{ x_{t}, \alpha_{t} \} \mid \hat{D_{t}}] - E[c x_{t} \mid \hat{D_{t}}] + E[p \mid \hat{D_{t}}] - E[c \mid \hat{D_{t}}] ) \\ &= \hat{D_{t}} ( p \cdot E[\min \{ x_{t}, \alpha_{t} \} \mid \hat{D_{t}}] - c x_{t} + p - c ) \\ \\ &= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}} u f(u) du + p x_{t} (1 - F(x_{t})) - c x_{t} + p - c \end{aligned} \end{Bmatrix} \\ \\ &= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}} u f(u) du - p x_{t} F(x_{t}) + (p - c) x_{t} + p - c \end{aligned} \end{Bmatrix} \\ \\ &= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}} ( u - x_{t} ) f(u) du + p x_{t} - c x_{t} + p - c \end{aligned} \end{Bmatrix} \\ \\ &= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}} ( u - x_{t} ) f(u) du + (p - c)(x_{t} + 1) \end{aligned} \end{Bmatrix} \end{aligned} }$

今回は、この粗利期待値 $E[\mathrm{profit}_{t} \mid \hat{D_{t}}]$ を最大化する $x_{t}$ の算出を行います。

なお、これは次のような最適化問題として定式化することができます。

$\displaystyle{ \begin{aligned} \left \Vert \quad \begin{aligned} &\textbf{maximize} \quad && E[\mathrm{profit}_{t} \mid \hat{D_{t}}] \\ &&&= E[p s_{t} - c q_{t} \mid \hat{D_{t}}] \\ &&&= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}} u f(u) du - p x_{t} F(x_{t}) + (p - c) x_{t} + p - c \end{aligned} \end{Bmatrix} \\ &&&= \hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}} ( u - x_{t} ) f(u) du + (p - c)(x_{t} + 1) \end{aligned} \end{Bmatrix} \\ \\ &\textbf{subject to} \quad && x_{t} = \frac{q_{t} - \hat{D_{t}}}{\hat{D_{t}}} \in \mathbb{R} \end{aligned} \right. \end{aligned} }$

ここから、粗利期待値を最大化する $x_{t}$ を解析的に求めます。荒い議論になりますが、簡単に粗利期待値を $x_{t}$ で偏微分すると

$\displaystyle{ \begin{aligned} \dfrac{\partial}{\partial x_{t}} E[\mathrm{profit}_{t}] &= \dfrac{\partial}{\partial x_{t}} \hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}} u f(u) du - p x_{t} F(x_{t}) + (p - c) x_{t} + p - c \end{aligned} \end{Bmatrix} \\ \\ &= \hat{D_{t}} \{ p x_{t} f(x_{t}) - p F(x_{t}) - p x_{t} F' (x_{t}) + p - c \} \\ &= \hat{D_{t}} \{ p x_{t} f(x_{t}) - p F(x_{t}) - p x_{t} f(x_{t}) + p - c \} \\ &= \hat{D_{t}} \{ - p F(x_{t}) + p - c \} \end{aligned} }$

であるため、最適解があるとすれば少なくとも

$\displaystyle{ \begin{aligned} \dfrac{\partial}{\partial x_{t}} E[\mathrm{profit}_{t}] &= \hat{D_{t}} \{ - p F(x_{t}) + p - c \} \\ &= 0 \end{aligned} }$

であることが必要となります。これを変形すると、

$\displaystyle{ \begin{aligned} F (x_{t}) = \frac{p - c}{p} \in (0, 1) \end{aligned} }$

となります。ここで $F$ は0以上1以下の値を取り連続で単調増加するため、この式を満たす $x_{t}$ はただ1つ存在しそれを $x_{t}^*$ とすると、

$\displaystyle{ \begin{aligned} x_{t}^* &= F ^ {-1} \left( \frac{p - c}{p} \right) \end{aligned} }$

と書けます。 ( $F ^ {-1}$ は逆像ではなく, 逆写像)

増減表を考えれば粗利期待値は結局 $x_{t}^*$ で最大値を取るため、これが最適解だと分かります。

次に、ここまでに書いた内容を実際にPythonプログラムを用いてサンプルデータに適用し、粗利期待値が最大化されることを確認します。

Pythonプログラムによる在庫最適化の実装

問題設定

次のような問題を具体的に考えることにします。

前提
- ある商品について、日別の販売数データがあります。
- 全体で1200日間を考え、最初の1100日分を学習期間、最後の100日分を検証期間とします。
  - 「1100日分のデータが溜まっていて、今後100日の在庫を最適化したいような状況」を想定しています。
- 商品の販売価格は100円、原価は30円とします。
- 売れ残った商品は翌日には持ち越されず、廃棄されるとします。
問題
- 検証期間(最後の100日間)の 「粗利の期待値」を最大化する最適な在庫数量(100日分の計画)を考えます。

また、変数を「機械学習モデルによる将来の需要推定量を用いる場合」と同様に次のように定義しておきます。

$t$ : 何日かを表す添字
$D_{t}$ : $t$ 日の需要量(demand)。確率変数。
$\hat{D_{t}}$ : $t$ 日の需要量の推定値、需要推定量
$q_{t} \in \mathbb{Z}_{\geq 0}$ : $t$ 日の期初在庫数(quantity)
$s_{t} \in \mathbb{Z}_{\geq 0}$ : $t$ 日の販売数(sales)
$p = 100 \in \mathbb{R}$ : 1個あたりの売価(price)
$c = 30 \in \mathbb{R}$ : 1個あたりの原価(cost)
$\alpha_{t} = \frac{D_{t} - \hat{D_{t}}}{\hat{D_{t}}}$ : $t$ 日目の需要量が需要推定量からどの程度ぶれるかを表す量。確率変数。
$x_{t} = \frac{q_{t} - \hat{D_{t}}}{\hat{D_{t}}} \in \mathbb{R}$ : $t$ 日の需要推定量に対して在庫をどの程度多く用意しておくべきか
$E[\cdot]$ : 期待値
$F$ : $\alpha_{t}$ の累積分布関数
$F^{-1}$ : $F$ の逆関数

初期設定

使用するライブラリのimportと再現性を持たせるためseed値の設定をしておきます。

# Load modules
import numpy as np
import pandas as pd
import scipy

import matplotlib.pyplot as plt
from sklearn.metrics import mean_absolute_error as mae_val
from sklearn.metrics import mean_absolute_percentage_error as mape_val
from sklearn.metrics import mean_squared_error as mse_val

from prophet import Prophet
import japanize_matplotlib

# seed
np.random.seed(57)

サンプルデータの作成

平日は売上が高く土日は売上が低いような商品を想定し、日別の販売数サンプルデータを1200日分作成します。今回は120日周期の緩やかなサインカーブに従っています。

def calc_sample_base(row):
    """
    指定された行のデータから売上の基礎値を計算。

    Parameters:
    row (dict): 入力データを表す辞書

    Returns:
    dict: 基礎売上に関連する値を含む辞書。
        - 'base_sum' (float): 全体の基礎売上値。
        - 'base1' (int): 曜日に基づく基礎売上値。
        - 'base2' (float): 周期性（120日周期）による調整値。
    """
    day_num = int(row['date_idx'])
    weekday = row['weekday']
    base1 = 100
    
    # 土日は売上減少
    if weekday in ['Saturday', 'Sunday']:
        base1 -= 40

    # 長期的なトレンド
    base2 = 20 * np.sin(2 * np.pi * (day_num / 120))

    base_sum = base1 + base2

    ret = {
        'base_sum': base_sum,
        'base1': base1,
        'base2': base2
    }
    return ret

# 1200レコードのサンプルデータ作成
sales_base_df = pd.DataFrame({'date': pd.date_range(start='2022-01-01', periods=1200, freq='D')})
sales_base_df['date_idx'] = sales_base_df.index + 1
sales_base_df['weekday_num'] = sales_base_df['date'].dt.weekday
sales_base_df['weekday'] = sales_base_df['date'].dt.day_name()
sales_base_df['sales_base'] = sales_base_df.apply(lambda x: calc_sample_base(x)['base_sum'], axis=1)

# ポアソン分布に従う販売数データを作成
sales_sample_df = sales_base_df.copy()
sales_sample_df['sales'] = sales_sample_df['sales_base'].apply(lambda x: np.random.poisson(lam=x))
sales_sample_df = sales_sample_df[['date_idx', 'date', 'weekday', 'sales']]

作成したテーブルの中身は次のような内容になっています。

時系列として可視化すると次の通りです。

最初の40日分を拡大すると次の通りで、週末に売上が下がっているのが分かります。

シミュレーション用データの作成

次に、粗利が最適化できているか確かめるためのシミュレーション用データ(検証用データ)を作成します。今回は最後の100日分について、1000回シミュレーションを行います。

# シミュレーション回数
n_sim_iters = 1000

# 検証期間(シミュレーションで使用する日数)
n_valid_dates = 100

# 検証用データ
valid_df = sales_base_df.tail(n_valid_dates).reset_index(drop=True)

sim_sales_cols = []
for ii in range(n_sim_iters):
    sim_sales_col = f'sim_sales_{ii+1}'
    sim_sales_cols += [sim_sales_col]
    valid_df[sim_sales_col] = valid_df['sales_base'].apply(
        lambda x: np.random.poisson(lam=x)
    )

valid_df = valid_df[['date_idx', 'date', 'weekday'] + sim_sales_cols]

1000回分のシミュレーション用データを時系列として可視化すると次の通りです。

需要予測

次に需要予測を行います。今回は、需要予測をMeta社が開発した時系列予測モデルであるProphetを用いて行います。次のようにデータを分割し、1000日分で学習、合計200日分を予測します。

# trainとtestに分割
n_train = 1000
train_df = sales_sample_df[:n_train].reset_index(drop=True)
test_df = sales_sample_df[n_train:].reset_index(drop=True)
n_future = len(test_df) # 計200日

Prophetの学習と推論を行います。時系列予測(需要予測)自体は今回本質的な部分ではないので詳細は割愛します。

# prophet用trainデータの作成
train_prophet_df = train_df[['date', 'sales']].rename(
    columns={'date': 'ds', 'sales': 'y'}
)

prophet_model = Prophet()

# カスタム周期性 (120日) を追加
prophet_model.add_seasonality(name='120days_cycle', period=120, fourier_order=5)

# 学習を行う
prophet_model.fit(train_prophet_df)

# 推論を行う
future = prophet_model.make_future_dataframe(
    periods=n_future,
    freq='D',
    include_history=False
)
forecast_df = prophet_model.predict(future)

# 予測結果を格納
pred_df = test_df.copy()
pred_df['pred_sales'] = forecast_df['yhat']

予測結果を格納したテーブルpred_dfは次のような内容になっています。

なお、200日分のtestデータの予測に対して精度(MAPE)を計算すると約10%(0.0996)となりました。

prophet_test_mape = mape_val(
    pred_df['sales'],
    pred_df['pred_sales']
)
print('MAPE:', prophet_test_mape)
# MAPE: 0.09962763975531277

予測結果を時系列として可視化すると次の通りです。最後の100日分はシミュレーション用データに差し替えています。

拡大すると次の通りです。

需要推定量の誤差分布を取得

新聞売り子問題では需要の誤差分布を使用するため、1001日目から1100日目までの100日間分データを使用して、「 $t$ 日目の需要量が需要推定量からどの程度ぶれるかを表す量」 $\alpha_{t} = \frac{D_{t} - \hat{D_{t}}}{\hat{D_{t}}}$ の分布(密度関数)を取得します。

まずは、次の関数で各日付( $t$ 日)の $\alpha_{t}$ を計算します。

def calc_demand_diff_ratio_from_pred_demand(demand=None, pred_demand=None):
    """
    calc demand_diff_ratio_from_pred_demand
    
    Args:
        demand: float
            demand
        
        pred_demand: float
            demand
    
    Returns:
        alpha: float
            demand_diff_ratio_from_pred_demand
    """
    alpha = (demand - pred_demand) / pred_demand
    return alpha

# 誤差分布取得には最初の100日分を使用
calc_dist_df = pred_df.head(100)

# alphaを計算する
calc_dist_df['alpha'] = calc_dist_df.apply(
    lambda row: calc_demand_diff_ratio_from_pred_demand(row['sales'], row['pred_sales']),
    axis=1
)

次に、 $\alpha_{t}$ の分布(確率密度関数)をKDE(カーネル密度推定法)を用いて取得します。今回はscipy.statsのgaussian_kdeを用います。

※通常の新聞売り子問題ではこの誤差が正規分布に従うと仮定することも多いですが、実際のデータは正規分布に従うとも限らないため、今回は分布を仮定せずノンパラメトリックに行います。

# get alpha probability density function
sp_kde_model = scipy.stats.gaussian_kde(calc_dist_df['alpha'])
alpha_density_func = sp_kde_model.pdf
alpha_cumulative_dist_func = (lambda x: sp_kde_model.integrate_box_1d(-np.inf, x))

今回KDEによって取得した分布を可視化すると次の通りです。

新聞売り子問題を応用した在庫の最適化

取得した分布から新聞売り子問題を用いて「粗利の期待値」を最大化するような最適な「需要推定量に対する在庫数」を算出します。

ここでは「 $t$ 日の需要推定量に対して在庫をどの程度多く用意しておくべきか」を表す量 $x_{t} = \frac{q_{t} - \hat{D_{t}}}{\hat{D_{t}}}$ の最適解 $x_{t}^*$ を算出することにします。

$x_{t}^*$ は解析的に次の式で求められるのでした。

$\displaystyle{ \begin{aligned} x_{t}^* &= F ^ {-1} \left( \frac{p - c}{p} \right) \end{aligned} }$

これを実際にプログラムで実装すると次の通りです。

def calc_best_excess(
    x_min=-1,
    x_max=2,
    x_step=1e-3,
    alpha_cumulative_dist_func=None,
    alpha_density_func=None,
    price=None,
    cost=None,
    epsilon=1e-2
):
    """
    calc best excess stocking quantity.
    
    Args:
        x_min: float, optional(default=-1)
            min of x
        
        x_max: float, optional(default=2)
            max of x
        
        x_step: float, optional(default=1e-3)
            step of x
        
        alpha_cumulative_dist_func: function
            alpha cumulative distribution function.
            alpha is demand diff ratio from pred demand.
        
        alpha_density_func: function
            alpha probablity density function.
            alpha is demand diff ratio from pred demand.
        
        price: float or array-like of float
            price.
        
        cost: float or array-like of float
            cost.
        
        epsilon: float, optional(default=1e-2)
    
    Returns:
        best_x: float
            best_x
    """
    try:
        x_arr = np.arange(x_min, x_max + x_step, x_step)
        profit_ratio = (price - cost) / price
        
        result_list = []
        best_x = x_min
        
        for x in x_arr:    
            if alpha_cumulative_dist_func is None:
                area_val = calc_integrate(alpha_density_func, -np.inf, x, mode='quad')
            else:
                area_val = alpha_cumulative_dist_func(x)
            
            # 絶対残差
            abs_res = np.abs(profit_ratio - area_val)
            result_list += [[x, abs_res]]
            if abs_res <= epsilon:
                best_x = x
                break
        result_arr = np.array(result_list)
        best_x = result_arr[result_arr.argmin(axis=0)[1]][0]
        return best_x
    except Exception as e:
        print('error:')
        print(e)

# 需要推定量に対してどの程度多く在庫を用意するのが最適か算出する
best_x = calc_best_excess(
    x_min=-0.5,
    x_max=0.5,
    x_step=1e-3,
    alpha_cumulative_dist_func=alpha_cumulative_dist_func,
    price=price, # 100
    cost=cost, # 30
    epsilon=1e-3
)
print('best_x:', best_x)
# best_x: 0.05300000000000049

この結果からここでは特に「 $t$ 日の需要推定量に対して在庫をどの程度多く用意しておくべきか」の最適解 $x_{t}^*$ の値は約5.3%となります。

つまり、需要推定量の1.053倍の在庫を用意しておくのが最適であるという結果になります。

最後に、ここで出した最適解が正しいかどうかをあらかじめ用意しておいたシミュレーション用データを用いて検証します。

シミュレーション

今回は、「需要推定量に対する在庫の比率」を様々に変えて(0.7から1.3まで0.025刻み)、その時の「検証期間(100日間)の粗利合計(gp_sum)」の期待値(全シミュレーションの平均)を計算します。

def calc_gp_sum(rate=1, price=100, cost=30):
    """
    ある「需要推定量に対する在庫の比率」を指定したときの、
    各シミュレーションにおける「検証期間(100日間)の粗利合計(gp_sum)」を計算して
    listに格納して返す
    
    gross_profit = price * sales - cost * quantity
    
    Args:
        rate: float
            需要推定量に対する在庫の比率
        
        price: float
            販売価格
        
        cost: float
            原価
    """
    gp_sum_list = []
    temp_df = valid_df.copy()
    
    # 指定した比率に基づいて在庫を用意
    temp_df['quantity'] = temp_df['pred_sales'].apply(lambda x: round(x * rate))
    
    for ii in range(n_sim_iters):
        sim_sales_col = f'sim_sales_{ii+1}'
        temp_sales_col = f'true_sales_{ii+1}'
        temp_gp_col = f'gp_{ii+1}'
        
        # 実際の販売数は需要と在庫の内、小さい方の値
        temp_df[temp_sales_col] = temp_df[[sim_sales_col, 'quantity']].min(axis=1)
        temp_df[temp_gp_col] = price * temp_df[temp_sales_col] - cost * temp_df['quantity']
        
        # 検証期間(100日間)の粗利合計(gp_sum)を計算
        gp_sum = float(temp_df[temp_gp_col].sum())
        gp_sum_list += [gp_sum]
    return gp_sum_list

# シミュレーション用データに予測値をmerge
valid_df = pd.merge(
    valid_df,
    pred_df[['date_idx', 'pred_sales']],
    on='date_idx',
    how='left'
)

# 検証結果格納list
rate_gp_list = []

# 「需要推定量に対する在庫の比率」を様々に変えて
# 「検証期間(100日間)の粗利合計(gp_sum)」の期待値(全シミュレーションの平均)を計算
# 比率のパターンを0.025刻みで用意
for _rate in list(np.arange(0.7, 1.3, 0.025)):
    temp_gp_sum_list = calc_gp_sum(
        rate=_rate,
        price=price, # 100
        cost=cost # 30
    )
    
    exp_gp = np.mean(temp_gp_sum_list)
    rate_gp_list += [[_rate, exp_gp]]

# 検証結果をDataFrameに格納
rate_gp_df = pd.DataFrame(
    rate_gp_list,
    columns=['rate', 'expected_gross_profit']
).sort_values(
    by=['expected_gross_profit'],
    ascending=False # 降順にsort
).reset_index(drop=True)