<jupyter_start><jupyter_text>Load modules<jupyter_code>import sys, os
import gc
import time
import random
import scipy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# original modules
from annin_dofu import calc, matrix, parallel, stats, utils

# news vendor modelのimport
from annin_dofu.stocking_quantity_optimization import news_vendor_model as nvm<jupyter_output><empty_output><jupyter_text>Configuration  random seed<jupyter_code>np.random.seed(57)
random.seed(57)<jupyter_output><empty_output><jupyter_text>warningの非表示<jupyter_code># warningの非表示
import warnings
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>左寄せにするマジックコマンド<jupyter_code>%%html
<style>
    table{float:left}
    .MathJax{float: left;}
</style><jupyter_output><empty_output><jupyter_text>データフレームの表示設定<jupyter_code># データフレームの表示行数, 表示列数設定.
pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 200)

# カラム内の文字数設定. デフォルトは50.
pd.set_option('display.max_colwidth', 100)<jupyter_output><empty_output><jupyter_text>カレントディレクトリの設定<jupyter_code># 念の為カレントディレクトリをファイルの場所に変更しておく
os.chdir(os.getcwd())<jupyter_output><empty_output><jupyter_text>経過時間取得のための開始時間保存<jupyter_code>start_time = time.time()<jupyter_output><empty_output><jupyter_text>Functions  norm_dist_prob 正規分布$$\begin{aligned}\frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right)\end{aligned}$$<jupyter_code>stats.norm_dist_prob(1)
stats.sp_norm_dist_prob(1)
x_arr = np.linspace(-5, 5, 1000)
y_arr = stats.norm_dist_prob(x_arr)

plt.figure(facecolor='white')
plt.plot(x_arr, y_arr)
plt.title('normal distribution')
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>calc_expected_gross_profit 需要推定量 $\hat{D_{t}}$ が与えられた時の粗利(gross profit)の期待値は次式で与えられる.$$\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}$$  calc_best_expected_gross_profit - 最適解における粗利期待値$$\begin{aligned}E[\mathrm{profit}_{t} \mid x_{t} = 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}}\begin{Bmatrix}\begin{aligned}p \int_{u \leq x_{t}^*} u f(u) du - p x_{t}^* \left( \frac{p - c}{p} \right) + (p - 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 - c) 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 f(u) du + p - c\end{aligned}\end{Bmatrix}\end{aligned}$$ - 増減表粗利期待値の増減表は次の通り.| $x_{t}$ | $\cdots$ | $x_{t}^*$ | $\cdots$ ||:-:|:-:|:-:|:-:|| $\dfrac{\partial}{\partial x_{t}} E[\mathrm{profit}_{t}]$ | $+$ | $0$ | $-$ || $E[\mathrm{profit}_{t}]$ | $\nearrow$ | $$\hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}^*} u f(u) du + p - c \end{aligned} \end{Bmatrix}$$ | $\searrow$ |  make sample data  variables- 変数の値を設定    - 売価    - 原価<jupyter_code># 売価
price = 100

# 原価
cost = 30

# sample size
n_samples = 100

# 周期性
cycle = 7<jupyter_output><empty_output><jupyter_text>時間(time)<jupyter_code># 時間(time)
t_list = np.arange(1, n_samples+1, 1)<jupyter_output><empty_output><jupyter_text>需要量(demand), 在庫量(stocking quantity), 販売量(sales)<jupyter_code># 時系列サンプル
ts_samples_list = np.round(30 + 0.5 * t_list + 10 * np.sin(2 * np.pi * t_list / cycle))

# 需要量(demand)(今回求めたい分布)
# d_list = np.round(np.random.normal(100, 20, n_samples)) # 正規分布
# d_list = np.round(np.random.poisson(100, n_samples)) # ポアソン分布
d_list = ts_samples_list + np.round(np.random.poisson(30, n_samples)) # ポアソン分布

# 在庫量(stocking quantity)
# q_list = np.round(np.random.normal(110, 15, n_samples))
q_list = ts_samples_list + np.round(np.random.normal(35, 10, n_samples))

# 販売量(sales)
# s_list = np.round(np.min(np.vstack([d_list, q_list]), axis=0))
s_list = np.min([q_list, d_list], axis=0) # 高速<jupyter_output><empty_output><jupyter_text>粗利(gross profit)- $t$ 便の粗利(gross profit)を次式で定義する.$$\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}$$- nvm.calc_gross_profit()で計算する.<jupyter_code># 粗利(gross_profit)
gp_list = nvm.calc_gross_profit(
    sales=s_list,
    quantity=q_list,
    price=price,
    cost=cost
)<jupyter_output><empty_output><jupyter_text>pandas.DataFrameに格納<jupyter_code>df = pd.DataFrame({
    'time': t_list,
    'demand': d_list,
    'quantity': q_list,
    'sales': s_list,
    'price': price,
    'cost': cost,
    'profit': gp_list
}).astype(int)
df.head(10)
df.shape<jupyter_output><empty_output><jupyter_text>Preprocessing  pred_demand- 需要予測値を格納する.<jupyter_code># 時系列予測が上手くいった場合
df['pred_demand'] = ts_samples_list + 30

# 平均値を予測値とする場合
# df['pred_demand'] = df['demand'].mean()
df.head(3)<jupyter_output><empty_output><jupyter_text>alpha- $\alpha$を次式で定義する.$$\begin{aligned}&\alpha_{t} = \frac{D_{t} - \hat{D_{t}}}{\hat{D_{t}}} : \text{$t$便の需要推定量に対して需要量が需要推定量からどの程度ぶれるかを表す量. 確率変数.}\end{aligned}$$- nvm.calc_demand_diff_ratio_from_pred_demand()で計算する.<jupyter_code>df['alpha'] = df.apply(lambda row: nvm.calc_demand_diff_ratio_from_pred_demand(row['demand'], row['pred_demand']), axis=1)<jupyter_output><empty_output><jupyter_text>x- $x$を次式で定義する.$$\begin{aligned}&x_{t} = \frac{q_{t} - \hat{D_{t}}}{\hat{D_{t}}} \in \mathbb{R} : \text{$t$便の需要推定量に対する期初総在庫数と需要推定量の差分の割合. どの程度多く在庫を用意しておくか. 推定される販売あたり廃棄率.}\end{aligned}$$- nvm.calc_quantity_diff_ratio_from_pred_demand()で計算する.<jupyter_code>df['x'] = df.apply(lambda row: nvm.calc_quantity_diff_ratio_from_pred_demand(row['quantity'], row['pred_demand']), axis=1)<jupyter_output><empty_output><jupyter_text>check DataFrame  head and tail<jupyter_code>df.head()
df.tail()
df.shape<jupyter_output><empty_output><jupyter_text>describe<jupyter_code>df.describe()<jupyter_output><empty_output><jupyter_text>Stocking Quantity Optimization  get alpha probability density function- $\alpha$の確率密度関数と確率分布を取得する.- 確率分布の関数をalpha_cumulative_dist_funcという変数に格納する.<jupyter_code># get alpha probability density function
sp_kde_model = scipy.stats.gaussian_kde(df['alpha'])
alpha_density_func = sp_kde_model.pdf
alpha_cumulative_dist_func = (lambda x: sp_kde_model.integrate_box_1d(-np.inf, x))<jupyter_output><empty_output><jupyter_text>calc best_excess- $\alpha$の確率分布関数を元に, 最適な$x$を取得する.$F$ は単調増加より この式を満たす $x_{t}$ はただ1つ存在し,それを $x_{t}^*$ とする.($F ^ {-1}$ は逆像ではなく, 逆写像)$$\begin{aligned}F (x_{t}^*) &= \int_{-\infty}^{p} f(u) du \\&= \frac{p - c}{p}\end{aligned}$$$$\begin{aligned}x_{t}^* &= F ^ {-1} \left( \frac{p - c}{p} \right)\end{aligned}$$- nvm.calc_best_excess()で計算する.<jupyter_code>best_x = nvm.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,
    cost=cost,
    epsilon=1e-3
)
print('best_x:', best_x)<jupyter_output><empty_output><jupyter_text>calc best_gross_profit- 最適解における粗利期待値$$\begin{aligned}E[\mathrm{profit}_{t} \mid x_{t} = 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}}\begin{Bmatrix}\begin{aligned}p \int_{u \leq x_{t}^*} u f(u) du - p x_{t}^* \left( \frac{p - c}{p} \right) + (p - 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 - c) 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 f(u) du + p - c\end{aligned}\end{Bmatrix}\end{aligned}$$- 増減表粗利期待値の増減表は次の通り.| $x_{t}$ | $\cdots$ | $x_{t}^*$ | $\cdots$ ||:-:|:-:|:-:|:-:|| $\dfrac{\partial}{\partial x_{t}} E[\mathrm{profit}_{t}]$ | $+$ | $0$ | $-$ || $E[\mathrm{profit}_{t}]$ | $\nearrow$ | $$\hat{D_{t}} \begin{Bmatrix} \begin{aligned} p \int_{u \leq x_{t}^*} u f(u) du + p - c \end{aligned} \end{Bmatrix}$$ | $\searrow$ |- nvm.calc_best_expected_gross_profit()で計算する.<jupyter_code># 平均を予測値とする
pred_demand = df['pred_demand'].mean()

print('pred demand:', pred_demand)
# 最適な粗利期待値
# nvm.calc_best_expected_gross_profit()で直接算出する.
best_egp1 = nvm.calc_best_expected_gross_profit(
    best_x=best_x,
    alpha_density_func=alpha_density_func,
    pred_demand=pred_demand,
    price=price,
    cost=cost
)
print('best expected_gross_profit:', best_egp1)
# 最適な粗利期待値
# nvm.calc_expected_gross_profit()に最適なxを代入して算出する.
best_egp2 = nvm.calc_expected_gross_profit(
    x=best_x,
    alpha_density_func=alpha_density_func,
    pred_demand=pred_demand,
    price=price,
    cost=cost
)
print('best expected_gross_profit:', best_egp2)
# 2通りの算出で結果がほぼ同じになることが確認できる.
print('best expected_gross_profit')
print('直接算出:', best_egp1)
print('最適値を代入して算出:', best_egp2)<jupyter_output><empty_output><jupyter_text>Visualization  time series  需要量(demand), 在庫量(stocking quantity), 販売量(sales)<jupyter_code>plt.figure(
    figsize=(20, 5),
    facecolor='white',
    dpi=150
)

plt.title('sample data\ndemand, quantity, sales', fontsize=20)
plt.xlabel('time', fontsize=15)
plt.ylabel('value', fontsize=15)

# demand
plt.plot(
    df['time'],
    df['demand'],
    label='demand',
    lw=3
)

# quantity
plt.plot(
    df['time'],
    df['quantity'],
    label='quantity',
    lw=3
)

# sales
plt.plot(
    df['time'],
    df['sales'],
    label='sales',
    lw=2
)

plt.ylim([0, 160])
plt.grid()
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>alpha, x<jupyter_code>plt.figure(
    figsize=(20, 5),
    facecolor='white',
    dpi=150
)

plt.title('sample data\nalpha, x', fontsize=20)
plt.xlabel('time', fontsize=15)
plt.ylabel('value', fontsize=15)

# alpha
plt.plot(
    df['time'],
    df['alpha'],
    label='alpha',
    lw=3
)

# x
plt.plot(
    df['time'],
    df['x'],
    label='x',
    lw=3
)

plt.ylim([-0.5, 0.5])
plt.grid()
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>粗利(gross profit)<jupyter_code>fig, ax1 = plt.subplots(
    figsize=(40, 10),
    facecolor='white',
    dpi=150
)
ax2 = ax1.twinx()

# demand
g1 = ax1.plot(
    df['time'],
    df['demand'],
    label='demand',
    lw=2,
    alpha=0.5
)

# quantity
g2 = ax1.plot(
    df['time'],
    df['quantity'],
    label='quantity',
    lw=2,
    alpha=0.5
)

# sales
g3 = ax1.plot(
    df['time'],
    df['sales'],
    label='sales',
    lw=2,
    alpha=0.5
)

# gross_profit
g4 = ax2.plot(
    df['time'],
    df['profit'],
    label='gross_profit',
    lw=4,
    color='red'
)

_title = r'gross_profit'
ax1.set_title(_title, fontsize=20)
ax1.set_xlabel(r'time', fontsize=15)
ax1.set_ylabel(r'sales', fontsize=15)
ax2.set_ylabel(r'gross_profit', fontsize=15)

g_list = g1 + g2 + g3 + g4
plt.legend(
    g_list,
    [g.get_label() for g in g_list],
    loc=0,
    fontsize=15
)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>histogram and distribution  需要量(demand)<jupyter_code>var = 'demand'
plt.figure(figsize=(6, 5), facecolor='white')
sns.distplot(
    df[var],
    label=var
)
plt.title(var)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>在庫量(stocking quantity)<jupyter_code>var = 'quantity'
plt.figure(figsize=(6, 5), facecolor='white')
sns.distplot(
    df[var],
    label=var
)
plt.title(var)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>販売量(sales)<jupyter_code>var = 'sales'
plt.figure(figsize=(6, 5), facecolor='white')
sns.distplot(
    df[var],
    label=var
)
plt.title(var)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>需要量(demand), 在庫量(stocking quantity), 販売量(sales)<jupyter_code>vars_list = [
    'demand',
    'quantity',
    'sales'
]
plt.figure(
    figsize=(10, 5),
    facecolor='white',
    dpi=150
)

for var in vars_list:
    sns.distplot(
        df[var],
        label=var
    )

plt.title('hist and density\ndemand, quantity, sales')
plt.grid()
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>alpha<jupyter_code>grid_x = np.linspace(-0.2, 0.2, 1000)
m_inf_arr = np.array([-np.inf for x in grid_x])

prob_y = alpha_density_func(grid_x)
integral_y = np.array([alpha_cumulative_dist_func(x) for x in grid_x])
# plt.figure(figsize=(10, 5), facecolor='white')
fig, ax1 = plt.subplots(
    figsize=(10, 5),
    facecolor='white',
    dpi=150
)
ax2 = ax1.twinx()

g1 = ax1.plot(
    grid_x,
    prob_y,
    color='blue',
    label='Probablity density function'
)

ax1.fill_between(
    x=grid_x,
    y1=prob_y,
    y2=0,
    color='orange',
    alpha=0.5
)

ax1.hist(
    df['alpha'],
    bins=20,
    density=True,
    alpha=0.3
)

g2 = ax2.plot(
    grid_x,
    integral_y,
    color='red',
    label='Cumulative distribution function'
)

_title = r'$\alpha$ - distribution'
ax1.set_title(_title, fontsize=15)
ax1.set_xlabel(r'$\alpha = \dfrac{D - \hat{D}}{\hat{D}}$', fontsize=15)

ax1.set_xlim([-0.21, 0.21])
ax1.set_ylim([-0.1, 7])

g_list = g1 + g2
plt.legend(
    g_list,
    [g.get_label() for g in g_list],
    loc=0
)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>x<jupyter_code>var = 'x'
plt.figure(figsize=(6, 5), facecolor='white')
sns.distplot(
    df[var],
    label=var
)
plt.title(var)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>expected gross profit- xを様々に変化させた時の粗利期待値を可視化し, 最適解best_xが正しいことを確認する.<jupyter_code># 平均を予測値とする
pred_demand = df['pred_demand'].mean()
pred_demand
# 様々なxの値
x_arr = np.arange(-0.5, 0.51, 1e-3)
x_arr[:5], x_arr[-5:]
# 粗利期待値を計算
egp_arr = np.array([
    nvm.calc_expected_gross_profit(
        x=x,
        alpha_density_func=alpha_density_func,
        pred_demand=pred_demand,
        price=price,
        cost=cost
    ) for x in x_arr])
egp_arr[:3]
x_egp_arr = np.vstack([x_arr, egp_arr]).T
x_egp_df = pd.DataFrame(x_egp_arr, columns=['x', 'expected_gross_profit'])
x_egp_df.head()
x_arr.shape, egp_arr.shape, x_egp_arr.shape
# 粗利期待値が最大となるxの値の前後
ii = x_egp_df['expected_gross_profit'].idxmax(axis=0)

tmp_best_x = x_egp_df.iloc[ii]['x']
tmp_best_egp = x_egp_df.iloc[ii]['expected_gross_profit']
print('best_x:', tmp_best_x)
print('best_egp:', tmp_best_egp)

x_egp_df.iloc[ii-2:ii+3]
# xを様々に変化させた時の粗利期待値を可視化し, 最適解best_xが正しいことを確認する.
plt.figure(
    figsize=(12, 8),
    facecolor='white',
    dpi=150
)

# x and egp
plt.plot(
    x_arr,
    egp_arr,
    color='blue'
)

# best_x
plt.vlines(
    tmp_best_x,
    -10000, 10000,
    color='red',
    ls='--',
    label='best_x'
)

# best_egp
plt.hlines(
    tmp_best_egp,
    -10, 10,
    color='red',
    label='best_expected_gross_profit'
)

_title = 'expected gross profit\n'
_title += 'pred_demand: {}'.format(pred_demand)

plt.title(_title, fontsize=15)
plt.xlabel('x', fontsize=15)
plt.ylabel('expected gross profit', fontsize=15)

plt.xticks(fontsize=10)
plt.yticks(fontsize=10)

plt.xlim([-0.55, 0.55])
plt.ylim([0, egp_arr.max() + 1000])

plt.grid()
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Check  get elapsed time<jupyter_code>utils.get_elapsed_time(start_time)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20065363507430475, 'fraction_numerical': 0.018992415366590614, 'mean_word_length': 4.903895158354569, 'pattern_counts': {'":': 0, '<': 131, '<?xml version=': 0, '>': 131, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '559482', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 7960, 'n_tokens_neox': 7595, 'n_words': 1072, 'path': 'examples/annin_dofu_stocking_quantity_optimization_sample.ipynb', 'repo_name': 'laplaciannin102/annin_dofu', 'size': '16217'}