Solved-Kalman Filtering (part 1) –Solution

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

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.

2. 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 .

2. 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 .

4. 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 _XA^T +A _X;YB^T +B _Y;XA^T +B _YB^T

Similarly, nd expressions for _Z ; _W;Z ; _Z in terms of _X ; _X;Y ; _Y .

5. Let _1M 2 R^M denote the M 1 column vector, whose elements are all the number 1. Suppose A 2 R^{n M} . Let a_k 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

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 R^{n n} is invertible, then W ^T W 0″, try begin code

• W = rand(5,5);

• eig(W’*W)

end code

3 of 8

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

8. Two important commands are rand and randn.

8. 1. Use hist to visualize the distribution of the numbers generated by A = rand(1,500).

8. 2. 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,

13. 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 _M¹ .

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 : ! Rⁿ.

Using Xmat = rand(n, M);, the mean of X will be approximately 0:5 _1n, and the variance will be approximately ₁₂¹ I_n.

Using Xmat = randn(n, M);, the mean of X will be approximately 0 _1n, and the variance will be approximately I_n.

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 I_n.

4 of 8

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.

5. 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

6. 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 = AA^T . 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

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.

9. 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 n_x = 5; n_d = 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(x₀) 0_{5 1}; E(d_k) 0_{3 1}; E(d_kd^T_k ) I₃; E(d_kd^T_j ) 0_{3 3}

6 of 8

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 d₀; d₁; : : : ; d₉₉ 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(x₀) 0_{5 1};

and for all k and all j 6= k

E(d_k) 0_{3 1}; E(d_kd^T_k ) I₃; E(d_kd^T_j ) 0_{3 3}

(a) We wish to simulate the system

x_k+1 = Ax_k + Ed_k; y_k = Cx_k + F d_k

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

Write code that (via for loop) picks 8 outcomes from a probability model (for x₀ 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, y_k^[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 x₀ and d.