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

10

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.

mat kelcey
źródło

Odpowiedzi:

11

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} )
Cam.Davidson.Pilon
źródło
2
Bezwstydna promocja: Właśnie napisałem artykuł na blogu o Bayesie i pyMC dosłownie 1 minutę przed opublikowaniem tego, więc zapraszam do sprawdzenia. The Awesome Power of Bayes - Part 1
Cam.Davidson.Pilon
niesamowite! to podejście do mieszania tych dwóch środków jest dokładnie tym, co próbowałem uzyskać.
mat kelcey
Nie jestem pewien, czy w pełni rozumiem prawdziwą korzyść z modelowania, mówiącą, że średnia1 i średnia2 są zwykle dystrybuowane zamiast jednolitego (to samo dotyczy precyzji, jeśli mam być szczery, używam Gammy od „kogoś innego”). Muszę się wiele nauczyć :)
mat kelcey
Używanie munduru, jak w twoim oryginalnym przykładzie, oznacza, że ​​wiesz z absolutną pewnością, że średnia nie przekracza pewnej wartości. To jest trochę patologiczne. Lepiej jest użyć normalnej, ponieważ pozwala na uwzględnienie wszystkich liczb rzeczywistych.
Cam.Davidson.Pilon
1
Wybór gamma ma matematyczny powód. Gamma jest sprzężona przed precyzją, patrz tabela tutaj
Cam.Davidson.Pilon
6

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 .

Fonnesbeck
źródło
1
Ah Fonnesbeck jest jednym z głównych twórców pymc! Czy możesz nam pokazać przykład kodowania punktu 2?
Cam.Davidson.Pilon
dzięki Fonnesbeck i tak, proszę! do szybkiego np. punktu 2 :)
mat kelcey
1
tak naprawdę zgaduję, że masz na myśli coś w stylu ... gist.github.com/4404631 ?
mat kelcey
Tak, dokładnie. Możesz dokonać transformacji nieco bardziej zwięźle:tau = std_dev**-2
fonnesbeck,
jakie byłoby właściwe miejsce, aby przeczytać o tym, skąd bierze się ta relacja między precyzją a std_dev?
user979