Как да моделирам сумата от две разпределения на Поасон с pymc 2?

Опитвам се да моделирам прост вероятностен пример за програмиране с помощта на pymc 2. Играя си с други езици като Church и Anglican и мога да моделирам този проблем без затруднения. Изглежда обаче не мога да го разбера в 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 пъти от 10 000 проби.   -  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