Binomial probit regression, sample code

Sharing
1

Hi
Here is sample code to do probit regression for Binomial or Bernoulli processes.
Based of python - Probit regression using PyMC 3 - Stack Overflow
Updated to use pytensor

import pymc as pm
import pytensor
from pytensor import tensor as pt

def probit_phi(x):
    """ Probit transform assuming 0 mean and 1 sd """
    return (1 + pt.erf(x / pt.sqrt(2)) )/2

# data. Set of N coin-flip observations, with k heads
N_ar=[1000,10]
k_ar=[500,0]

with pm.Model() as model:
    theta_p = pm.Normal('theta_p', mu=0, sigma =100)
    theta = probit_phi(theta_p)

    # likelihood
    y = pm.Binomial('y', p=theta, n=N_ar, observed= k_ar)
    trace= pm.sample()
    summary = pm.stats.summary(trace)
    print(summary)
2

Thanks for sharing!

You can actually skip the custom probit_phi and use pm.math.invprobit, as in:

N_ar=[1000,10]
k_ar=[500,0]

with pm.Model() as model:
    theta_p = pm.Normal('theta_p', mu=0, sigma =100)
    theta = pm.math.invprobit(theta_p)

    # likelihood
    y = pm.Binomial('y', p=theta, n=N_ar, observed= k_ar)
    idata = pm.sample()
3

Excellent! Thanks