Pymc3: NUTS Sampling, Multiple Hierarchical Model, Wishart Distribution

Question

Pymc3: NUTS Sampling, Multiple Hierarchical Model, Wishart Distribution

I am having problems fetching NUTS in a hierarchical model in pymc3. I made some simple models in pymc3 without any problems, but with this assignment I cannot move. Can someone help me?

Data can be obtained here: cars2004.csv . First you need to do some work with the data set:

import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as T

data = pd.read_csv("cars2004.csv", sep=",", index_col= 0)

interest = ["price.retail", "cons.city", "cons.highway", "engine.size", 
            "horsepower","weight", "wheel.base", "length", "width", "ncylinder"]

df1 = data[data.ncylinder != -1]
df1 = df1[df1.fhybrid != "Yes"]
df = df1[interest].copy()
n, p = len(df), 9
df["price.retail"] = df["price.retail"].apply(lambda x: x/1000)

Model in python:

with pm.Model() as model:
    lambd = pm.Gamma("lambd", alpha=1, beta=.005)
    beta0 = pm.Normal("beta0", mu = 80, sd = 100)
    beta = pm.MvNormal("beta", mu = np.zeros(p), tau = lambd*np.eye(p), shape = p)

    mu = pm.MvNormal("mu", mu = np.array([10,10,3,200,1500,250,500,10,5]),
                     cov = np.diag([100,100,100,1000**2,1000**2,100**2,100**2,100**2,100]),
                     shape=p)

    sigmainvert = pm.WishartBartlett("sigmainvert", nu=p,S=.001*np.eye(p),
                                     testval=.001*np.eye(p))

    tau = pm.Gamma("tau", 1, .005)
    #sigma = pm.Uniform("sigma", .1, 100)

    X = pm.MvNormal("X", mu = mu, tau = sigmainvert, shape = (n,p),
                    observed = df[interest[1:]])

    Y = pm.Normal("Y", mu = beta0 + T.dot(X, beta), tau = tau,
                  observed = df["price.retail"])
    dev = pm.Deterministic("dev",-2*Y.logpt)

Values in previous distributions are obtained from assignment. "X" have some missing data, and "sigmainvert" should be Wishart (9, diag (0.001, ..., 0.001).

Using ADVI for NUTS steps and initial values (at first I had some problems with NaN elbo and NaN values, now elbo = -inf). Sample:

with model:
    means, sds, elbo = pm.variational.advi(n=50000, learning_rate=0.1)
    step = pm.NUTS(scaling = model.dict_to_array(sds)**2, is_cov=True)
    trace = pm.sample(1000, step, start = means)  

plot1 = pm.traceplot(trace, varnames=["beta0","beta","tau","mu","dev"],alpha=1)