Как смоделировать сумму двух распределений Пуассона с помощью pymc 2?

Я пытаюсь смоделировать простой пример вероятностного программирования, используя pymc 2. Я играл с другими языками, такими как церковный и англиканский, и могу без труда смоделировать эту проблему. Однако я не могу понять это в Python.

Вот код на англиканском, я думаю, это довольно очевидно:

[assume a (- (poisson 100) 100)]
[assume b (- (poisson 100) 100)]
[observe (normal (+ a b) .00001) 7]
[predict (list a b)]

Используя сэмплер Metropolis-Hastings, я получаю:

   1 (10 1)
   2 (10 8)
9977 (7 0)
  20 (7 1)

С Particle Gibbs я получаю:

 669 (-1 8)
  71 (-10 17)
  66 (-11 18)
 208 (-12 19)
  19 (-13 20)
  84 (-14 21)
  72 (-15 22)
 441 (-2 9)
...and so on...

Я пытаюсь смоделировать это в pymc так:

def make_model():
    a = (pymc.Poisson("a", 100) - 100)
    b = (pymc.Poisson("b", 100) - 100)

    precision = pymc.Uniform('precision', lower=.0001, upper=1.0)

    @pymc.deterministic
    def mu(a=a, b=b):
        return a+b

    y = pymc.Normal("y", mu=mu, tau=precision, observed=True, value=7)

    return pymc.Model(locals())

def run_mcmc(model):
    mcmc = pymc.MCMC(model)
    mcmc.sample(5000, burn=1000, thin=2)
    return mcmc

result = run_mcmc(make_model())
pymc.Matplot.plot(result)

Я получаю трассировки, где a и b составляют около 100. Однако, если я запускаю (pymc.Poisson("a", 100) - 100).value, я получаю числа ближе к 0.

Я что-то упустил здесь? Я в восторге от возможностей, но в данный момент очень сбит с толку! Спасибо за любую помощь!


python probability pymc
person stratospark    schedule 10.02.2015    source источник
comment
Не могли бы вы подробнее описать вывод англиканских моделей? Я не знаком с этой системой.   -  person Abraham D Flaxman    schedule 10.02.2015
comment
Да, это в основном подсчет частоты. Например, 9977 (7 0) означает, что a=7, b=0 появлялись 9977 раз из 10000 образцов.   -  person stratospark    schedule 11.02.2015

Ответы (1)


arrow_upward
2
arrow_downward

Если я правильно понимаю, это хороший пример, демонстрирующий некоторые различия в мышлении между англиканцами и PyMC.

Вот измененная версия вашего кода PyMC, которая, я думаю, отражает ваше намерение:

def make_model():
    a = pymc.Poisson("a", 100)  # better to have the stochastics themselves available
    b = pymc.Poisson("b", 100)

    precision = 1e-4**-2 #  Seems like precision is fixed in Anglican model (and different from the meaning of precision in PyMC)

    @pymc.deterministic
    def mu(a=a, b=b):
        return (a-100) + (b-100)

    y = pymc.Normal("y", mu=mu, tau=precision, observed=True, value=7)

    return pymc.Model(locals())

def run_mcmc(model):
    mcmc = pymc.MCMC(model)
    mcmc.use_step_method(pymc.AdaptiveMetropolis, [mcmc.a, mcmc.b])
    mcmc.sample(20000, burn=10000, thin=10)
    return mcmc

result = run_mcmc(make_model())
pymc.Matplot.plot(result)

Вот основные отличия моего кода:

  • a и b являются стохастиками. PyMC делает что-то слишком умное для своего же блага, когда вы используете (stochastic - 100) в своей версии.
  • precision — это число, а не стохастик, и большое, а не маленькое число. Это связано с тем, что PyMC использует точность для обозначения 1/дисперсии в нормальном распределении, но в англиканском языке (я думаю) точность означает, насколько точно вам нужен оператор равенства.
  • mcmc использует адаптивный пошаговый метод Metropolis с более длительным периодом приработки. Это важно, потому что совместное апостериорное распределение a и b имеет экстремальную корреляцию, и шаги MCMC никуда не денутся, пока это не будет выяснено.

Вот блокнот IPython, который показывает немного больше деталей.

person Abraham D Flaxman    schedule 12.02.2015