Перенос кода pymc2 в pymc3: пользовательская функция правдоподобия

Я пытаюсь реализовать пример цензурированных данных в книге Ли и Вагенмакерса (глава 5.5, стр. 70). В pymc2 у меня есть следующая модель:

nattempts = 950   
nfails = 949   
n = 50    # Number of questions
y = np.zeros(nattempts)
y[nattempts-1] = 1
z = 30
unobsmin = 15
unobsmax = 25
unobsrange = np.arange(unobsmin,unobsmax+1)

theta = pymc.Uniform("theta",lower = .25, upper = 1)

@pymc.observed
def Ylike(value=z, theta = theta, n=n, censorn=nfails, unobs=unobsrange):
    ylikeobs = pymc.binomial_like(x=value, n=n, p=theta)
    ylikeunobs = np.array([])
    for i in unobs:
        ylikeunobs = np.append(pymc.binomial_like(x=i, n=n, p=theta),ylikeunobs)
    return ylikeobs+sum(ylikeunobs)*censorn

testmodel = pymc.Model([theta,Ylike])
mcmc = pymc.MCMC(testmodel)
mcmc.sample(iter = 20000, burn = 50, thin = 2)

который включал декоратор @pymc.observed.
Я думаю, что мне нужно выразить вероятность с помощью pm.DensityDist, однако я не мог понять, как это сделать.


python pymc3 pymc
person Junpeng Lao    schedule 22.02.2016    source источник
comment
Этот ответ может помочь указать вам правильное направление. stackoverflow.com/questions/22015055/< /а>   -  person inversion    schedule 26.02.2016

Ответы (1)


arrow_upward
0
arrow_downward

Хорошо, я узнал, как это сделать:

with pm.Model():
    theta = pm.Uniform("theta",lower = .25, upper = 1)
    def logp(value,n,p):
        return pm.dist_math.bound(
        pm.dist_math.binomln(n, value) 
      + pm.dist_math.logpow(p, value) 
      + pm.dist_math.logpow(1 - p, n - value),
        0 <= value, value <= n,
        0 <= p, p <= 1)
    def Censorlike(value=z, n=n, censorn=nfails, unobs=unobsrange):
        ylikeobs = logp(value=value, n=n, p=theta)
        ylikeunobs = 0
        for i in unobs:
            ylikeunobs += logp(value=i, n=n, p=theta)
        return ylikeobs+ylikeunobs*censorn

    ylike = pm.DensityDist('ylike', Censorlike, observed={'value':z,'n':n,'censorn':nfails,'unobs':unobsrange})
    trace = pm.sample(3e3)
person Junpeng Lao    schedule 15.03.2016