Figure 5.16:

The change in 95% confidence intervals for \hat {x}(k\vert k) versus time for a stable, optimal estimator. We start at k=0 with a large initial variance P(0), which gives a large confidence interval.

Code for Figure 5.16

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
%
% Show 95 percent confidence ellipses versus time
% jbr
% 3/30/2013

ntsteps = 10;
U=0.5*[sqrt(2), sqrt(2); sqrt(2), -sqrt(2)];

theta = (-25)*360/(2*pi);
A = [cos(theta), sin(theta); -sin(theta), cos(theta)];
n = size(A, 1);
C = [0.5 0.25];
G = eye(2);
m = 2;
B = eye(m);
g = size(G, 2);
p = size(C, 1);
alpha = 0.95;
nell = 281;
bn = chi2inv(alpha, n);

prior = 3;
% Q0 = prior*U*diag([3, 1/2])*U';
Q0 = prior*eye(n);
Q  = 0.01*eye(g);
R  = 0.01*eye(p);

% steady-state solution
%[Ls, Pmins, Ps, Es] = dlqe(A, G, C, Q, R);

% transient solution
Pminus(1) = {Q0};

for i = 1:ntsteps
  L{i} = Pminus{i}*C'*inv(R+C*Pminus{i}*C');
  P{i} = Pminus{i} - L{i}*C*Pminus{i};
  [exm{i}(:,1), exm{i}(:,2)] = ellipse(inv(Pminus{i}),bn, nell);
  [ex{i}(:,1), ex{i}(:,2)] = ellipse(inv(P{i}),bn, nell);
  if (i == ntsteps)
    break
  end
  Pminus{i+1} = A*P{i}*A' + G*Q*G';
end

figure()
hold on
for i = 1:ntsteps
  plot (ex{i}(:,1),ex{i}(:,2))
end
axis ([-1,1,-2,2]);

figure()
hold on
for i = 1:ntsteps
  plot (exm{i}(:,1),exm{i}(:,2))
end
axis();

exx = cell2mat(ex);
save KFpayoff.dat  exx