<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Step1:  We will use the same data as in the 01_RandomWalkPlusObservation notebook.
Step2:  Next, we create all the tensors required to describe our model
Step3:  We will also create some actual values for them
Step4:  Let's calculate the likelihood of the observed values, given the parameters above
Step5:  Time required for the log-likelihood calculation
Step6:  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
Step7:  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
<ASSISTANT_TASK:>
Python Code:
import numpy              as np
import theano
import theano.tensor      as tt
import kalman
Explanation: 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:
End of explanation
# 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
Explanation: We will use the same data as in the 01_RandomWalkPlusObservation notebook.
End of explanation
# 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')
Explanation: Next, we create all the tensors required to describe our model:
End of explanation
ɛ_σ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]]))
Explanation: We will also create some actual values for them:
End of explanation
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()
Explanation: Let's calculate the likelihood of the observed values, given the parameters above:
End of explanation
print('Measuring time...')
%timeit f(**args)
Explanation: Time required for the log-likelihood calculation:
End of explanation
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()
Explanation: 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:
End of explanation
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()
Explanation: 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:
End of explanation
<END_TASK>