Let's admit it, bayesian modeling on time series is slow. In pymc3, it typically implies using theano scan operation. Here, we will show how to profile one step of the kalman filter, as well as the scan operation over the time series.
First, load the required packages:
In [1]:
import numpy as np
import theano
import theano.tensor as tt
import kalman
We will use the same data as in the 01_RandomWalkPlusObservation notebook.
In [2]:
# True values
T = 500 # Time steps
sigma2_eps0 = 3 # Variance of the observation noise
sigma2_eta0 = 10 # Variance in the update of the mean
# Simulate data
np.random.seed(12345)
eps = np.random.normal(scale=sigma2_eps0**0.5, size=T)
eta = np.random.normal(scale=sigma2_eta0**0.5, size=T)
mu = np.cumsum(eta)
y = mu + eps
Next, we create all the tensors required to describe our model:
In [3]:
# Upon using pymc3, the following theano configuration flag is changed,
# leading to tensors being required to have test values
#theano.config.compute_test_value = 'ignore'
# Tensors for the measurement equation
Z = tt.dmatrix(name='Z')
d = tt.dvector(name='d')
H = tt.dmatrix(name='H')
# Tensors for the transition equation
T = tt.dmatrix(name='T')
c = tt.dvector(name='c')
R = tt.dmatrix(name='R')
Q = tt.dmatrix(name='Q')
# Initial position and uncertainty
a0 = tt.dvector(name='a0')
P0 = tt.dmatrix(name='P0')
We will also create some actual values for them:
In [4]:
ɛ_σ2 = 3.
η_σ2 = 10.
args = dict(Z = np.array([[1.]]),
d = np.array([0.]),
H = np.array([[ɛ_σ2]]),
T = np.array([[1.]]),
c = np.array([0.]),
R = np.array([[1.]]),
Q = np.array([[η_σ2]]),
a0 = np.array([0.]),
P0 = np.array([[1e6]]))
Let's calculate the likelihood of the observed values, given the parameters above:
In [5]:
kalmanTheano = kalman.KalmanTheano(Z, d, H, T, c, R, Q, a0, P0)
(at, Pt, lliks), updates = kalmanTheano.filter(y[:,None])
f = theano.function([Z, d, H, T, c, R, Q, a0, P0], lliks)
llik = f(**args)
llik[1:].sum()
Out[5]:
Time required for the log-likelihood calculation:
In [6]:
print('Measuring time...')
%timeit f(**args)
Profiling a non-scan operation is relatively simple. As an example, let's create a function to calculate the first time step of the Kalman filter:
In [7]:
Y0 = tt.dvector(name='Y0')
_,_,llik = kalman.core._oneStep(Y0, Z, d, H, T, c, R, Q, a0, P0)
profiler = theano.compile.ScanProfileStats()
f = theano.function([Y0, Z, d, H, T, c, R, Q, a0, P0], llik, profile=profiler)
f(y[0,None], **args);
profiler.summary()
Repeating the procedure with a scan procedure, we can see that the code inside it is not profiled. It took me a while to make it work (not even stackoverflow helped!!!). In the end, this is how I made it work:
In [8]:
profiler = theano.compile.ScanProfileStats()
(_,_,llik),_ = kalmanTheano.filter(y[:,None], profile=profiler)
f = theano.function([Z, d, H, T, c, R, Q, a0, P0], llik, profile=profiler)
f(**args);
# Select the node corresponding to the scan operation
scan_op = next(k for k in profiler.op_nodes()
if isinstance(k, theano.scan_module.scan_op.Scan))
scan_op.profile.summary()