Model dopasowania dla dwóch normalnych rozkładów w PyMC

Ponieważ jestem inżynierem oprogramowania i próbuję dowiedzieć się więcej statystyk, musisz mi wybaczyć, zanim zacznę, dlatego jest to poważna nowość ...

Uczę się PyMC i pracuję nad kilkoma naprawdę (naprawdę) prostymi przykładami. Jednym z problemów, których nie mogę zabrać do pracy (i nie mogę znaleźć żadnych powiązanych przykładów), jest dopasowanie modelu do danych wygenerowanych z dwóch normalnych dystrybucji.

Powiedz, że mam 1000 wartości; 500 wygenerowanych z a Normal(mean=100, stddev=20)oraz kolejnych 500 wygenerowanych z Normal(mean=200, stddev=20).

Jeśli chcę dopasować do nich model, tj. Określić dwa średnie i pojedyncze odchylenie standardowe, używając PyMC. Wiem, że to coś w stylu ...

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)

tzn. proces generowania jest normalny, ale mu jest jedną z dwóch wartości. Po prostu nie wiem, jak przedstawić „decyzję” między tym, czy wartość pochodzi, m1czy też m2.

Być może po prostu całkowicie niewłaściwie podchodzę do modelowania tego? Czy ktoś może wskazać mi przykład? Potrafię czytać BŁĘDY i JAGI, więc wszystko jest w porządku.

Czy jesteś absolutnie pewien, że połowa pochodzi z jednej dystrybucji, a druga połowa z drugiej? Jeśli nie, możemy modelować proporcję jako zmienną losową (co jest bardzo bayesowską czynnością).

Oto, co bym zrobił, niektóre wskazówki są osadzone.

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

Kilka punktów związanych z powyższą dyskusją:

1. Wybór rozproszonej normalnej vs. mundurowej jest dość akademicki, chyba że (a) martwisz się koniugacją, w którym to przypadku użyłbyś normalnej lub (b) istnieje uzasadniona szansa, że ​​prawdziwa wartość może znajdować się poza punktami końcowymi munduru . Dzięki PyMC nie ma powodu, aby martwić się o koniugację, chyba że specjalnie chcesz użyć samplera Gibbs.

2. Gamma nie jest tak naprawdę dobrym wyborem dla nieinformacyjnych przed parametrem wariancji / precyzji. Może to być bardziej pouczające, niż myślisz. Lepszym wyborem jest umieszczenie munduru przed odchyleniem standardowym, a następnie przekształcenie go odwrotnym kwadratem. Szczegółowe informacje można znaleźć w Gelman 2006 .