Как я могу включить категориальные распределения в качестве наблюдаемых данных в PyMC3?

Question

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

Необработанные наблюдаемые данные у:

[[  1.93542253e-01   1.39657327e-04   7.53918636e-01   5.23994535e-02]
 [  6.44964587e-02   8.50087384e-01   1.09894665e-02   7.44266910e-02]
 [  1.68387463e-02   5.38121456e-01   6.98554551e-02   3.75184342e-01]
 ..., 
 [  5.79786789e-01   1.47417427e-02   3.15395731e-01   9.00757372e-02]
 [  8.66796124e-02   8.66999904e-02   4.47848127e-02   7.81835584e-01]
 [  8.18765043e-01   3.23448859e-03   5.61247840e-04   1.77439220e-01]]

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

Выборка из исходных данных y_choice:

[ 2.  3.  3.  1.  2.  2.  2.  2.  3.  3.  1.  2.  3.  0.  2.  0.  3.  1.
  3.  0.  2.  0.  3.  0.  2.  0.  1.  0.  3.  0.  2.  0.  0.  0.  3.  0.
  2.  0.  0.  3.  3.  1. ...

Однако, похоже, что я теряю информацию, потому что мой выборочный сэмплер находится за пределами модели PyMC. Как я могу предоставить фактические данные наблюдений y непосредственно в модель? Я предполагаю, что это как-то связано с другим параметром модели, основанным на Дирихле, но я не смог обернуть голову, как это работает.

Пример кода, с которым я работаю, приведен ниже. Я хочу быть в состоянии предоставить y модели и вернуть true_probs, но пока мне удалось заставить его работать только с y_choice.

import numpy as np
from pymc3 import *
import pymc3 as pm
import pandas as pd
print 'pymc3 version: ' + pm.__version__


def generate_noisy_observations():
    y = np.ones((sample_size,k))
    for i in range(sample_size):
        #print("Iteration %d" % i)
        true_category = np.random.choice(k, size=1, p=true_probs)
        true_distribution = np.zeros(k)
        true_distribution[true_category] = 1
        noise_distribution = np.random.dirichlet(np.ones(k))

        noise = np.random.normal(0,1,k)

        distribution_weights = [0.9, 0.1]
        raw_distribution = (true_distribution*distribution_weights[0] + noise**2*distribution_weights[1] )/\       (np.sum(true_distribution*distribution_weights[0])+np.sum(noise**2*distribution_weights[1]))
        y[i] = raw_distribution
return y

def generate_choices_from_noisy_observations(y):
    y_choice = np.ones(sample_size)
    for i in range(sample_size):
        y_choice[i] = np.random.choice(k, size=1, p=y[i])
    return y_choice

sample_size = 1000

true_probs = [0.2, 0.1, 0.3, 0.4]
k = len(true_probs)


y = generate_noisy_observations()

y_choice = generate_choices_from_noisy_observations(y)

with pm.Model() as multinom_test:
    probs = pm.Dirichlet('a', a=np.ones(k))
    #data = Categorical('data',p = probs, observed = y)
    data = Categorical('data',p = probs, observed = y_choice)
    start = pm.find_MAP()
    trace = pm.sample(50000, pm.Metropolis())

pm.traceplot(trace[500:])