0

私は単純なモデルを解決するためにpymcを使用しようとしています:

  • パレート分布から引き出されたことがわかっている N=1000 フラックスがあります: フラックス ~ パレート(アルファ, 1)

  • パレートのアルファパラメータを計算しようとしています: alpha ~ Uniform(1, 3)

  • フラックスの測定値がガウス ノイズによって汚染されています。flux_meas ~ Gaussian(flux, tau)

現在、ノイズの量を変更するとどうなるかをシミュレートしようとしています。

問題は、非常に少量の (つまり無視できる) 量のノイズでモデルを実行すると、アルファの平均値が実行ごとに根本​​的に変化し、アルファの真の値とはまったく関係がないように見えることです。しかし、ガウス ノイズを完全に無視して、パレート分布を直接観察しただけだとすると、期待どおりに機能します。

私は何を間違っていますか？

動作するコード スニペットはこちら: (より複雑な古いものはこちら)

キービットは次のとおりです。

import pymc
N = 1000
true_alpha = 2
noise = 0.001 # This noise is much smaller than the signal

# Simulated fluxes
s_arr = pymc.rpareto(true_alpha, 1, size=N)

# the unknown alpha
alpha = pymc.Uniform('alpha', 1, 3)

# fluxes are drawn from a Pareto distribution
flux = pymc.Pareto('flux', alpha, 1, size=N)

# My observed fluxes are contaminated by Gaussian noise
flux_meas = pymc.Normal('flux_meas', mu=flux, tau=noise**-2, observed=True,
                         value=pymc.rnormal(s_arr, tau=noise**-2, size=N))

model = pymc.MCMC([alpha, flux, flux_meas])

# If I run this model several times, the mean of alpha will be somewhere between
# 1 and 3. The variance of alpha is pretty small
model.sample(5000, 1000, 5)
4

1 に答える 1