Contents Previous Next Subchapters Current Chapters-> kalman itrsmo AffineKbs Parent Chapters-> Omatrix6 filtering kalman_bucy AffineKbs Search Tools-> contents reference index search

The Affine Kalman Bucy Smoother and Marginal Likelihood
 Syntax [Obj, L, Xhat] = ...      AffineKbs(function meas, function tran, N, X, y) See Also itrsmo , kalman

Bell:00
The marginal likelihood for parameters in a discrete Gauss-Markov process
, IEEE Transactions On Signal Processing, March 2000.

Description
Calculates the Kalman-Bucy smoother solution and the marginal likelihood for the affine case using the formula in Theorem 8 of [Bell:00] .

Notation
 Math Description n number of components in the state vector at one time point m number of components in the data vector at k-th time index N number of time points z_k value of the data vector at k-th time index R_k(y) variance of v_k; i.e., noise in the measurement equation h_k(x, y) expected value of z_k given x at k-th time index Q_k(y) variance of w_k; i.e., the noise in the transition equation g_k(x, y) expected value of x at next time index given value at k-th time index
Measurement Model
[zkhkdhkRk] = meas(kxy)
The input value k is an integer between 1 and N. The input value x is a the real or double-precision column vector of length n. The input value y is a the real or double-precision column vector equal to the corresponding argument of AffineKbs. The output values are all real or double-precision and are specified by
 Notation Dimension Description zk (m, 1) measurement at k-th time index hk (m, 1) h_k(x, y) dhk (m, n) d/dx [h_k(x, y)] Rk (m, m) R_k(y)

Transition Model
[gkdgkQk] = tran(kxy)
The input value k is an integer between 0 and N. The input value x is a the real or double-precision column vector of length n. The input value y is a the real or double-precision column vector equal to the corresponding argument of AffineKbs. The output values are all real or double-precision and are specified by
 Notation Dimension Description gk (n, 1) g_k(x, y) dgk (n, n) d/dx [g_k(x, y)] Qk (n, n) Q_k(y)
Note that as specified in [Bell:00] : if k is zero, g_k(xy) is defined to be the initial state estimate, and if k is N, g_k(xy) is defined to be zero. Also note that dgk must be zero for these cases and that Qk can be any invertible matrix.

Other Inputs
The integer scalar N specifies the number of time index values. The real or double-precision column vector X is a column vector of length N * n containing a sequence of state values for all time around which the affine approximations of the model functions are made. The real or double-precision vector y specifies the value of the parameter vector.

Outputs
The output value Obj is a real or double-precision scalar containing the affine approximation for the marginal likelihood as specified in Theorem 8 of [Bell:00] . The output value L is a real or double precision scalar containing the value of the joint likelihood corresponding to Xhat and y; i.e., the term L(Xhat(y), y) in Theorem 8 of [Bell:00] . The output value Xhat is a real or double-precision column vector of length N * n containing a sequence of state values for all time. This is the solution to the affine smoother problem corresponding to X and y.

Example
The following example uses AffineKbs to estimate the initial dose, noise levels, and decay rate corresponding to a one compartment stochastic model. The elements of the parameter vector y are
 y(1) initial dose y(2) decay rate y(3) process variance y(4) measurement variance

clear                    # initialize O-Matrix
N     = 50               # number of data points
T     = 1.:1.:N          # grid of time points between 1 and N
Zsim  = fill(0., N, 1)   # vector that will hold simulated data
B     = 1.               # Basal value of concentration
#
# define the transition function
function [gk, dgk, Qk] = tran(k, x, y) begin
global N, T, B                       # make variables visible in function
if k == 0 then begin
Phi = exp(-y(2) * T(1))        # transition matrix for differential equation
gk  = (y(1) - B) * Phi + B     # solution corresponding to initial x = y(1)
dgk = 0.                       # derivative of solution w.r.t x
Qk  = y(3)                     # variance in transition equation
end else if k == N then begin
gk  = 0.                       # values for this case not used
dgk = 0.
Qk  = 1.                       # Qk must be invertible
end else begin
dt  = T(k+1) - T(k)            # time to next measurement
Phi = exp(-y(2) * dt)          # transition matrix for differential equation
gk  = (x - B) * Phi + B        # solution corresponding to previous x value
dgk = Phi                      # derivative of solution w.r.t x
Qk  = y(3)                     # variance in transition equation
end
return
end
#
# define the measurement function
function [zk, hk, dhk, Rk] = meas(k, x, y) begin
global Zsim, T         # make visible in function
zk  = Zsim(k)          # data value for this time index
hk  = x                # model for the data
dhk = 1.               # derivative of the model w.r.t. x
Rk  = y(4)             # variance of the model
return
end
#
# simulation
ranseed                                       # seed the random number generator
ysim = {10., 5./N, .1, .1}                    # parameter values for simulation
[gk, dgk, Qk] = tran(0, novalue, ysim)        # initial state estimate and variance
for k = 1 to N begin
wk         = snormal(1, 1) * sqrt(Qk)   # process noise
xk         = gk + wk                    # next state value for this realization
[zk, hk, dhk, Rk] = meas(k, xk, ysim)   # measurement information
vk         = snormal(1, 1) * sqrt(Rk)   # measurement noise
Zsim(k)    = hk + vk                    # measurement value
[gk, dgk, Qk] = tran(k, xk, ysim)       # transition information
end
#
# arguments to optimizer
X0    = fill(0d0, N, 1)                   # initial estimate of state for all indices
yin   = ysim / 2                          # initial value of parameters
step  = ysim * 1e-3                       # convergence criteria and step size
mitr  = 100                               # maximum number of iterations
level = 1                                 # level of tracing
function f = AffineKbs(meas, tran, N, X0) # objective function
#
# optimization
yout = conjdir(function f, yin, step, mitr, level)
yest = yout.col(coldim(yout))
#
print "Parameter vector used for simulation  =", ysim'
print "Estimated parameter vector            =", yest'
#
[Obj, L, Xhat] = AffineKbs(function meas, function tran, N, X0, yest)
gplot(T, Zsim, "square")    # simulated data values
gplot(T, Xhat, "solid")     # affine smoother estimate for state values