Implementation of hierarchical bayesian model produces errors with PyMC

I'm trying to solve an exercise from this book in which I'm supposed to fit data on temperature and elevation in Colorado to this model:

$$\boldsymbol{Y} = \boldsymbol{\mu} + \boldsymbol{\alpha} \otimes \boldsymbol{1}_{n\times 1} + \boldsymbol{1}_{T\times 1} \otimes \boldsymbol{w} + \boldsymbol{\epsilon}$$

$$p(\boldsymbol{Y}| \boldsymbol{\beta},\boldsymbol{\alpha},\boldsymbol{w},\sigma^{2}_\epsilon) \sim N(\boldsymbol{\mu} + \boldsymbol{\alpha} \otimes \boldsymbol{1}_{n\times 1} + \boldsymbol{1}_{T\times 1} \otimes \boldsymbol{w}, \sigma^{2}_\epsilon I_{T n \times T n})$$

$\boldsymbol{Y}$ is a vector $(Y_{1}, \dots , Y_{T})$ of observations (temperatures) at each time $t \in \{1,2,\dots, T\}$ where each $Y_{t} = (Y(s_{1}, t), \dots, Y(s_{n},t))$ describes a vector of observations in some location $s_{i}$ at time $t$. Therefore, we can join all observations for every single time and location in a $nT\times 1$ vector $\boldsymbol{Y}$:

$$\boldsymbol{Y} = \left(Y(s_{1}, t_{1}), Y(s_{2}, t_{1}), \dots, Y(s_{n}, t_{1}), \dots ,Y(s_{1}, t_{T}), Y(s_{2}, t_{T}), \dots Y(s_{n}, t_{T})\right)^{T}$$

$\boldsymbol{\mu}$ is similar to $\boldsymbol{Y}$ and contains means that depend on a covariate in the form $\boldsymbol{\mu} = X(s_i, t)\boldsymbol{\beta}(t)$. On the contrary, $\boldsymbol{\alpha}$ only depends on time $\boldsymbol{\alpha} = (\alpha(1),\dots, \alpha(T))$ and $\boldsymbol{\omega}$ only depends on space $\boldsymbol{\omega} = (\omega(1),\dots, \omega(s_{n}))$. That is why we need to use kronecker products to obtain a $nT\times 1$ vector.

I wrote some code with PyMC to reproduce this WinBUGS code in which $\mu = \beta_{1} + \beta_{2}\text{Elevation}$. I think I obtained reasonable results. Now I want to extend this and use a $\boldsymbol{\beta}$ for each time $t$. This means the equation above should be:

$$\left(\begin{array}{ccc}Y(s_{1}, t_{1}) \\ Y(s_{2}, t_{1}) \\ \vdots \\Y(s_{n}, t_{1}) \\ \vdots \\ Y(s_{1}, t_{T}) \\ Y(s_{2}, t_{T}) \\ \vdots \\Y(s_{n}, t_{T})\end{array}\right) = \left(\begin{array}{ccc}\beta_{0t_{1}} + \beta_{1t_{1}}\text{Elevation}(s_{1}) \\ \beta_{0t_{1}} + \beta_{1t_{1}}\text{Elevation}(s_{2}) \\ \vdots \\\beta_{0t_{1}} + \beta_{1t_{1}}\text{Elevation}(s_{n}) \\ \vdots \\ \beta_{0t_{T}} + \beta_{1t_{T}}\text{Elevation}(s_{1}) \\ \beta_{0t_{T}} + \beta_{1t_{T}}\text{Elevation}(s_{2}) \\ \vdots \\ \beta_{0t_{T}} + \beta_{1t_{T}}\text{Elevation}(s_{n})\end{array}\right) + \left(\begin{array}{ccc}\alpha(1) \\ \alpha(1) \\ \vdots \\ \alpha(1) \\ \vdots \\ \alpha(T) \\ \alpha(T) \\ \vdots \\ \alpha(T)\end{array}\right) + \left(\begin{array}{ccc}\omega(s_{1}) \\ \omega(s_{2}) \\ \vdots \\ \omega(s_{n}) \\ \vdots \\ \omega(s_{1}) \\ \omega(s_{2}) \\ \vdots \\ \omega(s_{n})\end{array}\right) + \left(\begin{array}{ccc}\epsilon \\ \epsilon \\ \vdots \\ \epsilon \\ \vdots \\ \epsilon \\ \epsilon \\ \vdots \\ \epsilon\end{array}\right)$$

The code is the following:

import pandas as pd
import numpy as np
import pymc as pm
from scipy.spatial.distance import pdist, squareform
from scipy.linalg import inv

# data

elevation = np.array(pd.concat([df[['Elevation']]]*12))
temperature = np.array(df.ix[:, 3:]).reshape(-1)

ntime = 12 # number of months
nsite = 50 # number of measurement sites

# error term epsilon
error_prec = pm.Gamma('error_prec', 0.10, 0.10, value=1.0)
error_prec_m = pm.Lambda('error_prec_m', lambda error_prec=error_prec: np.diag([error_prec]*nsite*ntime))

# intercept and coefficient of covariate in mu.
beta_1 = pm.Normal('beta_1', 0.0, 0.0001, value=[0.0]*ntime, size=ntime)
beta_2 = pm.Normal('beta_2', 0.0, 0.0001, value=[0.0]*ntime, size=ntime)

# precision and variance for covariance matrix associated to time
temporal_prec = pm.Gamma('temporal_prec', 0.10, 0.10, value=1.0)
sigmasq_temp = pm.Lambda('sigmasq_temp', lambda temporal_prec=temporal_prec: 1.0/temporal_prec)
# precision matrix associated to time
phi_t = pm.Uniform('phi_t', 0, 10, value=1.0)
diff_time = squareform(pdist(np.arange(1,13).reshape(12,1), metric='minkowski', p=1))
omega = pm.Lambda('omega', lambda diff_time=diff_time, sigmasq_temp=sigmasq_temp, phi_t=phi_t:
inv(sigmasq_temp*np.exp(-phi_t*np.abs(diff_time))))
# calculating value of alpha
V = pm.MvNormal('V', [0]*ntime, omega, value=[0]*ntime)
ones = np.ones((nsite, 1))
# kronecker product of alpha and 1_{nX1} vector
a = pm.Lambda('a', lambda V=V, ones=ones: np.kron(V.reshape(-1,1), ones))

# powered exponential associated to space
diff_space = squareform(pdist(np.array(zip(df['Lon'], df['Lat']))))
spat_prec = pm.Gamma('spat_prec', 0.10, 0.10, value=1.0)
phi = pm.Uniform('phi', 0.0, 10.0, value=1.0)
precision = pm.Lambda('exp', lambda u=diff_space, tau=spat_prec, phi=phi, kappa=1: inv((1/tau)*np.exp(-(phi*u)**kappa)))
# calculating omega from the equation above
w = pm.MvNormal('w', [0]*nsite, precision, value=[0]*nsite)
ones_ = np.ones((ntime, 1))
# kronecker product of 1_{nX1} vector and omega
W = pm.Lambda('W', lambda w=w, ones=ones: np.kron(ones_, w.reshape(-1,1)))

# mean of probability distribution of Y. There is a beta_1 and beta_2 for each time point.
mu = pm.Lambda('mu', lambda beta_1=beta_1, beta_2=beta_2,
a=a, W=W: (np.kron(beta_1.reshape(-1,1), np.ones((nsite,1))) +
np.kron(beta_2.reshape(-1,1), np.ones((nsite,1))) * elevation + a + W).flatten())
# observations
Y = pm.MvNormal('Y', mu, error_prec_m, observed=True, value=temperature)

M = pm.MCMC([Y, mu, error_prec_m, beta_1, beta_2, a, w, phi, phi_t, spat_prec])
M.sample(20000,1000,150)

However, after running this, I encounter errors such as:

Error:

Traceback (most recent call last):
File "/home/user/anaconda/envs/sc/lib/python2.7/site-packages/pymc/database/base.py", line 272, in tally
self._traces[name].tally(chain)
File "/home/user/anaconda/envs/sc/lib/python2.7/site-packages/pymc/database/ram.py", line 99, in tally
self._trace[chain][self._index[chain]] = value
IndexError: index 126 is out of bounds for axis 0 with size 126

%s""" % (name, ''.join(traceback.format_exception(cls, inst, tb))))
/home/user/anaconda/envs/sc/lib/python2.7/site-packages/pymc/database/base.py:282: UserWarning:
Error tallying spat_prec, will not try to tally it again this chain.
Did you make all the samevariables and step methods tallyable
as were tallyable last time you used the database file?

Error:

Traceback (most recent call last):
File "/home/user/anaconda/envs/sc/lib/python2.7/site-packages/pymc/database/base.py", line 272, in tally
self._traces[name].tally(chain)
File "/home/user/anaconda/envs/sc/lib/python2.7/site-packages/pymc/database/ram.py", line 97, in tally
self._trace[chain][self._index[chain]] = value.copy()
IndexError: index 126 is out of bounds for axis 0 with size 126

%s""" % (name, ''.join(traceback.format_exception(cls, inst, tb))))

The errors vary depending on the number of samples I run and the type of sampling (e.g. pm.AdaptiveMetropolis is faster but produces more errors.) Sometimes, I don't get any errors but the results are not satisfactory. Furthermore, M.summary() results in another error:

error_prec_m:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-de6c24f0e4aa> in <module>()
----> 1 M.summary()

/home/user/anaconda/envs/sc/lib/python2.7/site-packages/pymc/Model.pyc in summary(self, variables, alpha, start, batches, chain, roundto)
495             variable.summary(
496                 alpha=alpha, start=start, batches=batches, chain=chain,
--> 497                 roundto=roundto)
498
499     # Property --- status : the sampler state.

/home/user/anaconda/envs/sc/lib/python2.7/site-packages/pymc/Node.pyc in summary(self, alpha, start, batches, chain, roundto)
329         for index in indices:
330             # Extract statistics and convert to string
--> 331             m = str(round(statdict['mean'][index], roundto))
332             sd = str(round(statdict['standard deviation'][index], roundto))
333             mce = str(round(statdict['mc error'][index], roundto))

TypeError: only length-1 arrays can be converted to Python scalars

What is the cause of these errors?

I'm not able to reproduce the first part of your issue, as the model runs just fine. Would it be possible for you to try building from the 2.3 branch?

pip install git+git://github.com/pymc-devs/pymc.git@2.3

If this fixes it for you, then I will try and get a PyMC 2.3.4 build up on Binstar as soon as I can.

I can repeat the second part of your issue, having to do with the summary call, this does appear to be a bug. It is clearly not happy summarizing a 50x50 array. I will fix this.

• Thank you. Sure, I will pip install in a moment and I let you know what happens. – Robert Smith Oct 1 '14 at 3:21
• Unfortunately, I still get the same error in branch 2.3.4. Are there other differences between environments? I have numpy 1.9.0, pandas 0.14.1, scipy 0.14.0 and ipython notebook 2.2.0. – Robert Smith Oct 1 '14 at 3:47
• By the way, did you use a large number of samples? I have seen that the error doesn't appear with something like M.sample(2000,500,150). – Robert Smith Oct 1 '14 at 3:52
• Good news is, I fixed the summary() bug -- thanks for spotting that. I will push it to the 2.3 branch in a few minutes. As for the other error, I still cannot replicate this even by sampling for thousands of iterations. We have similar setups (PyMC does not use Pandas). Can you double check that your build was installed by checking pymc.__version__? – fonnesbeck Oct 1 '14 at 22:06
• Great! I checked the version and it is PyMC 2.3.4. I noticed that I start getting errors when I use M.sample(20000,1000,150). If a remember correctly, even M.sample(15000,1000,150) doesn't result in errors. By the way, I ran this in a clean anaconda environment. Another unusual thing is that variable V outputs an array of zeros only when I pass the parameter value=[0]*ntime as above. – Robert Smith Oct 1 '14 at 22:17