Создание смеси биномиальных распределений

Question

Я хочу создать смесь биномиального распределения. Почему мне это нужно, потому что Я хочу иметь нормальную дискретную смесь гауссовских дистрибутивов. Есть ли какая-нибудь библиотека scipy, доступная для нее, или можете ли вы, пожалуйста, направить меня на алгоритм.

Я знаю, что для предопределенных распределений можно использовать ppf. Но для этой функции я не думаю, что есть простой способ использования ppf.

Отбор проб из каждого и их смешение также кажется проблематичным, потому что я не знаю, как Множество экземпляров Мне нужно выбирать из разных дистрибутивов.

В конце концов, что я хочу, чтобы иметь что-то вроде этого:

Answer 1

Благодаря @ sega_sai, @ askewchan и @Zhenya, я сделал код сам и я believ из-за реализации это было бы самым эффективным. Есть две функции, первая из которых делает смеси биномиальных распределений «binoNumber» одинаковыми с N = максимальным минимальным параметром и одинаковыми p = 0,5, но сдвигаются в зависимости от созданных для них случайных центров I.

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

Answer 2

Если вы не найти умный способ вычисления обратной CDF (! В этом случае дайте нам знать), отбор проб отказ является надежный способ. А wikipedia entry содержит описание модели. То, что я нашел на практике, вам нужно быть немного осторожным с «инструментальным» распределением: в частности, он не должен распадаться намного быстрее, чем целевой дистрибутив, - если это так, вы, вероятно, потеряете вклад хвостов ,

То, как я бы это сделать, я бы начать с плоским инструментальным распределения: генерировать пару однородных случайных чисел x и y, где y находятся в диапазоне от [0, 1) и x от [0, L), где L достаточно большой. Затем сравните y и cdf(x), повторите до конвергенции. Если это сработает, вы все настроены. Если это не достаточно хорошо, используйте не плоское инструментальное распределение: если хвост смеси гауссовский, вам, вероятно, лучше всего использовать гауссовский.

В качестве побочного примечания, если вы имеете дело с биномиальным распределением, вам нужно следить за переходом/недополнением --- в зависимости от параметров вам может понадобиться гауссовское приближение.

Answer 3

Вот простой способ генерации произвольных смесей биномиальных (и других) распределений. Он полагается на то, что если вы хотите получить образцы (Nsamp) из смеси P (x) = sum (w [i] * P_i (x), i = 1..Nmix), то вы можете сделать это путем выборки Nsamp от каждого из P_i (x). Затем получим еще один образец Nsamp другой случайной величины, который равен i с вероятностью w [i]. Эта случайная переменная может быть использована для выбора, какой из P_i (х) данный образец будет приходить:

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