Модель подгонки для двух нормальных распределений в PyMC

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

Я изучал PyMC и работал над некоторыми (очень) простыми примерами. Одна проблема, которую я не могу получить (и не могу найти связанных примеров), - это подгонка модели к данным, сгенерированным из двух нормальных распределений.

Скажем, у меня есть 1000 значений; 500 генерируется из а Normal(mean=100, stddev=20)и еще 500 генерируется из а Normal(mean=200, stddev=20).

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

mean1 = Uniform('mean1', lower=0.0, upper=200.0)
mean2 = Uniform('mean2', lower=0.0, upper=200.0)
precision = Gamma('precision', alpha=0.1, beta=0.1)

data = read_data_from_file_or_whatever()

@deterministic(plot=False)
def mean(m1=mean1, m2=mean2):
    # but what goes here?

process = Normal('process', mu=mean, tau=precision, value=data, observed=True)

процесс генерации нормальный, но мю - одно из двух значений. Я просто не знаю, как представить «решение» между значением m1или значением m2.

Возможно, я просто совершенно неверно подходил к моделированию этого? Кто-нибудь может указать мне на пример? Я могу читать ОШИБКИ и ЯБЛОКИ, так что все в порядке на самом деле.

Вы абсолютно уверены, что половина пришла из одного дистрибутива, а другая - из другого? Если нет, мы можем смоделировать пропорцию как случайную величину (что очень байесовская вещь).

Ниже я бы сделал несколько советов.

from pymc import *

size = 10
p = Uniform( "p", 0 , 1) #this is the fraction that come from mean1 vs mean2

ber = Bernoulli( "ber", p = p, size = size) # produces 1 with proportion p.

precision = Gamma('precision', alpha=0.1, beta=0.1)

mean1 = Normal( "mean1", 0, 0.001 ) #better to use normals versus Uniforms (unless you are certain the value is  truncated at 0 and 200 
mean2 = Normal( "mean2", 0, 0.001 )

@deterministic
def mean( ber = ber, mean1 = mean1, mean2 = mean2):
    return ber*mean1 + (1-ber)*mean2


#generate some artificial data   
v = np.random.randint( 0, 2, size)
data = v*(10+ np.random.randn(size) ) + (1-v)*(-10 + np.random.randn(size ) )


obs = Normal( "obs", mean, precision, value = data, observed = True)

model = Model( {"p":p, "precision": precision, "mean1": mean1, "mean2":mean2, "obs":obs} )

Пара моментов, связанных с обсуждением выше:

1. Выбор диффузного нормального или униформного является довольно академическим, если (а) вы не беспокоитесь о сопряженности, в этом случае вы бы использовали нормальное или (б) есть некоторый разумный шанс, что истинное значение может быть за пределами конечных точек униформы , С PyMC нет причин беспокоиться о сопряженности, если только вы не хотите использовать сэмплер Gibbs.

2. Гамма на самом деле не лучший выбор для неинформативного до параметра дисперсии / точности. Это может оказаться более информативным, чем вы думаете. Лучше выбрать стандартное отклонение, а затем преобразовать его в обратный квадрат. См. Гельман 2006 для деталей.