CODEWITHSUNDEEP
pythonbayesianpymc3
aloctavodia
aloctavodia

Reputation: 2070

How to define a custom prior in PyMC3

I would like to know if it is possible to define a custom prior in PyMC3 (and how to do it). From here it seems that in PyMC2 is relatively easy to do (without the need to modified the source code), but in PyMC3 is not that easy (or I am not understanding something). I am trying to replicate a prior from the book "Doing Bayesian Data Analysis", that is implemented in BUGS:

model {
# Likelihood. Each flip is Bernoulli.
for ( i in 1 : N1 ) { y1[i]  ̃ dbern( theta1 ) }
for ( i in 1 : N2 ) { y2[i]  ̃ dbern( theta2 ) }
# Prior. Curved scallo not ps!
x  ̃ dunif(0,1)
y  ̃ dunif(0,1)
N <- 4
xt <- sin( 2*3.141593*N * x ) / (2*3.141593*N) + x
yt <- 3 * y + (1/3)
xtt <- pow( xt , yt )
theta1 <- xtt
theta2 <- y
}

The prior does not make too much sense, it just an example of how to define a custom prior and the versatility of BUGS.

My attempt to implement the above custom prior is:

from __future__ import division
import numpy as np
import pymc as pm
from pymc import Continuous
from theano.tensor import sin, log

# Generate the data
y1 = np.array([1, 1, 1, 1, 1, 0, 0])  # 5 heads and 2 tails
y2 = np.array([1, 1, 0, 0, 0, 0, 0])  # 2 heads and 5 tails

class Custom_prior(Continuous): 
"""
custom prior
"""
    def __init__(self, y, *args, **kwargs):
        super(Custom_prior, self).__init__(*args, **kwargs)
        self.y = y
        self.N = 4
        self.mean = 0.625  # FIXME
    def logp(self, value):
        N = self.N
        y = self.y
        return -log((sin(2*3.141593*N * value)
                     / (2*3.141593*N) + value)**(3 * y + (1/3)))

with pm.Model() as model:
    theta2 = pm.Uniform('theta2', 0, 1)  # prior
    theta1 = Custom_prior('theta1', theta2)  # prior
    # define the likelihood
    y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
    y2 = pm.Bernoulli('y2', p=theta2, observed=y2)
    # Generate a MCMC chain
    start = pm.find_MAP()  # Find starting value by optimization
    trace = pm.sample(5000, pm.NUTS(), progressbar=False)

EDIT

Following the answer of chris-fonnesbeck

I think I need something like:

with pm.Model() as model:
    theta2 = pm.Uniform('theta2', 0, 1)  # prior
    N = 4
    theta1 = pm.DensityDist('theta1', lambda value: -log((sin(2*3.141593*N * value)
                       / (2*3.141593*N) + value)**(3 * theta2 + (1/3))))
    # define the likelihood
    y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
    y2 = pm.Bernoulli('y2', p=theta2, observed=y2)

    # Generate a MCMC chain
    start = pm.find_MAP()  # Find starting value by optimization
    trace = pm.sample(10000, pm.NUTS(), progressbar=False) # Use NUTS sampling

The only problem is that I get the same value for all the posterior samples of theta1 and theta2, I guess there is some problem with my custom prior or the combination of priors and likelihoods. A successful definition of a custom prior can be found in this example

Upvotes: 3

Views: 2619

Answers (1)

Chris Fonnesbeck
Chris Fonnesbeck

Reputation: 4203

Can you post the full BUGS model? In the above, it just looks like a series of deterministic transformations in BUGS after the priors for x and y, rather than the definition of a prior.

Assuming the logp above is what you want, though, you can implement it in PyMC much more simply as:

def logp(value, y):
    N  = 4
    return -log((sin(2*3.141593*N * value)
                 / (2*3.141593*N) + value)**(3 * y + (1/3)))

theta1 = pm.DensityDist('theta1', logp, value, y=theta2)

Upvotes: 4

Related Questions