pythonbayesianpymcpymc3

How to set Bernoulli distribution parameters in pymc3


I have a model described in pymc3 using the following:

from pymc3 import * 
basic_model = Model()

with basic_model:
    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=18)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2 + beta[2]*X3

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

However, my Ys are not normally distributed but are binary (so, Bernoulli, I think). I can't figure out how to change the Normal distrubtion of Y to Bernoulli though because I can't figure out what the params would be of Y_obs in that case.


Solution

  • What you are looking for is logistic regression. Here you use the logistic function to convert the output of your linear model to a probability.

    In your example this could be specified as follows:

    import pyMc3 as pm
    import theano.tensor as T
    basic_model = pm.Model()
    
    def logistic(l):
        return 1 / (1 + T.exp(-l))
    
    with basic_model:
        # Priors for unknown model parameters
        alpha = pm.Normal('alpha', mu=0, sd=10)
        beta = pm.Normal('beta', mu=0, sd=10, shape=18)
    
        # Expected value of outcome
        mu = alpha + beta[0]*X1 + beta[1]*X2 + beta[2]*X3
    
        # Likelihood (sampling distribution) of observations
        Y_obs = pm.Bernoulli('Y_obs', p=logistic(mu), observed=Y)