首先,我相信你遇到的问题是因为你错误地对概率进行了归一化处理。以下这行代码是错误的:
a = np.exp(l) / scipy.misc.logsumexp(l)
你正在将一个概率除以一个对数概率,这是没有意义的。相反,你可能想要
a = np.exp(l - scipy.misc.logsumexp(l))
如果您这样做,您会发现a = [1, 0]
,并且您的多项式采样器可以按照第二个概率的浮点精度正常工作。
针对小N的解决方案:直方图
话虽如此,如果您仍需要更高的精度,并且性能不是太大的问题,您可以通过从头开始实现一个多项式采样器,然后修改它以在更高精度上工作。
NumPy的多项式函数是Cython实现的,本质上是对一定数量的二项式样本进行循环,并将它们组合成一个多项式样本。您可以像这样调用它:
np.random.multinomial(10, [0.1, 0.2, 0.7])
(注意,这里以及以下的精确输出值都是随机的,并且会在每次调用时发生变化。)
另一种实现多项式采样器的方法是生成N个均匀随机值,然后使用由累积概率定义的直方图来计算它们:
def multinomial(N, p):
rand = np.random.uniform(size=N)
p_cuml = np.cumsum(np.hstack([[0], p]))
p_cuml /= p_cuml[-1]
return np.histogram(rand, bins=p_cuml)[0]
multinomial(10, [0.1, 0.2, 0.7])
有了这个方法,我们可以考虑通过保持 一切 在对数空间中来实现更高的精度。主要的技巧是意识到均匀随机偏差的对数等于指数随机偏差的负数,因此您可以在不离开对数空间的情况下完成上述所有操作:
def multinomial_log(N, logp):
log_rand = -np.random.exponential(size=N)
logp_cuml = np.logaddexp.accumulate(np.hstack([[-np.inf], logp]))
logp_cuml -= logp_cuml[-1]
return np.histogram(log_rand, bins=logp_cuml)[0]
multinomial_log(10, np.log([0.1, 0.2, 0.7]))
产生的多项式随机抽样将即使在p数组中有非常小的值时仍保持精度。
不幸的是,这些基于直方图的解决方案将比本地的numpy.multinomial
函数慢得多,因此如果性能是一个问题,你可能需要另一种方法。其中一种选择是将上面链接的Cython代码适应为在对数空间中工作,使用类似于我在这里使用的数学技巧。
大N的解决方案:泊松近似
上述解决方案的问题是,随着N的增长,它变得非常缓慢。
我思考了一下,意识到有一种更有效的方法,尽管np.random.multinomial
在概率小于1E-16
左右时会失败。
这是一个失败的例子:在64位机器上,由于代码实现的方式,第一个条目总是为零,而实际上它应该接近10:
np.random.multinomial(1E18, [1E-17, 1])
如果您深入挖掘源代码,您可以追溯到构建多项式函数的二项式函数上出现的这个问题。Cython代码在内部执行类似于以下操作:
def multinomial_basic(N, p, size=None):
results = np.array([np.random.binomial(N, pi, size) for pi in p])
results[-1] = int(N) - results[:-1].sum(0)
return np.rollaxis(results, 0, results.ndim)
multinomial_basic(1E18, [1E-17, 1])
问题在于binomial
函数在非常小的p
值上会出现问题——这是因为该算法计算了值(1-p)
,所以p
的值受到浮点精度的限制。
那么我们该怎么办呢?嗯,事实证明,对于小的p
值,Poisson分布是二项式分布的极好近似,并且该实现没有这些问题。因此,我们可以基于一个稳健的二项式采样器构建一个强大的多项式函数,在小的p
时切换到泊松采样器:
def binomial_robust(N, p, size=None):
if p < 1E-7:
return np.random.poisson(N * p, size)
else:
return np.random.binomial(N, p, size)
def multinomial_robust(N, p, size=None):
results = np.array([binomial_robust(N, pi, size) for pi in p])
results[-1] = int(N) - results[:-1].sum(0)
return np.rollaxis(results, 0, results.ndim)
multinomial_robust(1E18, [1E-17, 1])
array([ 12, 999999999999999988])
第一个条目不为零且接近于10,正如预期!请注意,我们不能使用大于1E18
的N
,因为这将使长整数溢出。但是,我们可以通过使用size
参数对较小概率进行确认,并对结果进行平均来验证我们的方法:
p = [1E-23, 1E-22, 1E-21, 1E-20, 1]
size = int(1E6)
multinomial_robust(1E18, p, size).mean(0)
我们可以看到即使对于这些非常小的概率,多项式值也以正确比例出现。结果是对于小的p
,多项式分布的近似非常稳健且非常快速。