Giải quyết các vấn đề nghịch đảo với PyMC

Question

Giả sử chúng tôi được đưa ra trước trên X (ví dụ: X ~ Gaussian) và toán tử chuyển tiếp y = f (x). Giả sử chúng tôi đã quan sát thêm y bằng thử nghiệm và thử nghiệm này có thể được lặp lại vô thời hạn. Đầu ra Y được giả định là Gaussian (Y ~ Gaussian) hoặc không có nhiễu (Y ~ Delta (quan sát)).

Làm cách nào để cập nhật thường xuyên mức độ hiểu biết chủ quan của chúng tôi về X với các quan sát? Tôi đã thử các mô hình sau đây với PyMC, nhưng có vẻ như tôi là thiếu một cái gì đó:

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 

Các sau không được hội tụ để xtrue.

Answer 1

Chức năng được nhắm mục tiêu bởi @ user1572508 giờ đây là một phần của PyMC dưới tên stochastic_from_data() hoặc Histogram(). Giải pháp cho chủ đề này sau đó trở thành:

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

Vấn đề là chức năng của bạn, $ y = x^2 $, không phải là một đối một. Cụ thể, bởi vì bạn mất tất cả thông tin về dấu X khi bạn tạo hình vuông, không có cách nào để nói từ các giá trị Y của bạn cho dù ban đầu bạn đã sử dụng 2 hoặc -2 để tạo dữ liệu. Nếu bạn tạo biểu đồ theo dõi của mình cho X chỉ sau lần lặp đầu tiên, bạn sẽ thấy điều này:

Phân phối này có 2 chế độ, một tại 2 (giá trị thực) và một ở -2. Ở lần lặp tiếp theo, x.mean() sẽ gần bằng không (trung bình trên phân phối hai mặt), rõ ràng không phải là thứ bạn muốn.