Risoluzione problemi inversi con PyMC

Question

Supponiamo siamo in una precedente su X (ad esempio X ~ gaussiana) e un attaccante operatore y = f (x). Supponiamo di aver osservato y per mezzo di un esperimento e che questo esperimento possa essere ripetuto indefinitamente. Si presume che l'uscita Y sia gaussiana (Y ~ gaussiana) o priva di disturbi (Y ~ Delta (osservazione)).

Come aggiornare costantemente il nostro grado di conoscenza soggettiva su X date le osservazioni? Ho provato il seguente modello con PyMC, ma sembra che mi manca qualcosa:

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 

Il posteriore non sta convergendo a Xtrue.

Answer 1

La funzionalità di @ user1572508 è ora parte di PyMC con il nome stochastic_from_data() o Histogram(). La soluzione a questo thread diventa:

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

Il problema è che la funzione, $ y = x^2 $, non è uno-a-uno. In particolare, poiché perdi tutte le informazioni sul segno di X quando lo piazza, non c'è modo di dire dai tuoi valori Y se hai originariamente usato 2 o -2 per generare i dati. Se crei un istogramma della traccia per X dopo la prima iterazione, vedrai:

Questa distribuzione ha 2 modalità, una a 2 (il tuo valore reale) e una a -2. Alla successiva iterazione, x.mean() sarà vicino a zero (facendo una media sulla distribuzione bimodale), che ovviamente non è quello che vuoi.