Zenn
OPENLOGI Tech BlogPublicationへの投稿
📦

ガンマ分布に基づく発注計算

akabe
2024/12/17に公開
math
統計
物流
発注
idea

OPENLOGI Advent Calendar 2024 の記事です。

物流系テック企業のオープンロジで機械学習エンジニア・データエンジニアとして働いていてる阿部です。最近、発注や在庫関連の仕事をしていて、自由研究として新しい発注計算の方法を思いついたので、紹介します。

TL;DR

商品の発注量・発注点を計算する際、欠品する確率をコントロールするために、安全在庫という安全マージンを使うことが多いです。安全在庫は「商品の1日の販売個数が正規分布する」ことを仮定して計算することが一般的ですが、安全在庫が小さく見積もられることが多く、欠品のリスクが上がります。より正確に欠品確率を調整できる手段として、ガンマ分布を採用して発注量・発注点を計算する方法を考えました。以下の方法で計算できます（四則演算だけ計算できるロジックにしています）。

ガンマ分布を使うことで欠品許容率に近い欠品率に収めることができます。一方で、発注量がやや大きく見積もられる傾向にあり、保管料が上がるデメリットがあるため、過剰在庫の抑制が今後の課題です。

まえおき

発注計算

物流において、店舗や倉庫に保管されている商品の在庫量を考えることは非常に重要です。在庫が少なすぎると、欠品による機会損失が発生し、逆に在庫が多すぎると、余剰在庫の保管料が経営を圧迫することになります。適切な在庫量を保つためには売れ行き（＝販売数量の傾向）に応じて発注量や発注のタイミングをコントロールする必要があり、いくつかの方法が存在します。ここでは、定期発注と定量発注を紹介します。

この 2 つの式に共通しているのは、適当な時間間隔 T（= リードタイム or 発注間隔 + リードタイム）について

T × 1日の平均販売個数 + 安全在庫

という計算を行っていることです。この式を理解するために、次の節で安全在庫について説明します。

安全在庫

結論から言うと、T × 1日の平均販売個数 + 安全在庫 という式は「T 日間で在庫切れを起こす確率を一定に抑える在庫量」を意味しています。例えば、T = 1 ヶ月とすると、1 ヶ月間に販売できる商品の平均個数は 30日 × 1日の平均販売個数 ですね。仮に月初の時点で 30日 × 1日の平均販売個数 の数量の商品が在庫として確保されている場合、月末まで在庫切れ（欠品）を起こす確率は 50% になります。欠品確率 50% はとても大きい数字で、2 ヶ月に 1 回は欠品を起こし、機会損失が発生することになります。欠品確率を 5% や 1% などの一定の値に制御することが一般的ですが、そのためには 1 ヶ月の平均的な販売個数 (30日 × 1日の平均販売個数) に加えて、少し余分に在庫を確保する必要があります。この余分な在庫量が 安全在庫 に該当します。

以前、「安全在庫の考え方を理解する」という記事で一般的によく使われる安全在庫の計算式を紹介しました。

安全在庫 = 安全係数 × 1日の販売個数の標準偏差 × sqrt(T)

安全係数は欠品許容率から決まる係数で、例えば欠品許容率 5% では安全係数に 1.65 を使います。

また、上記の計算式を導出する上で

  • 日々の販売数は独立な正規分布に従う
  • T は常に一定

という仮定を置いていることを述べました。

販売数は正規分布に従うか？

安全在庫について解説している多くの文献では、特に根拠なく正規分布を仮定していることが多いです。しかし、私はずっと正規分布するのか疑問に思っていました。まず、正規分布する仮定を疑って、別な確率分布の可能性がないか探ってみます。

既存研究

個人的な予想だけで話をしても説得力に欠けるので、既存研究を調べてみます。日々の販売量（＝一定期間での商品の販売個数）の確率分布について調査した文献を調べてみると、意外にもどのような確率分布を使うか意見が分かれているようです。しかし、いずれの文献でも左に偏った分布を採用しています。

実際のデータで確認してみる

オープンロジの商品出荷実績データを使って、実際にどの確率分布が近しいか確認してみます。結論から言うと、ガンマ分布が近しいと思われます。今回は以下の手順で、各確率分布と実データの KL (Kullback–Leibler) 情報量 を計算して、どの確率分布が実測値によく当てはまるか確認しました。

  1. 商品の 1 日ごとの出荷数のデータを集計します。
    • オープンロジの中でも一定の販売実績がある商品約 2500 種類について 6 ヶ月分の販売データを利用して集計しました。
  2. 正規分布・ポアソン分布・ガンマ分布について商品ごとにパラメータを最尤推定します。
  3. 実際の出荷数のヒストグラムと正規分布・ポアソン分布・ガンマ分布の KL 情報量を計算します。
    • KL 情報量：D_\mathrm{KL}(P \| Q) = \sum_i P(i) \log(P(i) / Q(i))
    • 直感的には 2 種類の分布 P, Q の類似度みたいなもので、数字が小さいほどよく似た分布形状となります。
    • P をヒストグラムから計算した実際の分布、Q はパラメータを推定した各確率分布の理論値を当てはめています。

※ ワイブル分布はパラメータ推定の方法がよくわかっていないので、除外しています。詳しい人はやり方を教えてください。。。

確率分布ごとに、KL 情報量の分布と中央値をヴァイオリンプロットしたものが以下の図です。

この図をみる限り、ガンマ分布の KL 情報量が小さく、実際のデータに当てはまりが良いことがわかります。

実際のヒストグラムと最尤推定された確率分布をプロット [1] してみると、ガンマ分布が一番近い形をしているように見えます。ガンマ分布よりもやや減衰が緩やかな曲線に見えますが、正規分布やポアソン分布よりは近い形状になっています。

「一定期間での商品の販売個数」がガンマ分布に従う理由やメカニズムについて、私の中では仮説を持てていません。ここで述べた根拠だけでは「商品の販売個数はガンマ分布に従う」と断言することはできませんが、正規分布よりもマシな確率分布を使うという雑な考え方のもと、今回はガンマ分布を仮定したいと思います。商品の販売数はトレンドの変化やセールなど様々な要因に影響を受けるため、おそらくどの確率分布を持ってきてもある程度の差異は生じると思います。

ガンマ分布による発注計算

ガンマ分布とは？

今更ですが、ガンマ分布 (Gamma distribution) とは指数分布を一般化した確率分布で

  • 人の体重
  • ウイルスの潜伏期間
  • 電子部品の寿命
  • トラフィックの待ち時間
  • 所得

などがガンマ分布に従うと言われています。形状パラメータ k > 0、尺度パラメータ \theta > 0 を用いて、確率分布関数は

\mathit{Ga}(x | k, \theta) = \frac{1}{\Gamma(k) \theta^k} x^{k-1} e^{-x / \theta} \qquad (x > 0)

ガンマ分布の確率密度関数
ガンマ分布の確率密度関数

と表されます。ガンマ関数 \Gamma(k) は階乗を連続量に拡張したもので、自然数について k! = \Gamma(k+1) が成立します。

\Gamma(z) = \int^\infty_0 t^{z-1} e^{-t} \mathrm{d}t

発注計算式の導出過程

日付を t として、1 日単位での販売数量を x_t とおきます。これまでの議論より、x_t は独立なガンマ分布に従うものとします。また、パラメータ k, \theta は t に依らず一定であることを仮定します。

x_t \sim \mathit{Ga}(x | k, \theta)

ここでは、例として発注点方式を採用し、発注時に安全在庫を計算することを考えます。発注してから販売可能な在庫に計上されるまで T 日間かかるとすると、T 日間の合計販売数量 x_{1\dots T} の分布は、ガンマ分布の再生性により、

x_{1\dots T} = x_1 + \cdots + x_{T} \sim \mathit{Ga}(x | T k, \theta)

となります。

発注数量を z 個とすると、z < x_{1\dots T} である場合に販売数が発注数量を超過して欠品が発生します。欠品が発生する確率は

\begin{align*} P(z < x_{1\dots T}) & = 1 - P(z > x_{1\dots T}) \\ & = 1 - \int^z_0 \mathit{Ga}(x | T k, \theta) \mathrm{d}x \\ & = 1 - \hat{\gamma}(T k, z / \theta) \end{align*}

となります。\hat\gamma は正則化第1種不完全ガンマ関数であり、ガンマ分布の累積分布関数です。

\hat{\gamma}(a, x) = \left( \int^x_0 t^{a-1} e^{-t} \mathrm{d} t \right) / \Gamma(a)

正則化第1種不完全ガンマ関数
正則化第1種不完全ガンマ関数

事前にビジネス上許容される欠品率 p（欠品許容率）を欠品の機会損失や保管料を元に決めているものとします。一般的には 5% がよく使われますが、任意に決めることができます。理論上の欠品率 P(z < x_{1\dots T}) と欠品許容率が一致しなければいけないので、

p = P(z < x_{1\dots T}) = 1 - \hat{\gamma}(T k, z / \theta)

となるような z を求めれば良いことがわかります。\hat\gamma の第 2 引数に対する逆関数を \hat\gamma^{-1} とすると、

z = \hat{\gamma}^{-1}(T k, 1 - p) \cdot \theta

で発注量・発注点を計算することができます。この式は冒頭で紹介した T × 1日の平均販売個数 + 安全在庫 の代わりに使うことを想定しています。

ただし、不完全ガンマ関数の逆関数は解析的に計算できないので、数値的な近似解を求めるしかありません。この方法についてはガンマ分布に関連した数値計算のあれこれで説明しているので、ご覧ください。また、ありがちな Tk, p の組み合わせで \hat{\gamma}^{-1}(T k, 1 - p) の値を計算した表を作ったので、プログラムを組むことが面倒な人はご利用ください（Appendix: 安全係数表）。

ここまでの計算を踏まえると、最初に紹介した以下のロジックで発注に関する計算ができることがわかります。

  1. 欠品許容率 p を決める。
  2. 1日単位で商品の販売個数の実績値を集計する。
  3. 日次販売個数からガンマ分布のパラメータ k, \theta を推定する。
  4. 発注方式に合わせて計算する。
    • 定期発注：\mathtt{発注量} = \hat{\gamma}^{-1}(T k, 1 - p) \cdot \theta - \mathtt{現在の在庫量} - \mathtt{到着残数}
      • ただし、T = \mathtt{リードタイム} + \mathtt{発注間隔}
    • 定量発注：\mathtt{発注点} = \hat{\gamma}^{-1}(T k, 1 - p) \cdot \theta
      • ただし、T = \mathtt{リードタイム}

欠品率のシミュレーション

シミュレーション方法

十分な販売実績がある商品およそ 2500 種類について、6ヶ月間の販売実績を用いて以下のシミュレーションを行います。

  1. 販売実績から正規分布 (\mu, \sigma)・ガンマ分布 (k, \theta) のパラメータをそれぞれ最尤推定する。
  2. パラメータを元に、発注量を計算する。
    • 欠品許容率 p = 0.05 (5%)
    • 正規分布ベースの発注量: z_\mathrm{norm} = T \mu + \alpha \sigma \sqrt{T} (安全係数 \alpha = 1.65)
    • ガンマ分布ベースの発注量: z_\mathrm{gamma} = \hat{\gamma}^{-1}(T k, 1-p) \cdot \theta
  3. 連続した T 日間の期間の合計販売数量を移動平均のロジックで計算し、発注量を超過している日数を調べて、欠品率を計算する。
  4. シミュレーションの欠品率と予め決めた欠品率 5% を比較する。

ガンマ分布を使う方が欠品率と欠品許容率が近づく

T = 1, 7, 14, 30, 60 について、シミュレーションの欠品率をプロットすると以下のようになりました。実線は欠品率の平均値で、薄く色をつけた部分は 10%tile〜90%tile を可視化しています。赤い線は欠品率 5% のラインで、一致していれば狙った欠品率になっていることを示しています。

シミュレーションの欠品率の比較
シミュレーションの欠品率の比較

また、欠品率の平均値を表にすると、以下の通りです。

T 平均欠品率（正規分布） 平均欠品率（ガンマ分布）
1 3.7 % 4.5 %
7 7.7 % 3.3 %
14 10.6 % 4.5 %
30 14.5 % 6.7 %
60 16.2 % 7.6 %

いずれの T でもガンマ分布を使う方が 5% (0.05) に近い数値が得られていることがわかります。正規分布を使うと、欠品許容率 5% を設定したにも関わらず、実際の欠品率が 16% 程度になっており、欠品リスクが大きいことがわかります。このことから、正規分布を元にした発注量（安全在庫）は過剰に小さく見積もられていると言えます。

ガンマ分布を使うと在庫量が過剰になるケースがある

日次販売個数のヒストグラムがファットテールを持つ曲線になる場合に、正規分布では欠品率が 50% を超えます。このようなケースにおいて、ガンマ分布では欠品率が 0% になっていました。欠品率が 0% ということは発注量が多すぎて、過剰な在庫を抱えてしまったことを意味します。前半でも述べましたが、ヒストグラム（実際の分布）とガンマ分布を見比べると、ガンマ分布よりもやや減衰が緩やかな曲線に見えます。このため、ガンマ分布を使う場合、商品の販売個数が多い日が発生する確率を実際よりも多めに見積もっている可能性があります。


正規分布で欠品率が高くなる商品の販売個数ヒストグラム

今後の課題

ガンマ分布を使うことで、事前に決めた欠品許容率に近い欠品率に制御することができ、なかなか性能よく発注量を計算できるようになりました。しかし、一部で発注量が過剰になり、欠品率が 0% になってしまう事例もありました。欠品率は低ければ良いというものではありません。低すぎると過剰在庫により保管料が増えてしまうため、欠品許容率と実際の欠品率はほぼ同じになることが望ましいです。今後は過剰在庫を抑制する方法を考えたいところです。

Appendix

安全係数表

正則化第1種不完全ガンマ関数の逆関数 \hat{\gamma}^{-1}(T k, 1 - p) を数値計算で求めた表です。列が p（欠品許容率）、行が T k で別れています。

T k ＼ p 0.1 0.05 0.01 0.005 0.001
0.01 1.50364e-05 0.00336279 0.265055 0.556058 1.50928
0.02 0.00294966 0.0459111 0.558911 0.93414 2.0057
0.03 0.0174561 0.116137 0.774528 1.18936 2.31793
0.04 0.0433841 0.1919 0.945289 1.3847 2.5497
0.05 0.0763176 0.265934 1.08758 1.54456 2.73637
0.06 0.112877 0.336272 1.21039 1.68113 2.89382
0.07 0.151042 0.40264 1.31906 1.80066 3.03056
0.08 0.189687 0.465207 1.41679 1.90772 3.15271
0.09 0.228174 0.524353 1.50616 2.00505 3.26269
0.1 0.266155 0.580436 1.58847 2.09449 3.36294
0.2 0.604902 1.03052 2.20233 2.75461 4.10119
0.3 0.884808 1.37235 2.63947 3.22049 4.6198
0.4 1.12984 1.66195 3.00008 3.60337 5.04274
0.5 1.35277 1.92074 3.31751 3.93964 5.41371
0.6 1.5605 2.15902 3.60658 4.24564 5.75087
0.7 1.75713 2.38263 3.87573 4.52978 6.06384
0.8 1.94526 2.59514 4.12993 4.79779 6.35865
0.9 2.12667 2.79896 4.3723 5.05328 6.63886
1 2.30259 2.99572 4.6051 5.29831 6.90686
2 3.88973 4.74388 6.63837 7.4301 9.23295
3 5.32232 6.29577 8.40592 9.27355 11.2295
4 6.6808 7.75366 10.0451 10.9772 13.0629
5 7.99361 9.15352 11.6046 12.5942 14.7952
6 9.27467 10.513 13.1086 14.1498 16.4546
7 10.5321 11.8424 14.5706 15.6599 18.0619
8 11.7709 13.1481 15.9999 17.1335 19.626
9 12.9947 14.4346 17.4027 18.5781 21.1577
10 14.206 15.7052 18.7832 19.9987 22.657
11 15.4066 16.9622 20.1446 21.3978 24.1341
12 16.5981 18.2075 21.4899 22.7792 25.5908
13 17.7816 19.4426 22.8209 24.1449 27.0261
14 18.9579 20.6686 24.1393 25.4967 28.4464
15 20.128 21.8865 25.4461 26.8356 29.8517
16 21.2924 23.0971 26.7429 28.1642 31.2424
17 22.4516 24.3012 28.0306 29.4818 32.623
18 23.6061 25.4992 29.3095 30.7906 33.9908
19 24.7563 26.6917 30.5809 32.0904 35.3512
20 25.9025 27.8792 31.8454 33.3827 36.702
21 27.0451 29.062 33.1032 34.6677 38.0411
22 28.1843 30.2404 34.355 35.9465 39.3747
23 29.3203 31.4148 35.6006 37.2186 40.701
24 30.4533 32.5854 36.8412 38.4845 42.0186
25 31.5836 33.7524 38.077 39.7447 43.3301
26 32.7112 34.9161 39.3079 41.0001 44.634
27 33.8364 36.0766 40.5346 42.2506 45.936
28 34.9592 37.2342 41.7568 43.4973 47.2324
29 36.0799 38.3889 42.9753 44.7389 48.5198
30 37.1985 39.541 44.1898 45.9756 49.8041
31 38.3151 40.6905 45.4009 47.2088 51.0844
32 39.4298 41.8376 46.6085 48.4388 52.3597
33 40.5427 42.9825 47.8127 49.6652 53.629
34 41.6539 44.1251 49.0141 50.888 54.8965
35 42.7635 45.2656 50.2125 52.1074 56.1563
36 43.8715 46.4041 51.4084 53.3242 57.4181
37 44.978 47.5407 52.6011 54.5369 58.6733
38 46.0831 48.6754 53.7913 55.7481 59.9238
39 47.1868 49.8085 54.979 56.9553 61.1716
40 48.2891 50.9397 56.1643 58.1605 62.4189
41 49.3902 52.0694 57.3473 59.3632 63.6597
42 50.49 53.1974 58.5281 60.5627 64.9018
43 51.5886 54.3239 59.707 61.7608 66.1361
44 52.6861 55.449 60.8835 62.9562 67.3708
45 53.7825 56.5726 62.0581 64.1491 68.6025
46 54.8778 57.6949 63.2308 65.3407 69.8308
47 55.9721 58.8158 64.4016 66.5297 71.0613
48 57.0653 59.9355 65.5703 67.7164 72.2812
49 58.1577 61.0539 66.7378 68.9013 73.5027
50 59.249 62.1711 67.9033 70.0846 74.7253
脚注

  1. 横軸の数字は都合上あえて省略しています。 ↩︎

OPENLOGI Tech Blog
OPENLOGI Tech BlogPublication

私たちは身近な社会課題である物流の未来を変えるべく、「物流版AWS」をコンセプトとしたプラットフォームの開発、データの活用によって、業界全体の最適化を目指す会社です。

Discussion