一个最大化条件概率问题

背景

我们在《报童问题》、《报童问题的简单解法》等文中介绍了一种通过考虑需求的不确定性来最大化销售利润的商品采购模型：首先预测需求所服从的概率分布，然后取能使得期望收益最大的分位数作为预估的需求，据此来决定采购量，对应的分位值定义为服务水平。

在实际应用中，一次采购需要满足未来一段时间的总需求，具体是多长时间取决于商品的提前期和目标库转等因素。理论上我们可以直接预测这段时间的总需求所服从的概率分布。但站在甲方的角度，一段时间的总需求？还概率分布？没概念啊！你不告诉我每天的情况，直接丢一个最终结果给我，我怎么知道你靠不靠谱呢？

甲方有这样的需求无疑是十分合理的，作为乙方我们没有理由不去满足。为此，我们可以预测每一天的需求所服从的概率分布，然后计算总需求所服从的概率分布。如下图所示。

考虑到可控性，需要允许甲方人为调整服务水平。站在甲方的角度看，问题又来了，我调整服务水平的时候，只能看到预估的总需求量在变，我想知道对应的每天的需求量是怎么变的。这就不好讲了，您想啊，总需求多 10 件：有可能是第一天多了 5 件，第二天多了 3 件，第三天多了 2 件；也有可能是第一天多了 1 件，第二天多了 3 件，第三天多了 6 件……可能的情况多了去了。不出意外的话，甲方这时候就会问了，在这么多的情况中，你能不能告诉我哪一种的可能性最高呢？你不是概率分布预测吗？算算概率呗。

甲方有这样的需求无疑是十分合理的，作为乙方我们没有理由不去满足。为此，我们需要求解一个最大化条件概率的问题。

问题

考虑一组独立的随机变量 $X_1, X_2, \cdots, X_n$，它们各自所服从的概率分布已知。
令

$X = (X_1, X_2, \cdots, X_n)\\ Z = \sum_{i=1}^{n} X_i$

给定一个数 $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$ 的概率密度函数或概率质量函数，则

$\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}$

对于正态分布，有

$\eta(x) = \frac{\mu-x}{\sigma^2}$

如果随机变量 $X_1, X_2, \cdots, X_n$ 均服从正态分布，则只需要求解线性方程组

$\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^*$。

遗憾的是，销量通常并不服从正态分布。我们来考虑一下销量预测中常用的分布形式。

• 对于泊松分布，有

$\eta(x) = \ln\lambda - \psi(x+1)$

• 对于二项分布，有

$\eta(x) = \psi(n-x+1) - \psi(x+1) + \ln\frac{p}{1-p}$

• 对于负二项分布，有

$\eta(x) = \psi(x+r) - \psi(x+1) + \ln p$

其中 $\psi(x)$ 为 digamma 函数，定义为

$\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$

• 对于负二项分布，有

$\eta'(x) = \psi_1(x+r) - \psi_1(x+1) <0\qquad if \quad r > 1$

其中 $\psi_1(x)$ 为 trigamma 函数，定义为

$\psi_1(x) = \frac{\mathrm d^2}{\mathrm dx^2}\ln\left(\Gamma(x)\right)=\frac{\mathrm d}{\mathrm dx}\psi(x)$

也就是说这些分布的 $\eta(x)$ 都是单调递减的。因此只要令

$\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^*$。如下图所示。

这里的问题在于我们并不知道 $\eta(x)$ 的反函数 $\eta^{-1}(y)$ 的解析形式。好在同样可以使用二分法来求解（其实我原本用的是牛顿法，但实验中发现存在一些难以收敛的情况，故改用二分法）。

代码

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 分位数
        """
        raise NotImplementedError()
    
    @abstractmethod
    def eta(self, x):
        """
        $\eta(x)$
        """
        raise NotImplementedError()
    
    def ieta(self, y, x_min, x_max):
        """
        $\eta^{-1}(y)$
        """
        # 默认使用二分法求解
        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):
    """
    正态分布
    """
    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):
        # 用解析解法覆盖父类的数值解法
        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):
    """
    负二项分布
    """
    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):
        # 我们将负二项分布定义为成功概率为 p 的伯努利试验失败 r 次时成功次数所服从的分布
        # 而 scipy 中的定义则是成功概率为 p 的伯努利试验成功 r 次时失败次数所服从的分布
        # 因此在调用 scipy 中的相关函数时需要注意转换
        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
        总和的目标值 $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

我们用一个可以解析求解的例子来验证一下代码。设 $X_1\sim N(0, 1)$，$X_2\sim N(0, 4)$，$z=5$。根据前面的推导，只需要求解线性方程组：

$\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$。我们来看看数值解：

>>> 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

可以看到，数值解法给出了与解析解非常接近的结果。

最后用一个复杂的例子直观地感受一下效果：

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)
]

附录

附上 $\eta(x)$ 的推导供感兴趣的同学参考。

• 首先是正态分布，有

$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}$

• 最后是负二项分布，有

$\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}$

完。

参考文献

• Normal distribution - Wikipedia
• Poisson distribution - Wikipedia
• Binomial distribution - Wikipedia
• Negative binomial distribution - Wikipedia
• Gamma function - Wikipedia
• Digamma function - Wikipedia
• Trigamma function - Wikipedia