本文鏈接：個(gè)人站 | 簡(jiǎn)書(shū) | CSDN
版權(quán)聲明：除特別聲明外赃梧，本博客文章均采用 BY-NC-SA 許可協(xié)議峡钓。轉(zhuǎn)載請(qǐng)注明出處峦耘。
背景
我們?cè)?a href="http://www.reibang.com/p/dfac740faa34" target="_blank">《報(bào)童問(wèn)題》、《報(bào)童問(wèn)題的簡(jiǎn)單解法》等文中介紹了一種通過(guò)考慮需求的不確定性來(lái)最大化銷售利潤(rùn)的商品采購(gòu)模型：首先預(yù)測(cè)需求所服從的概率分布杆麸，然后取能使得期望收益最大的分位數(shù)作為預(yù)估的需求，據(jù)此來(lái)決定采購(gòu)量浪感，對(duì)應(yīng)的分位值定義為服務(wù)水平昔头。
在實(shí)際應(yīng)用中，一次采購(gòu)需要滿足未來(lái)一段時(shí)間的總需求影兽，具體是多長(zhǎng)時(shí)間取決于商品的提前期和目標(biāo)庫(kù)轉(zhuǎn)等因素揭斧。理論上我們可以直接預(yù)測(cè)這段時(shí)間的總需求所服從的概率分布。但站在甲方的角度峻堰，一段時(shí)間的總需求讹开？還概率分布？沒(méi)概念凹攵省萧吠！你不告訴我每天的情況，直接丟一個(gè)最終結(jié)果給我桐筏，我怎么知道你靠不靠譜呢纸型？
甲方有這樣的需求無(wú)疑是十分合理的，作為乙方我們沒(méi)有理由不去滿足梅忌。為此狰腌，我們可以預(yù)測(cè)每一天的需求所服從的概率分布，然后計(jì)算總需求所服從的概率分布牧氮。如下圖所示琼腔。
考慮到可控性，需要允許甲方人為調(diào)整服務(wù)水平踱葛。站在甲方的角度看丹莲，問(wèn)題又來(lái)了光坝，我調(diào)整服務(wù)水平的時(shí)候，只能看到預(yù)估的總需求量在變甥材，我想知道對(duì)應(yīng)的每天的需求量是怎么變的盯另。這就不好講了，您想啊洲赵，總需求多 10 件：有可能是第一天多了 5 件鸳惯，第二天多了 3 件，第三天多了 2 件叠萍；也有可能是第一天多了 1 件芝发，第二天多了 3 件，第三天多了 6 件……可能的情況多了去了苛谷。不出意外的話辅鲸，甲方這時(shí)候就會(huì)問(wèn)了，在這么多的情況中抄腔，你能不能告訴我哪一種的可能性最高呢瓢湃？你不是概率分布預(yù)測(cè)嗎？算算概率唄赫蛇。
甲方有這樣的需求無(wú)疑是十分合理的绵患，作為乙方我們沒(méi)有理由不去滿足。為此悟耘，我們需要求解一個(gè)最大化條件概率的問(wèn)題落蝙。
問(wèn)題
考慮一組獨(dú)立的隨機(jī)變量 ，它們各自所服從的概率分布已知暂幼。
令
給定一個(gè)數(shù) 筏勒，求使得條件概率
最大的
的取值，即
求解
我們用 表示隨機(jī)變量
的概率密度函數(shù)或概率質(zhì)量函數(shù)旺嬉，則
定義
令
得
也就是說(shuō)管行，原問(wèn)題的解 滿足
對(duì)于正態(tài)分布，有
如果隨機(jī)變量 均服從正態(tài)分布邪媳，則只需要求解線性方程組
即可得到 捐顷。
遺憾的是，銷量通常并不服從正態(tài)分布雨效。我們來(lái)考慮一下銷量預(yù)測(cè)中常用的分布形式迅涮。
對(duì)于泊松分布，有
對(duì)于二項(xiàng)分布徽龟，有
對(duì)于負(fù)二項(xiàng)分布叮姑，有
其中 為 digamma 函數(shù)，定義為
對(duì)于這些分布据悔， 都是非線性的传透，無(wú)法通過(guò)求解線性方程組的方式來(lái)計(jì)算
耘沼。考慮到
對(duì)于泊松分布旷祸，有
對(duì)于二項(xiàng)分布耕拷，有
對(duì)于負(fù)二項(xiàng)分布，有
其中為 trigamma 函數(shù)托享，定義為
也就是說(shuō)這些分布的 都是單調(diào)遞減的。因此只要令
就可以使用二分法求得 浸赫，使得
從而求得 闰围。如下圖所示。
這里的問(wèn)題在于我們并不知道 的反函數(shù)
的解析形式既峡。好在同樣可以使用二分法來(lái)求解（其實(shí)我原本用的是牛頓法羡榴，但實(shí)驗(yàn)中發(fā)現(xiàn)存在一些難以收斂的情況，故改用二分法）运敢。
代碼
from abc import ABC, abstractmethod
import numpy as np
from scipy.stats import norm, poisson, binom, nbinom
from scipy.special import digamma
class Distribution(ABC):
"""
概率分布基類
"""
@property
@abstractmethod
def mu(self):
"""
概率分布的期望
"""
raise NotImplementedError()
@abstractmethod
def quantile(self, q):
"""
概率分布的 q 分位數(shù)
"""
raise NotImplementedError()
@abstractmethod
def eta(self, x):
"""
$\eta(x)$
"""
raise NotImplementedError()
def ieta(self, y, x_min, x_max):
"""
$\eta^{-1}(y)$
"""
# 默認(rèn)使用二分法求解
while True:
x_mid = (x_min + x_max) / 2
diff = self.eta(x_mid) - y
if np.abs(diff) <= 1e-6:
break
if diff > 0:
x_min = x_mid
else:
x_max = x_mid
if x_max == x_min:
raise ValueError('Unable to solve')
return x_mid
class Normal(Distribution):
"""
正態(tài)分布
"""
def __init__(self, mu, sigma):
super().__init__()
self._mu = mu
self._sigma = sigma
@property
def mu(self):
return self._mu
def quantile(self, q):
return norm.ppf(q=q, loc=self._mu, scale=self._sigma)
def eta(self, x):
return (self._mu - x) / self._sigma ** 2
def ieta(self, y, x_mid=None, x_max=None):
# 用解析解法覆蓋父類的數(shù)值解法
return self._mu - y * self._sigma ** 2
class Poisson(Distribution):
"""
泊松分布
"""
def __init__(self, mu):
super().__init__()
self._mu = mu
@property
def mu(self):
return self._mu
def quantile(self, q):
return poisson.ppf(q=q, mu=self._mu)
def eta(self, x):
return np.log(self._mu) - digamma(x + 1)
class Binomial(Distribution):
"""
二項(xiàng)分布
"""
def __init__(self, n, p):
super().__init__()
self._n = n
self._p = p
@property
def mu(self):
return self._n * self._p
def quantile(self, q):
return binom.ppf(q=q, n=self._n, p=self._p)
def eta(self, x):
return digamma(self._n - x + 1) - digamma(x + 1) + np.log(self._p/(1 - self._p))
def ieta(self, y, x_min, x_max):
return super().ieta(y, x_min, min(x_max, self._n))
class NegativeBinomial(Distribution):
"""
負(fù)二項(xiàng)分布
"""
def __init__(self, r, p):
super().__init__()
self._r = r
self._p = p
@property
def mu(self):
return self._r * self._p / (1 - self._p)
def quantile(self, q):
# 我們將負(fù)二項(xiàng)分布定義為成功概率為 p 的伯努利試驗(yàn)失敗 r 次時(shí)成功次數(shù)所服從的分布
# 而 scipy 中的定義則是成功概率為 p 的伯努利試驗(yàn)成功 r 次時(shí)失敗次數(shù)所服從的分布
# 因此在調(diào)用 scipy 中的相關(guān)函數(shù)時(shí)需要注意轉(zhuǎn)換
return nbinom.ppf(q=q, n=self._r, p=1-self._p)
def eta(self, x):
return digamma(x + self._r) - digamma(x + 1) + np.log(self._p)
def max_posterior(distrs, z):
"""
用二分法求解使條件概率 $P(X=\vec x|Z=z)$ 最大的 $\vec_x^*$
Parameters
----------
distrs : List<Distribution>
概率分布列表
z : float
總和的目標(biāo)值 $z$
Returns
-------
List<float>
$\vec_x^*$
"""
n = len(distrs)
e_max = max(d.eta(z/n) for d in distrs)
e_min = max(d.eta(z) for d in distrs)
while True:
e_mid = (e_min + e_max) / 2
xs = [d.ieta(e_mid, 0, z) for d in distrs]
z_hat = sum(xs)
if np.abs(z_hat - z) <= 1e-2:
break
if z_hat < z:
e_max = e_mid
else:
e_min = e_mid
if e_max - e_min < 1e-6:
raise ValueError('Unable to solve')
return xs
我們用一個(gè)可以解析求解的例子來(lái)驗(yàn)證一下代碼校仑。設(shè) ，
传惠，
迄沫。根據(jù)前面的推導(dǎo)，只需要求解線性方程組：
用高斯消元法解得 卦方、
羊瘩。我們來(lái)看看數(shù)值解：
>>> distrs = [Normal(0, 1), Normal(0, 4)]
>>> print(max_posterior(distrs, 5))
[0.294189453125, 4.70703125]
>>>
>>> print(5/17)
0.29411764705882354
>>>
>>> print(80/17)
4.705882352941177
可以看到，數(shù)值解法給出了與解析解非常接近的結(jié)果盼砍。
最后用一個(gè)復(fù)雜的例子直觀地感受一下效果：
distrs = [
Poisson(20),
Poisson(19),
Poisson(18),
Poisson(19),
Binomial(40, 0.2),
Binomial(45, 0.25),
NegativeBinomial(60, 0.3),
NegativeBinomial(60, 0.25),
Poisson(21),
Poisson(20)
]
附錄
附上 的推導(dǎo)供感興趣的同學(xué)參考尘吗。
首先是正態(tài)分布，有
故
其次是泊松分布浇坐，有
故
接著是二項(xiàng)分布睬捶，有
故
最后是負(fù)二項(xiàng)分布，有
故
完近刘。
