модель гауссовой смеси с биномиальными весами - PullRequest
Новая
       5

модель гауссовой смеси с биномиальными весами

1 голос
/ 15 января 2020

В моем проекте вывода Байеса у меня есть вероятность, которая выглядит следующим образом:

equation1

с equation2

Я должен оценить параметры λ , σ нормалей и n, p биномиального распределения по данным. Если n фиксировано, я могу сделать вывод. 

import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import theano.tensor as tt
import scipy.special
import arviz as az

true_p = 0.35
n = 10
N_data = 75
true_mean = 10
true_sigma = 4
s = np.random.binomial(n,true_p,N_data)
mean = [true_mean * (i + 1) for i in range(n + 1)]
data = np.concatenate(
    [np.random.normal(mean[i], true_sigma, len(s[s==i])) for i in range(n + 1)])

def Bin_pmf(n, p):
    return [scipy.special.binom(n, i) * (p**i) * (1 - p)**(n - i) for i in range(n + 1)]

with pm.Model() as fixed_n_model:
    p = pm.Beta("p", 2, 2)
    mu = pm.Uniform('mu', 0, 20)
    sigma = pm.HalfCauchy('sigma', 2.5)
    mu1 = tt.stack([(i + 1) * mu for i in range(n + 1)])
    components = pm.Normal.dist(mu=mu1, sigma=sigma, shape=(n + 1))
    like = pm.Mixture('like', w=tt.stack(Bin_pmf(n, p)),
                      comp_dists=components, observed=data)

with fixed_n_model:
    start = pm.find_MAP()
    start["p"] = 0.6
    start["mu"] = 15
    step1 = pm.NUTS([sigma, p, mu], max_treedepth=15, target_accept=0.95)
    trace = pm.sample(10000, start=start, step=[step1])
pm.traceplot(trace, lines=[('p', {}, [true_p]), ("mu", {}, [
             true_mean]), ("sigma", {}, [true_sigma])])
az.plot_posterior(trace)
ppc = pm.sample_posterior_predictive(trace, samples=500, model=fixed_n_model)
fig1, ax1 = plt.subplots()
ax1.hist(np.asarray(ppc['like'])[:, 1],round(N_data/2))
ax1.hist(data,bins=round(N_data/2))

Но у этого решения есть параметр n в аргументе формы, который не может быть случайным в pymc3.

Я пытался чтобы обойти это, используя смесь распределений смеси, но это не кажется хорошей идеей. 

k = 12
n = 10
l = 2
n_dist = [n, 2]
with pm.Model() as fixed_n_model:
    a = [tt.vector() for _ in range(l)]
    b = [tt.vector() for _ in range(l)]
    p = pm.Beta("p", 2, 2)
    mu = pm.Uniform('mu', 0, 20)
    n = pm.DiscreteUniform("n", 0, 1, shape=2)
    for i in range(l):

        t = [(j + 1) for j in range(n_dist[i] + 1)] + \
            [0 for i in range(k - n_dist[i] - 1)]
        a[i] = tt.stack([j * mu for j in t])
        b[i] = tt.stack(Bin_pmf(n_dist[i], p) +
                        [0 for j in range(k - n_dist[i] - 1)])
    sigma = pm.HalfCauchy('sigma', 2)
    components_1 = pm.Normal.dist(mu=a[0], sigma=sigma, shape=k)
    components_2 = pm.Normal.dist(mu=a[1], sigma=sigma, shape=k)
    like_1 = pm.Mixture.dist(w=b[0],
                             comp_dists=components_1)
    like_2 = pm.Mixture.dist(w=b[1],
                             comp_dists=components_2)
    like = pm.Mixture('like', w=n,
                      comp_dists=[like_1, like_2], observed=data)
...