Contents Previous Next Subchapters

The Affine Kalman Bucy Smoother and Marginal Likelihood
Syntax [ObjLXhat] = ...
     AffineKbs(function 
meas, function tranNXy)
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(xy) 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(xy) 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(xy)
dhk (mn) d/dx [h_k(xy)]
Rk (mm) 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(xy)
dgk (nn) d/dx [g_k(xy)]
Qk (nn) 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