生成二项分布混合体

Question

我想生成一个二项分布的混合体。之所以需要这个，是因为我想要一个正常的离散高斯分布混合体。请问有没有现成的scipy库可以用，或者能不能给我讲讲这个算法。

我知道一般来说，对于预定义的分布，可以使用ppf函数。但对于这个功能，我觉得没有简单的方法来使用ppf。

从每个分布中抽样并进行混合似乎也有问题，因为我不知道应该从不同的分布中选择多少个实例。

最后，我想要的结果是这样的：

Answer 1

除非你找到了一种聪明的方法来计算逆累积分布函数（如果找到了，请告诉我们！），否则拒绝采样是一种可靠的方法。维基百科上有一个关于拒绝采样的介绍，可以了解基本概念。在实际操作中，我发现你需要对“工具”分布稍微小心：具体来说，它的衰减速度不应该比目标分布快太多——如果太快，你可能会失去尾部的贡献。

我会这样做：从一个平坦的工具分布开始：生成一对均匀随机数 x 和 y，其中 y 在 [0, 1) 范围内，x 在 [0, L) 范围内，L 要足够大。然后比较 y 和 cdf(x)，重复这个过程直到收敛。如果这样可以工作，那就没问题。如果效果不够好，可以使用一个非平坦的工具分布：如果混合分布的尾部是高斯分布，使用高斯分布可能是最好的选择。

另外，如果你在处理二项分布时，需要注意溢出和下溢——根据参数的不同，你可能需要使用高斯近似。

Answer 2

这里有一个简单的方法，可以生成任意组合的二项分布（以及其他分布）。这个方法基于一个事实：如果你想从一个混合分布中获取样本（称为Nsamp），这个混合分布可以表示为P(x)=sum(w[i]*P_i(x), i=1..Nmix)，那么你可以从每个P_i(x)中抽取Nsamp个样本。接着，再从一个随机变量中抽取Nsamp个样本，这个随机变量的值为i的概率是w[i]。这个随机变量可以用来选择你要从哪个P_i(x)中获取样本：

import numpy as np,numpy.random, matplotlib.pyplot as plt

#parameters of the binomial distributions: pairs of (n,p)
binomsP = np.array([.5, .5, .5])
binomsCen = np.array([15, 45, 95]) # centers of binomial distributions
binomsN = (binomsCen/binomsP).astype(int)

fractions = [0.2, 0.3, 0.5]
#mixing fractions of the binomials
assert(sum(fractions)==1)

nbinoms = len(binomsN)
npoints = 10000
cumfractions = np.cumsum(fractions)
def mapper(x):
    # convert the random number between 0 and 1 to
    # the ID of the distribution according to the mixing fractions
    return np.digitize(x, cumfractions)

x0 = np.random.binomial(binomsN[None, :],
        binomsP[None, :], size=(npoints, nbinoms))

x = x0[:, mapper(np.random.uniform(size=npoints))]
plt.hist(x, bin=150, range=(0, 150))

Answer 3

感谢@sega_sai、@askewchan和@Zhenya， 我自己写了这段代码，我相信由于实现方式，这将是最有效的。这里有两个函数，第一个函数是用“binoNumber”个二项分布来生成混合，这些分布都有相同的N=最大值-最小值参数，并且p=0.5，但它们是根据我为它们生成的随机中心进行偏移的。

global binoInitiated
binoInitiated=False;
def binoMixture(minimum,maximum,sampleSize):
    global centers
    binoNumber=10;
    if (not binoInitiated):
        centers=np.random.randint(minimum,maximum+1,binoNumber)
    sigma=maximum-minimum-2
    sam=np.array([]);
    while sam.size<sampleSize:
        i=np.random.choice(binoNumber);
        temp=np.random.binomial(sigma, 0.5,1)+centers[i]-sigma/2+1
        sam=np.append(sam,temp)
    return sam

这个函数是用来绘制之前生成的分布的近似概率密度函数（PDF）。感谢@EnricoGiampieri，我使用了他的代码来完成这一部分。

def binoMixtureDrawer(minimum,maximum):
    global binoInitiated
    global centers
    sam=binoMixture(minimum,maximum,50000)    
    # this create the kernel, given an array it will estimate the probability over that values
    kde = gaussian_kde( sam )
    # these are the values over wich your kernel will be evaluated
    dist_space = linspace( min(sam), max(sam), 500 )
    # plot the results
    fig.plot( dist_space, kde(dist_space),'g')