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:
\begin{equation} \boldsymbol{Y} = \boldsymbol{\mu} + \boldsymbol{\alpha} \otimes \boldsymbol{1}_{n\times 1} + \boldsymbol{1}_{T\times 1} \otimes \boldsymbol{w} + \boldsymbol{\epsilon} \end{equation}
\begin{equation} 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}) \end{equation}
$\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}$:
\begin{equation} \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} \end{equation}
$\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:
\begin{equation} \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)\end{equation}
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
url = 'http://www.biostat.umn.edu/~brad/data/ColoradoS-T.dat'
df = pd.read_csv(url, sep='\s+')
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?
