|
4 | 4 | "cell_type": "markdown",
|
5 | 5 | "metadata": {},
|
6 | 6 | "source": [
|
7 |
| - "### Derivation of the Kalman Filter\n", |
| 7 | + "# The Linear Dynamical System and the Kalman Filter\n", |
8 | 8 | "\n",
|
| 9 | + "The Kalman filter is a recursive algorithm for estimating the latent state of a linear dynamical system given the observations.\n", |
9 | 10 | "\n",
|
10 |
| - "Implement the Kalman filter with parameters $\\mu, P, A, C, Q, R$. \n", |
| 11 | + "The linear dynamical system is a key model and understanding it as a generative model is key in understanding the principles and limitations of the Kalman filter.\n", |
| 12 | + "\n" |
| 13 | + ] |
| 14 | + }, |
| 15 | + { |
| 16 | + "cell_type": "markdown", |
| 17 | + "metadata": {}, |
| 18 | + "source": [ |
| 19 | + "## The Linear Dynamical System (LDS)\n", |
| 20 | + "\n", |
| 21 | + "A (discrete time) linear dynamical system describes the evolution of the state of a system and\n", |
| 22 | + "the observations that can be obtained from the state.\n", |
| 23 | + "\n", |
| 24 | + "The system is described at times $t=1,2,\\dots$. At each time $t$, the state is denoted by an $N \\times 1$ state vector $x_t$. We don't directly observe the states but obtain observations that are related to the states. The observations are denoted by $M \\times 1$ observation vector $y_t$. The mathematical model is \n", |
11 | 25 | "\n",
|
12 | 26 | "\\begin{eqnarray}\n",
|
13 | 27 | "x_0 & \\sim & {\\mathcal N}(\\mu, P) \\\\\n",
|
14 |
| - "x_{t} & \\sim & {\\mathcal N}(A x_{t-1}, Q) \\\\\n", |
15 |
| - "y_{t} & \\sim & {\\mathcal N}(C x_{t}, R) \n", |
| 28 | + "x_{t} & = & A x_{t-1} + \\epsilon_t \\\\\n", |
| 29 | + "y_{t} & = & C x_{t} + \\nu_t \n", |
16 | 30 | "\\end{eqnarray}\n",
|
17 | 31 | "\n",
|
| 32 | + "here, $\\epsilon_t$ and $\\nu_t$ are assumed to be zero mean Gaussian random variables with variance $Q$ and $P$. \n", |
| 33 | + "\n", |
| 34 | + "Equivalently, we can express the system as a hierarchical generative model\n", |
| 35 | + "\\begin{eqnarray}\n", |
| 36 | + "x_0 & \\sim & {\\mathcal N}(\\mu, P) \\\\\n", |
| 37 | + "x_{t}|x_{t-1} & \\sim & {\\mathcal N}(A x_{t-1}, Q) \\\\\n", |
| 38 | + "y_{t}|x_{t} & \\sim & {\\mathcal N}(C x_{t}, R) \n", |
| 39 | + "\\end{eqnarray}\n", |
18 | 40 | "for $t=1\\dots T$.\n",
|
19 |
| - "The known Parameters:\n", |
20 | 41 | "\n",
|
21 |
| - "$\\mu$ is $N \\times 1$ \n", |
| 42 | + "Both formulations are equivalent; where the former is more often used in engineering and the latter in statistics. Later, we will describe a few small extensions to the model.\n", |
| 43 | + "\n", |
| 44 | + "We can try to understand the meaning of each parameter to understand qualitatively what the LDS is doing. The parameters reflect our knowledge about \n", |
22 | 45 | "\n",
|
23 |
| - "$P$ is $N \\times N$ and diagonal, positive semidefinite\n", |
| 46 | + "* the initial state: $\\mu, P$\n", |
| 47 | + "* the inner working of the system dynamics, the state transition model: $A, Q$\n", |
| 48 | + "* how the states are observed, the observation model $C, R$ \n", |
24 | 49 | "\n",
|
25 |
| - "$A$ is $N \\times N$\n", |
| 50 | + "The initial state mean $\\mu$ is a $N \\times 1$ vector, and the initial state covariance $P$ is an $N \\times N$ positive semidefinite matrix.\n", |
| 51 | + "These two parameters reflect our prior knowledge about the initial state of the process: the expected state and the uncertainty around it as characterized by the covariance matrix.\n", |
26 | 52 | "\n",
|
27 |
| - "$C$ is $M \\times N$\n", |
| 53 | + "Think of the pair $(\\mu, \\Sigma)$ defining an ellipse and the initial state $x_0$ is most likely located somewhere in this ellipse.\n", |
| 54 | + "\n", |
| 55 | + "The state transition matrix $A$ is $N \\times N$, and the state transition noise covariance matrix $Q$ is $N \\times N$, diagonal, positive semidefinite. The state transition matrix $A$ describes the dynamic behaviour of the linear dynamic system: the system would simply undergo a linear transformation if there was no uncertainty about the dynamics\n", |
| 56 | + "$$\n", |
| 57 | + "x_t = A x_{t-1}\n", |
| 58 | + "$$\n", |
| 59 | + "The linear dynamic system model also includes an additive noise term that enables us to model random deviations from an exactly deterministic linear transformation. The deviations are assumed to be zero mean and have the covariance matrix $Q$.\n", |
28 | 60 | "\n",
|
29 |
| - "$Q$ is $N \\times N$ and diagonal, positive semidefinite\n", |
| 61 | + "Given the previous state $x_{t-1}$ at time $t-1$, think of the pair $(A x_{t-1}, Q)$ defining an ellipse and the current state $x_t$ is most likely located somewhere in this ellipse.\n", |
30 | 62 | "\n",
|
31 |
| - "$R$ is $M \\times M$ and diagonal, positive semidefinite\n", |
| 63 | + "Finally, the observation matrix \n", |
32 | 64 | "\n",
|
| 65 | + "$C$ is $M \\times N$, observation matrix\n", |
| 66 | + "$R$ is $M \\times M$ and diagonal, positive semidefinite, observation noise covariance" |
| 67 | + ] |
| 68 | + }, |
| 69 | + { |
| 70 | + "cell_type": "markdown", |
| 71 | + "metadata": {}, |
| 72 | + "source": [ |
| 73 | + "Example: A point object moving with (almost) constant velocity.\n", |
| 74 | + "\n", |
| 75 | + "\n", |
| 76 | + "\n" |
| 77 | + ] |
| 78 | + }, |
| 79 | + { |
| 80 | + "cell_type": "markdown", |
| 81 | + "metadata": {}, |
| 82 | + "source": [ |
33 | 83 | "The input to the algorithm are the parameters and a sequence of \n",
|
34 | 84 | "observations $y_t$ for $t=1 \\dots T$\n",
|
35 | 85 | "\n",
|
|
49 | 99 | "\n",
|
50 | 100 | "-- Covariances: $\\Sigma_{t|t}$\n",
|
51 | 101 | "\n",
|
52 |
| - "-- A sequence of loglikelihoods $l_k$\n", |
| 102 | + "-- A sequence of loglikelihoods $l_k = p(y_k| y_{1:k-1})$\n", |
| 103 | + "\n", |
53 | 104 | "\n",
|
54 |
| - "## Kalman Filter\n", |
| 105 | + "## The algorithm\n", |
55 | 106 | "\n",
|
56 | 107 | "The Kalman filtering algorithm is as follows:\n",
|
57 | 108 | "\n",
|
|
86 | 137 | "End For"
|
87 | 138 | ]
|
88 | 139 | },
|
| 140 | + { |
| 141 | + "cell_type": "code", |
| 142 | + "execution_count": 24, |
| 143 | + "metadata": { |
| 144 | + "collapsed": false, |
| 145 | + "scrolled": false |
| 146 | + }, |
| 147 | + "outputs": [ |
| 148 | + { |
| 149 | + "name": "stdout", |
| 150 | + "output_type": "stream", |
| 151 | + "text": [ |
| 152 | + "mu=\n", |
| 153 | + " [[ 10.]\n", |
| 154 | + " [ 10.]]\n", |
| 155 | + "P=\n", |
| 156 | + " [[ 100. 0.]\n", |
| 157 | + " [ 0. 100.]]\n", |
| 158 | + "A=\n", |
| 159 | + " [[ 12. 4.]\n", |
| 160 | + " [ 1. -3.]]\n", |
| 161 | + "C=\n", |
| 162 | + " [[-3. 5.]\n", |
| 163 | + " [-4. 2.]\n", |
| 164 | + " [ 4. -6.]]\n", |
| 165 | + "Q=\n", |
| 166 | + " [[ 0.1 0. ]\n", |
| 167 | + " [ 0. 0.1]]\n", |
| 168 | + "R=\n", |
| 169 | + " [[ 2. 0. 0.]\n", |
| 170 | + " [ 0. 2. 0.]\n", |
| 171 | + " [ 0. 0. 2.]]\n", |
| 172 | + "observations=\n", |
| 173 | + " [[-1. -5. 6.]\n", |
| 174 | + " [ 3. -0. -5.]\n", |
| 175 | + " [ 1. -1. -8.]]\n" |
| 176 | + ] |
| 177 | + } |
| 178 | + ], |
| 179 | + "source": [ |
| 180 | + "import numpy as np\n", |
| 181 | + "\n", |
| 182 | + "N = 2\n", |
| 183 | + "M = 3\n", |
| 184 | + "\n", |
| 185 | + "A = np.matrix(np.ceil(5*np.random.randn(N,N)))\n", |
| 186 | + "C = np.matrix(np.ceil(5*np.random.randn(M,N)))\n", |
| 187 | + "R = np.matrix(2*np.eye(M))\n", |
| 188 | + "Q = np.matrix(0.1*np.eye(N))\n", |
| 189 | + "\n", |
| 190 | + "mu = np.matrix(10*np.ones((N,1)))\n", |
| 191 | + "P = np.matrix(100*np.eye(N))\n", |
| 192 | + "\n", |
| 193 | + "print('mu=\\n',mu)\n", |
| 194 | + "print('P=\\n',P)\n", |
| 195 | + "\n", |
| 196 | + "print('A=\\n',A)\n", |
| 197 | + "print('C=\\n',C)\n", |
| 198 | + "print('Q=\\n',Q)\n", |
| 199 | + "print('R=\\n',R)\n", |
| 200 | + "\n", |
| 201 | + "T = 3;\n", |
| 202 | + "\n", |
| 203 | + "y = np.matrix(np.ceil(5*np.random.randn(M,T)))\n", |
| 204 | + "\n", |
| 205 | + "print('observations=\\n',y)" |
| 206 | + ] |
| 207 | + }, |
89 | 208 | {
|
90 | 209 | "cell_type": "markdown",
|
91 | 210 | "metadata": {},
|
|
200 | 319 | "All the parameters will be specified when constructing the object and observations will be provided to the filter as a $M \\times T$ matrix."
|
201 | 320 | ]
|
202 | 321 | },
|
203 |
| - { |
204 |
| - "cell_type": "code", |
205 |
| - "execution_count": 24, |
206 |
| - "metadata": { |
207 |
| - "collapsed": false, |
208 |
| - "scrolled": false |
209 |
| - }, |
210 |
| - "outputs": [ |
211 |
| - { |
212 |
| - "name": "stdout", |
213 |
| - "output_type": "stream", |
214 |
| - "text": [ |
215 |
| - "mu=\n", |
216 |
| - " [[ 10.]\n", |
217 |
| - " [ 10.]]\n", |
218 |
| - "P=\n", |
219 |
| - " [[ 100. 0.]\n", |
220 |
| - " [ 0. 100.]]\n", |
221 |
| - "A=\n", |
222 |
| - " [[ 12. 4.]\n", |
223 |
| - " [ 1. -3.]]\n", |
224 |
| - "C=\n", |
225 |
| - " [[-3. 5.]\n", |
226 |
| - " [-4. 2.]\n", |
227 |
| - " [ 4. -6.]]\n", |
228 |
| - "Q=\n", |
229 |
| - " [[ 0.1 0. ]\n", |
230 |
| - " [ 0. 0.1]]\n", |
231 |
| - "R=\n", |
232 |
| - " [[ 2. 0. 0.]\n", |
233 |
| - " [ 0. 2. 0.]\n", |
234 |
| - " [ 0. 0. 2.]]\n", |
235 |
| - "observations=\n", |
236 |
| - " [[-1. -5. 6.]\n", |
237 |
| - " [ 3. -0. -5.]\n", |
238 |
| - " [ 1. -1. -8.]]\n" |
239 |
| - ] |
240 |
| - } |
241 |
| - ], |
242 |
| - "source": [ |
243 |
| - "import numpy as np\n", |
244 |
| - "\n", |
245 |
| - "N = 2\n", |
246 |
| - "M = 3\n", |
247 |
| - "\n", |
248 |
| - "A = np.matrix(np.ceil(5*np.random.randn(N,N)))\n", |
249 |
| - "C = np.matrix(np.ceil(5*np.random.randn(M,N)))\n", |
250 |
| - "R = np.matrix(2*np.eye(M))\n", |
251 |
| - "Q = np.matrix(0.1*np.eye(N))\n", |
252 |
| - "\n", |
253 |
| - "mu = np.matrix(10*np.ones((N,1)))\n", |
254 |
| - "P = np.matrix(100*np.eye(N))\n", |
255 |
| - "\n", |
256 |
| - "print('mu=\\n',mu)\n", |
257 |
| - "print('P=\\n',P)\n", |
258 |
| - "\n", |
259 |
| - "print('A=\\n',A)\n", |
260 |
| - "print('C=\\n',C)\n", |
261 |
| - "print('Q=\\n',Q)\n", |
262 |
| - "print('R=\\n',R)\n", |
263 |
| - "\n", |
264 |
| - "T = 3;\n", |
265 |
| - "\n", |
266 |
| - "y = np.matrix(np.ceil(5*np.random.randn(M,T)))\n", |
267 |
| - "\n", |
268 |
| - "print('observations=\\n',y)" |
269 |
| - ] |
270 |
| - }, |
271 | 322 | {
|
272 | 323 | "cell_type": "code",
|
273 | 324 | "execution_count": 18,
|
|
0 commit comments