Lösen inverser Probleme mit PyMC

Question

Angenommen, wir sind gegeben vorveröffentlichten auf X (z.B. X ~ Gaussian) und einen vorderen Operator y = f (x). Nehmen wir weiter an, wir haben y mittels eines Experiments beobachtet und dieses Experiment kann unbegrenzt wiederholt werden. Der Ausgang Y wird als Gauß (Y ~ Gauß) oder rauschfrei (Y ~ Delta (Beobachtung)) angenommen.

Wie aktualisiert man ständig unseren subjektiven Kenntnisgrad über X angesichts der Beobachtungen? Ich habe das folgende Modell mit PyMC versucht, aber es scheint, ich etwas fehlt:

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 

Der hintere nicht konvergierenden ist XTRUE.

Answer 1

Die von @ user1572508 vorgesehene Funktionalität ist jetzt Teil von PyMC unter dem Namen stochastic_from_data() oder Histogram(). Die Lösung für diesen Thread wird dann:

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')[:] 

Answer 2

Das Problem ist, dass Ihre Funktion, $ y = x^2 $, nicht eins zu eins ist. Insbesondere weil Sie alle Informationen über das Vorzeichen von X verlieren, wenn Sie es quadrieren, gibt es keine Möglichkeit, von Ihren Y-Werten zu unterscheiden, ob Sie ursprünglich 2 oder -2 verwendet haben, um die Daten zu generieren. Wenn Sie nach der ersten Iteration ein Histogramm Ihrer Kurve für X erstellen, sehen Sie Folgendes:

Diese Verteilung hat 2 Modi, einen bei 2 (Ihren wahren Wert) und einen bei -2. Bei der nächsten Iteration wird x.mean() nahe bei Null liegen (Mittelung über die bimodale Verteilung), was offensichtlich nicht das ist, was Sie wollen.