$29
Assignment: Kalman Filtering (part 1) Solution
Problem set due Tuesday, February 27, 8:00PM. You must submit a pdf le (from a typed document, or a neat scan of your handwritten/typed solution) via bCourses. The copier on 6th oor of Etcheverry is available 24/7, and can be used without a login to email scans to any berkeley.edu address.
Histograms in Matlab: If vecD is a row (or column) vector, the command hist(vecD) draws a histogram of the numerical values in vecD. Recall what a histogram is from high-school statistics. The command determines the minimum and maximum values contained within the data; splits that interval into several (default is 10) equally spaced bins; counts the number of data points within each bin; and plots the counts (vertical axis) versus the bin-intervals (horizontal) axis. Such a gure shows (in a coarse, but informative manner) how the numerical values within the data are distributed. Obviously the order of the numbers with the vector vecD does not a ect the plot, since the plot represents counts of values within each bin.
1
2
3
4
5
6
7
1
2
3
4
5
6
The command linspace generates points uniformly spaced between speci ed limits. Histograms of such data should be at, since there will be a uniform number of points within intervals of the same width. Using randperm (which creates a random permutation of the indices), the order of the data can be rearranged, but the histogram remains the same.
begin code
N = 500;
X = linspace(-3,2,N);
subplot(2,1,1); hist(X);
I = randperm(N);
Y = X(I);
isequal(X,Y)
subplot(2,1,2); hist(Y);
end code
For a collection of angles uniformly chosen around (and close to) 0, the cosine of the angles are more often close to 1 than to smaller values. A histogram makes this evident
begin code
N = 500;
X = linspace(-1,1,N);
subplot(2,1,1); hist(X);
xlabel(’Angle, rads’); ylabel(’Counts’)
subplot(2,1,2); hist(cos(X));
xlabel(’COS(Angle, rads)’); ylabel(’Counts’)
end code
The command logspace generates points logarithmically spaced between speci ed limits. His-tograms of such data should show a nonuniform distribution, since the horizontal axes in a his-
1 of 8
Available: Feb 20 Due: Feb 27
1
2
3
4
5
6
7
1
2
togram is (by default) linearly spaced. Creating 500 points logarithmically spaced between 0.1 and 10 will give 250 points from 0.1 to 1, and 250 points from 1 to 10. Hence there will be many more points in the bin [0 1] than there are in the bin [9 10]. Also, again, using randperm, the order of the data can be rearranged, but the histogram remains the same.
begin code
N = 500;
X = logspace(-1,1,N);
subplot(2,1,1); hist(X);
I = randperm(N);
Y = X(I);
isequal(X,Y)
subplot(2,1,2); hist(Y);
end code
A ner (or coarser) bin selection can be done with additional arguments. The simplest syntax is to alter the number of bins.
begin code
nBins = 25;
subplot(2,1,1); hist(X,nBins);
end code
Experiment more with histogram and be comfortable with its behavior to examine the distribution of scalar-valued data.
Suppose X is a random variable de ned on a probability model consisting of a sample space and probability law P. Let a; b 2 R be real numbers. De ne a new random variable Y := aX + b. In other words, for each ! 2 , de ne Y (!) := aX(!) + b. Find Y ; Y Y and XY in terms of a; b; X and XX .
Suppose X is a random variable de ned on a probability model consisting of a sample space and probability law P. Let a; b 2 R be real numbers. De ne a new random variable Y :=
a(X X ) + b. In other words, for each ! 2 , de ne Y (!) := a (X(!) X ) + b. Find Y ; Y Y and XY in terms of a; b; X and XX .
Suppose X and Y are vector-valued random variables, and matrices A; B; C; D are matrices of appropriate dimension. De ne
W :=AX+BY; Z :=CX+DY
Verify the expressions below for W and W in terms of X ; Y ; X ; X;Y ; Y ,
W =A X +B Y; W =A XAT +A X;YBT +B Y;XAT +B YBT
Similarly, nd expressions for Z ; W;Z ; Z in terms of X ; X;Y ; Y .
Let 1M 2 RM denote the M 1 column vector, whose elements are all the number 1. Suppose A 2 Rn M . Let ak denote the k’th column of A. Show that
1
M
1
1
M
1
X
ak =
A 1M;
X
akakT =
AAT
M
k=1
M
M
k=1
M
2 of 8
Available: Feb 20 Due: Feb 27
For an array A, how are the two assignments below related to the above expressions? begin code
B = (1/size(A,2))*A*ones(size(A,2),1)
2 C = (1/size(A,2))*A*A’
end code
If A 2 Rn M , then mean(A,1) computes the average of each column (ie., \down" the rst dimension), returning a 1 M array, while mean(A,2) computes the average of each row (ie., \across" the 2nd dimension), returning a n 1.
For an array A, what do the assignments below accomplish? begin code
B = A - repmat(mean(A,2),[1 size(A,2)])
C = A - repmat(mean(A,1),[size(A,1) 1]) end code
There are several commands in Matlab that produce random numbers. Without needing to know anything about the properties of the numbers that are produced, these are useful commands to quickly generate matrices populated with arbitrary numbers, and can be very useful for experi-mentation of calculations and learning behavior of basic Matlab commands.
For example, to test the claim that real, symmetric matrices have real eigenvalues, but real nonsymmetric matrices generally have complex eigenvalues, you can try this experiment using rand
1
2
3
4
5
6
begin code
A = rand(7,7);
isequal(A, A’)
eig(A)
Asym = A + A’; % make symmetric by adding transpose
isequal(Asym, Asym’)
eig(Asym)
end code
To test the imprecise claim that \square, randomly generated matrices are almost always invertible", try
1
2
3
begin code
A = rand(5,5);
shouldBeIdentity = A*inv(A)
norm(shouldBeIdentity-eye(5))
end code
To test the claim that \if W 2 Rn n is invertible, then W T W 0", try begin code
W = rand(5,5);
eig(W’*W)
end code
3 of 8
Available: Feb 20 Due: Feb 27
1
2
3
4
5
6
7
To test the claim that \sqrtm computes the unique positive-semide nite symmetric matrix square root of a positive-semide nite matrix," try
begin code
W = rand(5,5);
M = W’*W;
isequal(M, M’)
S = sqrtm(M);
isequal(S, S’)
eig(S)
S*S-M
end code
Two important commands are rand and randn.
Use hist to visualize the distribution of the numbers generated by A = rand(1,500).
The command randn generates random numbers with a di erent distribrution than that provided by rand. Use hist to visualize the distribution of the numbers generated by A = randn(1,500)
The properties of these random number generators are described by the distribution of the random variables that they produce, namely \uniform" for rand and \normal" (or ‘’gaussian") for randn. In this class, we did not discuss the distribution of a random variable (scalar or vector-valued). We only covered three properties: mean, variance and covariance. We talked about probability models with a countable number of outcomes, so these 3 properties were de ned in terms of in nite sums over the sample space. In that limited context, we can interpret the outputs of Xmat = rand(n,
and Xmat = randn(n, M) as follows:
Each column is associated with an outcome of the probability model, ( ; P ). There are M outcomes in this probability model, and the probabilities of each outcome are equal, each being M1 .
The value of a random variable X at the k’th outcome is the k’th column of Xmat. Since the columns are n-dimensional, the random variable X is vector-valued, namely X : ! Rn.
Using Xmat = rand(n, M);, the mean of X will be approximately 0:5 1n, and the variance will be approximately 121 In.
Using Xmat = randn(n, M);, the mean of X will be approximately 0 1n, and the variance will be approximately In.
Hence, in this interpretation, both commands create a probability model and random variables whose components are (approximately) uncorrelated, with (approximately) known means and variances.
(c) Verify the claims above for n = 3, using ideas from problems (5) and (6).
We know that a ne transformations (problems (2) and (3)) change the mean and variance, hence Ymat = sqrt(12)*(rand(n, M)-0.5) corresponds to a random variable Y with mean approximately 0, and the variance approximately In.
4 of 8
Available: Feb 20 Due: Feb 27
Run the code below, and then verify that the means and variances are approximately equal, as claimed
1
2
3
4
begin code
M = 500;
n = 1;
Xmat = sqrt(12)*(rand(n, M)-0.5);
Ymat = randn(n, M);
end code
Recall though, from parts (a) and (b), that the actual distributions are quite di erent.
Run the code below, and then verify that the means and variance are approximately as claimed above
1
2
3
begin code
M = 500;
n = 3;
Xmat = sqrt(12)*(rand(n, M)-0.5);
end code
Recall the results from (4) which show that if A is a matrix, and X is a vector-valued random variable with X = 0 and X = I, then Y := AX will have Y = 0 and Y = AAT . Explain how the following code generates a random variable Y (de ned on a nite Sample Space with 1000 outcomes) such that
Y 031;
Y
2
1
2
0:5 3
4
3
1
1
5
1
0:5
1
1
2
3
4
5
6
Also, verify that the construction of Y actually achieves the goals.
begin code
n = 3;
M = 1000;
X = randn(n, M);
S = [3, -1, 1; -1, 2, 0.5; 1, 0.5, 1];
A = sqrtm(S);
Y = A*X;
end code
(g) Somewhat surprisingly (perhaps), rearranging the numbers from rand and randn still pre-serves the approximate properties discussed above
begin code
n = 4;
M = 1000;
Xmat = randn(n, M);
muX = mean(Xmat,2)
D = Xmat - repmat(muX,[1 M]);
D*D’/M
7 I = randperm(n*M);
5 of 8
Available: Feb 20 Due: Feb 27
8
9
10
11
Ymat = reshape(Xmat(I),[n M]);
muY = mean(Ymat,2)
D = Ymat - repmat(muY,[1 M]);
D*D’/M
end code
Task: Verify the same rearrangment property for data produced by rand. Remark/Connection: In this class, we did not discuss the notions of independent events in a probability model, and independent random variables. These important concepts are why arbitrary rearrangements performed above have little e ect on the means and variances of the randome variables produced by rand and randn.
Generating random variables for initial conditions and disturbances: In this class, we began using random variables in the Kalman ltering module to represent a large collection of \scenarios" of initial conditions and disturbances for which the estimation scheme we design must work well, at least in an average sense. For use in simulation, we want to be able to quickly create a large collection of scenarios with appropriate properties (certain means and variances), and then be able to select one speci c scenario, and simulate the system under those conditions (and possibly repeat this for some of the other scenarios).
In order to create candidate random variables representing initial conditions and disturbances with certain means and variances, we can use rand and randn, and the rearrangment properties, and the a ne transformations to do so easily.
For example, take a system 5 states, and 3 disturbances, so nx = 5; nd = 3. In order to create a probability model of initial conditions and sequences of disturbances of length 100 (ie., for k = 0; 1; : : : ; 99), we can use rand and randn. Suppose we want a probability model with 5000 di erent outcomes, all with means of 0 and variance matrices equal to I.
1
2
3
4
5
begin code
M = 5000; % outcomes in probability model
nX = 5;
nD = 3;
duration = 100; % length of disturbance sequence R = randn(nX+nD*duration, M);
end code
The matrix R has all of the data representing the M(= 5000) outcomes, namely
8
2 d0
3
2 d0
3
2 d0
3
9
x0
x0
x0
d1
d1
d1
7
6
7
6
7
6
.
.
.
.
7
6
.
7
6
.
7
< 6
.
.
.
=
6
7
6
7
6
7
6
d
7
6
d
7
6
d
7
6
99
7 @!1
6
99
7 @!2
6
99
7 @!M
5
4
5
4
5
4
:
;
By construction, we know that over this probability model has for all k and k 6= j
E(x0) 05 1; E(dk) 03 1; E(dkdTk ) I3; E(dkdTj ) 03 3
6 of 8
Available: Feb 20 Due: Feb 27
Each column of R has dimension 305 1, and is an outcome of the random variable
3
x0
d0 7
7
d1 7
7
.
6 . 7
4 . 5
d99
For simulation studies, we can chose any such column, and use that in the simulation.
Because of the properties of randn, we can get the same results (in an average sense) by using randn to generate just the one single column we would use in a simulation. And, we can generate these in di erent rearrangements, perhaps more conducive to our needs, for example the single 305 1 column might be easier to manipulate if it was broken out into the initial condition vector, and the sequence of disturbances stored column-wise, as below.
1
2
begin code
x0vec = randn(nX, 1);
dSeq = randn(nD, duration);
end code
In this construction the x0vec vector is the initial condition and the columns of dSeq are the values of the disturbance d0; d1; : : : ; d99 associated with one outcome.
This data represents a single outcome from an underlying probability model with and by construc-tion, we know that over this probability model
E(x0) 05 1;
and for all k and all j 6= k
E(dk) 03 1; E(dkdTk ) I3; E(dkdTj ) 03 3
(a) We wish to simulate the system
xk+1 = Axk + Edk; yk = Cxk + F dk
under random initial conditions and random disturbances with the following properties
E(x0) =
2 1 3
;
x0 = 2
1
2
0:5 3
;
4
1
5
4
3
1
1
5
1
1
0:5
1
and for all k 0, and j 6= k
0:5
1
E(dk) = 0;x0;dk = 0;dk;dj
= 0;dk
0:5
0:5
with the following state-space data
A = 2
1:5
0:02
1:3 3
;
E =
2 0
0 3
4
0:8
0:6
0:15
1
0
0:26
0:16
1 5
4 0
0 5
7 of 8
Available: Feb 20
Due: Feb 27
C =
0
1
1
; F =
0
1
1
1
0
0
1
Write code that (via for loop) picks 8 outcomes from a probability model (for x0 and d) with those properties, and simulates the system from k = 0 to k = 30. The code should plot (on a single plot, all 8 outcomes) the second component of y, yk[2] versus k. Please document the code so that we can see how you are getting (in an average sense) the appropriate mean and variance for x0 and d.
8 of 8
