我想生成二项分布的混合。为什么我需要它是因为我想要高斯分布的正态离散混合。是否有任何可用的 scipy 库,或者请您指导我的算法。
我通常知道对于预定义的发行版,可以使用 ppf。但是对于这个功能,我认为没有任何直接的使用 ppf 的方法。
从每个样本中取样并混合它们似乎也有问题,因为我不知道我必须从不同的分布中选择多少个实例。
最后我想要的是这样的:
我想生成二项分布的混合。为什么我需要它是因为我想要高斯分布的正态离散混合。是否有任何可用的 scipy 库,或者请您指导我的算法。
我通常知道对于预定义的发行版,可以使用 ppf。但是对于这个功能,我认为没有任何直接的使用 ppf 的方法。
从每个样本中取样并混合它们似乎也有问题,因为我不知道我必须从不同的分布中选择多少个实例。
最后我想要的是这样的:
这是生成二项式(和其他)分布的任意混合的简单方法。它依赖于一个事实,如果你想从混合物 P(x)=sum(w[i]*P_i(x), i=1..Nmix) 中获取样本 (Nsamp),那么你可以通过采样来做到这一点来自每个 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))
除非您找到一种计算逆 cdf 的聪明方法(在这种情况下请告诉我们!),否则拒绝采样是一种万无一失的方法。维基百科条目给出了一般描述。我在实践中发现,你需要对“工具”分布有点小心:特别是它不应该比目标分布衰减得快得多——如果是这样,你可能会失去尾部的贡献.
我这样做的方式是从一个平坦的工具分布开始:生成一对均匀的随机数x
和y
,其中y
来自 [0, 1) 和x
来自[0, L)
,其中L
足够大。然后比较y
和cdf(x)
,重复直到收敛。如果这行得通,你就准备好了。如果这还不够好,请使用非平坦的工具分布:如果混合物的尾部是高斯分布,那么您最好使用高斯分布。
附带说明一下,如果您正在处理二项式分布,则需要注意上溢/下溢 --- 根据参数,您可能需要使用高斯近似。
感谢@sega_sai、@askewchan 和@Zhenya,我自己编写了代码,我相信由于实现,这将是最有效的代码。有两个函数,第一个函数使“binoNumber”二项式分布的混合都具有相同的 N=maximum-minimum 参数和相同的 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')