久久久久在线人妻,日本东京热色网站,网红主播亚洲情色

本文鏈接：個人站 | 簡書 | CSDN
版權(quán)聲明：除特別聲明外，本博客文章均采用 BY-NC-SA 許可協(xié)議。轉(zhuǎn)載請注明出處。

背景

我們在《報童問題》、《報童問題的簡單解法》等文中介紹了一種通過考慮需求的不確定性來最大化銷售利潤的商品采購模型：首先預(yù)測需求所服從的概率分布，然后取能使得期望收益最大的分位數(shù)作為預(yù)估的需求，據(jù)此來決定采購量，對應(yīng)的分位值定義為服務(wù)水平。

在實際應(yīng)用中，一次采購需要滿足未來一段時間的總需求，具體是多長時間取決于商品的提前期和目標(biāo)庫轉(zhuǎn)等因素。理論上我們可以直接預(yù)測這段時間的總需求所服從的概率分布。但站在甲方的角度，一段時間的總需求？還概率分布？沒概念??！你不告訴我每天的情況，直接丟一個最終結(jié)果給我，我怎么知道你靠不靠譜呢？

甲方有這樣的需求無疑是十分合理的，作為乙方我們沒有理由不去滿足。為此，我們可以預(yù)測每一天的需求所服從的概率分布，然后計算總需求所服從的概率分布。如下圖所示。

fig1.gif

考慮到可控性，需要允許甲方人為調(diào)整服務(wù)水平。站在甲方的角度看，問題又來了，我調(diào)整服務(wù)水平的時候，只能看到預(yù)估的總需求量在變，我想知道對應(yīng)的每天的需求量是怎么變的。這就不好講了，您想啊，總需求多 10 件：有可能是第一天多了 5 件，第二天多了 3 件，第三天多了 2 件；也有可能是第一天多了 1 件，第二天多了 3 件，第三天多了 6 件……可能的情況多了去了。不出意外的話，甲方這時候就會問了，在這么多的情況中，你能不能告訴我哪一種的可能性最高呢？你不是概率分布預(yù)測嗎？算算概率唄。

甲方有這樣的需求無疑是十分合理的，作為乙方我們沒有理由不去滿足。為此，我們需要求解一個最大化條件概率的問題。

問題

考慮一組獨立的隨機變量 $X_1, X_2, \cdots, X_n$ ，它們各自所服從的概率分布已知。
令
$X = (X_1, X_2, \cdots, X_n)\\ Z = \sum_{i=1}^{n} X_i$
給定一個數(shù) $z$ ，求使得條件概率 $P(X=\vec x|Z=z)$ 最大的 $\vec x$ 的取值，即
$\vec x^* = \max_{\vec x} P(X=\vec x|Z=z)$

求解

我們用 $f_i(x_i)$ 表示隨機變量 $X_i$ 的概率密度函數(shù)或概率質(zhì)量函數(shù)，則
$\begin{aligned} P(X=x|Z=z) &\propto f_n\left(z-\sum_{i=1}^{n-1}x_i\right)\cdot\prod_{i=1}^{n-1}f_i(x_i)\\ &\equiv g(x_1, x_2, \cdots, x_{n-1}) \end{aligned}$
定義
$\eta(x) =\frac{f'(x)}{f(x)}$
令
$\begin{aligned} \frac{\partial g}{\partial x_i} &= 0\\ &= -f_n'\left(z-\sum_{j=1}^{n-1}x_j\right)\cdot\prod_{j=1}^{n-1}f_j(x_j) + f_n\left(z-\sum_{j=1}^{n-1}x_j\right)\cdot\prod_{j=1, j\neq i}^{n-1}f_j(x_j)\cdot f_i'(x_i)\\ &= -\eta_n\left(z-\sum_{j=1}^{n-1}x_j\right)\cdot f_n\left(z-\sum_{j=1}^{n-1}x_j\right)\cdot\prod_{j=1}^{n-1}f_j(x_j) + f_n\left(z-\sum_{j=1}^{n-1}x_j\right)\cdot\prod_{j=1, j\neq i}^{n-1}f_j(x_j)\cdot f_i(x_i)\cdot\eta_i(x_i)\\ &= g\cdot(-\eta_n\left(z-\sum_{j=1}^{n-1}x_j\right) + \eta_i(x_i)) \end{aligned}$
得
$\eta_i(x_i) = \eta_n\left(z-\sum_{j=1}^{n-1}x_j\right)$
也就是說，原問題的解 $\vec x^*=(x_1, x_2, \cdots,x_n)$ 滿足
$\begin{cases} x_1 + x_2 + \cdots + x_n = z\\ \eta_1(x_1)=\eta_2(x_2)=\cdots=\eta_n(x_n) \end{cases}$
對于正態(tài)分布，有
$\eta(x) = \frac{\mu-x}{\sigma^2}$
如果隨機變量 $X_1, X_2, \cdots, X_n$ 均服從正態(tài)分布，則只需要求解線性方程組
$\left[ \begin{matrix} 1 & 1 & \cdots & 1 & 1\\ -1/\sigma_1^2 & 0 & \cdots & 0 & 1/\sigma_n^2\\ 0 & -1/\sigma_2^2 & \cdots & 0 & 1/\sigma_n^2\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & -1/\sigma_{n-1}^2 & 1/\sigma_n^2 \end{matrix} \right] \left[ \begin{matrix} x_1\\x_2\\x_3\\\vdots\\x_n \end{matrix} \right] = \left[ \begin{matrix} z\\ \mu_n/\sigma_n^2-\mu_1/\sigma_1^2\\ \mu_n/\sigma_n^2-\mu_2/\sigma_2^2\\ \vdots\\ \mu_n/\sigma_n^2-\mu_{n-1}/\sigma_{n-1}^2\\ \end{matrix} \right]$
即可得到 $\vec x^*$ 。

遺憾的是，銷量通常并不服從正態(tài)分布。我們來考慮一下銷量預(yù)測中常用的分布形式。

對于泊松分布，有
$\eta(x) = \ln\lambda - \psi(x+1)$
對于二項分布，有
$\eta(x) = \psi(n-x+1) - \psi(x+1) + \ln\frac{p}{1-p}$
對于負(fù)二項分布，有
$\eta(x) = \psi(x+r) - \psi(x+1) + \ln p$

其中 $\psi(x)$ 為 digamma 函數(shù)，定義為
$\psi(x) = \frac{\mathrm d}{\mathrm dx}\ln\left(\Gamma(x)\right)=\frac{\Gamma'(x)}{\Gamma(x)}$

對于這些分布， $\eta(x)$ 都是非線性的，無法通過求解線性方程組的方式來計算 $\vec x^*$ ?？紤]到

對于泊松分布，有
$\eta'(x) = -\psi_1(x+1) <0$
對于二項分布，有
$\eta'(x) = -\psi_1(n-x+1) - \psi_1(x+1) <0$
對于負(fù)二項分布，有
$\eta'(x) = \psi_1(x+r) - \psi_1(x+1) <0\qquad if \quad r > 1$
其中 $\psi_1(x)$ 為 trigamma 函數(shù)，定義為
$\psi_1(x) = \frac{\mathrm d^2}{\mathrm dx^2}\ln\left(\Gamma(x)\right)=\frac{\mathrm d}{\mathrm dx}\psi(x)$

也就是說這些分布的 $\eta(x)$ 都是單調(diào)遞減的。因此只要令
$\begin{aligned} \eta_{min} &= \max\left(\eta_1\left(z\right), \eta_2\left(z\right), \cdots, \eta_2\left(z\right)\right)\\ \eta_{max} &= \max\left(\eta_1\left(\frac zn\right), \eta_2\left(\frac zn\right), \cdots, \eta_2\left(\frac zn\right)\right) \end{aligned}$
就可以使用二分法求得 $\eta^*\in[\eta_{min}, \eta_{max}]$ ，使得
$\sum_{i=1}^{n}x_i^*=\sum_{i=1}^{n}\eta_i^{-1}(\eta^*) = z$
從而求得 $\vec x^*$ 。如下圖所示。

fig2.png

這里的問題在于我們并不知道 $\eta(x)$ 的反函數(shù) $\eta^{-1}(y)$ 的解析形式。好在同樣可以使用二分法來求解（其實我原本用的是牛頓法，但實驗中發(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):
    """
    二項分布
    """
    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ù)二項分布
    """
    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ù)二項分布定義為成功概率為 p 的伯努利試驗失敗 r 次時成功次數(shù)所服從的分布
        # 而 scipy 中的定義則是成功概率為 p 的伯努利試驗成功 r 次時失敗次數(shù)所服從的分布
        # 因此在調(diào)用 scipy 中的相關(guān)函數(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

我們用一個可以解析求解的例子來驗證一下代碼。設(shè) $X_1\sim N(0, 1^2)$ ， $X_2\sim N(0, 4^2)$ ， $z=5$ 。根據(jù)前面的推導(dǎo)，只需要求解線性方程組：
$\left[ \begin{matrix} 1 & 1\\ -1 & 1/16 \end{matrix} \right] \left[ \begin{matrix} x_1\\x_2 \end{matrix} \right] = \left[ \begin{matrix} 5\\0 \end{matrix} \right]$
用高斯消元法解得 $x_1=5/17$ 、 $x_2=80/17$ 。我們來看看數(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é)果。

最后用一個復(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)
]

fig3.gif

附錄

附上 $\eta(x)$ 的推導(dǎo)供感興趣的同學(xué)參考。

首先是正態(tài)分布，有
$f(x) = \frac{1}{\sigma\sqrt 2 \pi}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$
故
$\begin{aligned} \eta(x) &= \frac{f'(x)}{f(x)}\\ &= \frac{\mathrm d}{\mathrm dx}\ln f(x)\\ &= \frac{\mathrm d}{\mathrm dx}\left(\ln\frac{1}{\sigma\sqrt 2 \pi}-\frac{(x-\mu)^2}{2\sigma^2}\right)\\ &= \frac{\mathrm d}{\mathrm dx}\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\\ &= \frac{\mu-x}{\sigma^2} \end{aligned}$
其次是泊松分布，有
$f(x) = \frac{\mathrm e^{-\lambda}\lambda^x}{x!} = \frac{\mathrm e^{-\lambda}\lambda^x}{\Gamma(x+1)}$
故
$\begin{aligned} \eta(x) &= \frac{f'(x)}{f(x)}\\ &= \frac{\mathrm d}{\mathrm dx}\ln f(x)\\ &= \frac{\mathrm d}{\mathrm dx}\left(-\lambda+x\ln\lambda-\ln\Gamma(x+1)\right)\\ &= \ln\lambda-\frac{\mathrm d}{\mathrm dx}\ln\Gamma(x+1)\\ &= \ln\lambda-\psi(x+1) \end{aligned}$
接著是二項分布，有
$\begin{aligned} f(x) & = \tbinom{n}{x}p^x(1-p)^{n-x}\\ &=\frac{n!}{x!(n-x)!}p^x(1-p)^{n-x}\\ &= \frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)}p^x(1-p)^{n-x} \end{aligned}$
故
$\begin{aligned} \eta(x) &= \frac{f'(x)}{f(x)}\\ &= \frac{\mathrm d}{\mathrm dx}\ln f(x)\\ &= \frac{\mathrm d}{\mathrm dx}\left(\ln\Gamma(n+1)-\ln\Gamma(x+1)-\ln\Gamma(n-x+1)+x\ln p+(n-x)\ln (1-p)\right)\\ &= \frac{\mathrm d}{\mathrm dx}\left(-\ln\Gamma(x+1)-\ln\Gamma(n-x+1)+x\ln p-x\ln (1-p)\right)\\ &= -\frac{\mathrm d}{\mathrm dx}\ln\Gamma(n-x+1)-\frac{\mathrm d}{\mathrm dx}\ln\Gamma(x+1) + \ln\frac{p}{1-p}\\ &= \psi(n-x+1) - \psi(x+1) + \ln\frac{p}{1-p} \end{aligned}$
最后是負(fù)二項分布，有
$\begin{aligned} f(x) &= \tbinom{x+r-1}{x}p^x(1-p)^r\\ &=\frac{(x+r-1)!}{x!(r-1)!}p^x(1-p)^r\\ &=\frac{\Gamma(x+r)}{\Gamma(x+1)\Gamma(r)}p^x(1-p)^r\\ \end{aligned}$
故
$\begin{aligned} \eta(x) &= \frac{f'(x)}{f(x)}\\ &= \frac{\mathrm d}{\mathrm dx}\ln f(x)\\ &= \frac{\mathrm d}{\mathrm dx}\left(\ln\Gamma(x+r)-\ln\Gamma(x+1)-\ln\Gamma(r)+x\ln p+r\ln (1-p)\right)\\ &= \frac{\mathrm d}{\mathrm dx}\left(\ln\Gamma(x+r)-\ln\Gamma(x+1)+x\ln p\right)\\ &= \psi(x+r) - \psi(x+1) + \ln p \end{aligned}$