# 一个最大化条件概率问题

fig1.gif

## 求解

• 对于泊松分布，有

• 对于二项分布，有

• 对于负二项分布，有

• 对于泊松分布，有

• 对于二项分布，有

• 对于负二项分布，有

其中 为 trigamma 函数，定义为

fig2.png

## 代码

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


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

fig3.gif

• 首先是正态分布，有

• 其次是泊松分布，有

• 接着是二项分布，有

• 最后是负二项分布，有