How to use StudentT distribution in pymc3?

Question

I'm not sure whether this counts as a question or a bug report. I posted a GitHub gist here: https://gist.github.com/jbwhit/a9012e04b0f48e582c22

I found this question (pymc3: hierarchical model with multiple obsesrved variables) to be an excellent starting point for my own hierarchical model, but ran into difficulties as soon as I tried to modify it in any substantial way.

First, the model and setup that works:

import numpy as np
import pymc3 as pm

n_individuals = 200
points_per_individual = 10
means = np.random.normal(30, 12, n_individuals)
observed = np.random.normal(means, 1, (points_per_individual, n_individuals))

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu=0, sd=100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd=3)

    means = pm.Normal('means', mu=hyper_mean, sd=hyper_sigma, shape=n_individuals)
    sigmas = pm.HalfNormal('sigmas', sd=100)

    ye = pm.Normal('ye', mu=means, sd=sigmas, observed=observed)

    trace = pm.sample(10000)

All of the above works as expected (and the traces look nice). The next piece of code makes one change (swapping a T distribution for the Normal):

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu=0, sd=100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd=3)

    ### Changed to a T distribution ###
    means = pm.StudentT('means', nu=hyper_mean, sd=hyper_sigma, shape=n_individuals)

    sigmas = pm.HalfNormal('sigmas', sd=100)

    ye = pm.Normal('ye', mu=means, sd=sigmas, observed=observed)

    trace = pm.sample(10000)

The following is the output:

Assigned NUTS to hyper_mean
Assigned NUTS to hyper_sigma_log
Assigned NUTS to means
Assigned NUTS to sigmas_log
---------------------------------------------------------------------------
PositiveDefiniteError                     Traceback (most recent call last)
<ipython-input-12-69f59e2f3d47> in <module>()
     18     ye = pm.Normal('ye', mu=means, sd=sigmas, observed=observed)
     19 
---> 20     trace = pm.sample(10000)

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py in sample(draws, step, start, trace, chain, njobs, tune, progressbar, model, random_seed)
    121     """
    122     model = modelcontext(model)
--> 123 
    124     step = assign_step_methods(model, step)
    125 

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py in assign_step_methods(model, step, methods)
     66             selected_steps[selected].append(var)
     67 
---> 68     # Instantiate all selected step methods
     69     steps += [s(vars=selected_steps[s]) for s in selected_steps if selected_steps[s]]
     70 

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py in <listcomp>(.0)
     66             selected_steps[selected].append(var)
     67 
---> 68     # Instantiate all selected step methods
     69     steps += [s(vars=selected_steps[s]) for s in selected_steps if selected_steps[s]]
     70 

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/nuts.py in __init__(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model, profile, **kwargs)
     76 
     77 
---> 78         self.potential = quad_potential(scaling, is_cov, as_cov=False)
     79 
     80         if state is None:

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/quadpotential.py in quad_potential(C, is_cov, as_cov)
     33         return QuadPotential_SparseInv(C)
     34 
---> 35     partial_check_positive_definite(C)
     36     if C.ndim == 1:
     37         if is_cov != as_cov:

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/quadpotential.py in partial_check_positive_definite(C)
     56     if len(i):
     57         raise PositiveDefiniteError(
---> 58             "Simple check failed. Diagonal contains negatives", i)
     59 
     60 

PositiveDefiniteError: Scaling is not positive definite. Simple check failed. Diagonal contains negatives. Check indexes [202]

Any suggestion on how to get this to work?

Answer 1

As I mentioned in the comment, try running:

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu = 0, sd = 100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd = 3)
    nu = pm.Exponential('nu', 1./10, testval = 5.)

    ### Changed to a T distribution ###
    means = pm.StudentT('means', nu = nu, mu = hyper_mean, sd = hyper_sigma, shape = n_individuals)

    sigmas = pm.HalfNormal('sigmas', sd = 100)

    ye = pm.Normal('ye', mu = means, sd = sigmas, observed = observed)

    trace = pm.sample(10000)

In other words: use the mu argument of the pm.StudentT method for hyper_mean and nu for the degrees of freedom.

Once it starts working, you might also try to add the pm.find_MAP method (as suggested by @Chris Fonnesbeck).

Answer 2

Try finding the MAP estimate and use that as the starting point for the MCMC run:

start = pm.find_MAP()
trace = pm.sample(10000, start=start)