2013-07-01 12 views
7

Załóżmy, że otrzymaliśmy wcześniejszą wiadomość od X (na przykład X ~ Gaussian) i operatora przesyłania do przodu y = f (x). Przypuśćmy, że za pomocą eksperymentu zaobserwowaliśmy y i że ten eksperyment może być powtarzany w nieskończoność. Zakłada się, że wyjście Y ma charakter gaussowski (Y-Gaussian) lub wolne od szumów (Delta Y (obserwacja)).Rozwiązywanie odwrotnych problemów z PyMC

Jak konsekwentnie aktualizować nasz subiektywny poziom wiedzy o X biorąc pod uwagę obserwacje? Próbowałem następujący model z PyMC, ale wydaje mi czegoś brakuje:

from pymc import * 

xtrue = 2      # this value is unknown in the real application 
x = rnormal(0, 0.01, size=10000) # initial guess 

for i in range(5): 
    X = Normal('X', x.mean(), 1./x.var()) 
    Y = X*X      # f(x) = x*x 
    OBS = Normal('OBS', Y, 0.1, value=xtrue*xtrue+rnormal(0,1), observed=True) 
    model = Model([X,Y,OBS]) 
    mcmc = MCMC(model) 
    mcmc.sample(10000) 

    x = mcmc.trace('X')[:]  # posterior samples 

Niedrażniącym nie jest zbieżny do xtrue.

Odpowiedz

7

Funkcjonalność stworzona przez @ user1572508 jest teraz częścią PyMC pod nazwą stochastic_from_data() lub Histogram(). Rozwiązaniem tego wątku staje się:

from pymc import * 
import matplotlib.pyplot as plt 

xtrue = 2 # unknown in the real application 
prior = rnormal(0,1,10000) # initial guess is inaccurate 
for i in range(5): 
    x = stochastic_from_data('x', prior) 
    y = x*x 
    obs = Normal('obs', y, 0.1, xtrue*xtrue + rnormal(0,1), observed=True) 

    model = Model([x,y,obs]) 
    mcmc = MCMC(model) 
    mcmc.sample(10000) 

    Matplot.plot(mcmc.trace('x')) 
    plt.show() 

    prior = mcmc.trace('x')[:] 
5

Problem polega na tym, że twoja funkcja, $ y = x^2 $, nie jest jeden-do-jednego. W szczególności, ponieważ tracisz wszystkie informacje o znaku X podczas jego skracania, nie ma możliwości odróżnienia wartości Y od tego, czy pierwotnie użyłeś 2 czy -2 do wygenerowania danych. Jeśli utworzysz histogram śledzenia X po pierwszej iteracji, zobaczysz: histogram of trace after first iteration

Rozkład ten ma 2 tryby: jeden przy 2 (twoja prawdziwa wartość) i jedna przy -2. Przy następnej iteracji x.mean() będzie bliska zeru (uśredniając rozkład bimodalny), co oczywiście nie jest tym, czego potrzebujesz.

+0

Wiem, że f (x) nie jest biempionem, zostało wybrane z tego konkretnego powodu. Nie widzę żadnego argumentu, aby MCMC zawodziło z tą dystrybucją wejściową, o ile wiem, MCMC jest w stanie wypróbować złożone, wielomodalne dystrybucje. – juliohm

+0

Och, zrozumiałem twoją sprawę, problem polega na tym, jak aktualizuję stan sprzed. Używając po prostu pos.mean() i pos.var() zakładam rozwiązanie unimodalne. Jak rozwiązałbyś ten problem, znajdując zarówno 2, jak i -2? – juliohm

+1

Innymi słowy, jak reprezentować nieparametryczne dystrybucje za pomocą PyMC? Biorąc pod uwagę ślad, przygotuj plik PDF pasujący do histogramu. – juliohm

Powiązane problemy