|
|
Contents | Previous | Next | Subchapters |
| Syntax |
[Obj, L, Xhat] = ...meas, function tran, N, X, y) |
| See Also | itrsmo , kalman |
| 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 |
[zk, hk, dhk, Rk] = meas(k, x, y)
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)
|
[gk, dgk, Qk] = tran(k, x, y)
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)
|
g_k(x, y) is defined to be the initial state estimate, and
if k is N,
g_k(x, y) is defined to be zero.
Also note that dgk must be zero for these cases
and that Qk can be any invertible matrix.
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.
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.
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