Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Elizabeth Ramirez - Kalman Filters for non-rocket science

Elizabeth Ramirez - Kalman Filters for non-rocket science

Kalman Filters have been widely used for scientific applications. No wonder people often think they involve complex math, however you can actually introduce the Kalman Filter in your daily data processing work, without the complex math you would imagine. This talk will show how to implement the discrete Kalman Filter in Python using NumPy and SciPy.


PyCon 2016

May 29, 2016

More Decks by PyCon 2016

Other Decks in Programming


  1. Background Named after named after Rudolf Kalman Original paper: [

    ] http://www.cs.unc.edu/~welch/kalman/media/pdf/Kalman1960.pdf (http://www.cs.unc.edu/~welch/kalman/media/pdf/Kalman1960.pdf)
  2. Kalman Filters for rocket science Used for Apollo Space Program

    of NASA in early 1960's Transcription of the original code available at [ ] Implemented in AGC4 assembly language CCS: compare and skip TS: transfer to storage CA: clear and add http://www.ibiblio.org/apollo/ (http://www.ibiblio.org/apollo/)
  3. Kalman Filters for non-rocket science Used for some type of

    forecasting problems Generalization of least squares model Time series with varying mean and additive noise
  4. Least Squares Linear system If is square: But if is

    not square: System is overdetermined . Example: 100 points that t
  5. Solution : Find best estimate for state that minimizes: Solve

    for (an estimate) to minimize E Normal Equation :
  6. The Kalman Filter Time varying least squares problem: Estimate at

    each time 1. Recursion 2. Linear combination: innovation : 3. Reliability: covariance matrix
  7. Implementation Output Predicted mean and covariance of the state (before

    the measurement) Estimated mean and covariance of the state (after the measurement) Innovation Filter gain
  8. In [29]: def predict(u, P, F, Q): u = numpy.dot(F,

    u) P = numpy.dot(F, numpy.dot(P, F.T)) + Q return u, P
  9. Correction Predicted state vector Matrix in observation equations Vector of

    observations Predicted covariance matrix Process noise matrix Observations noise matrix
  10. In [30]: def correct(u, A, b, P, Q, R): C

    = numpy.dot(A, numpy.dot(P, A.T)) + R K = numpy.dot(P, numpy.dot(A.T, numpy.linalg.inv(C))) u = u + numpy.dot(K, (b - numpy.dot(A, u))) P = P - numpy.dot(K, numpy.dot(C, K.T)) return u, P
  11. In [124]: dt = 0.1 A = numpy.array([[1, 0], [0,

    1]]) u = numpy.zeros((2, 1)) # Random initial measurement centered at state value b = numpy.array([[u[0, 0] + randn(1)[0]], [u[1, 0] + randn(1)[0]]]) P = numpy.diag((0.01, 0.01)) F = numpy.array([[1.0, dt], [0.0, 1.0]]) # Unit variance for the sake of simplicity Q = numpy.eye(u.shape[0]) R = numpy.eye(b.shape[0])
  12. In [125]: N = 100 predictions, corrections, measurements = [],

    [], [] for k in numpy.arange(0, N): u, P = predict(u, P, F, Q) predictions.append(u) u, P = correct(u, A, b, P, Q, R) corrections.append(u) measurements.append(b) b = numpy.array([[u[0, 0] + randn(1)[0]], [u[1, 0] + randn(1)[0]]]) print 'predicted final estimate: %f' % predictions[-1][0] print 'corrected final estimate: %f' % corrections[-1][0] print 'measured state: %f' % measurements[-1][0] predicted final estimate: -23.417806 corrected final estimate: -22.995292 measured state: -22.720059
  13. In [126]: t = numpy.arange(50, 100) fig = plt.figure(figsize=(15,15)) axes

    = fig.add_subplot(2, 2, 1) axes.set_title("Kalman Filter") axes.plot(t, numpy.array(predictions)[50:100, 0], 'o', label='prediction') axes.plot(t, numpy.array(corrections)[50:100, 0], 'x', label='correction') axes.plot(t, numpy.array(measurements)[50:100, 0], '^', label='measurement') plt.legend() plt.show()
  14. Conclusions Kalman Filter is a viable forecasting technique for time

    series Computationally more ef cient Strang, Gilbert. Computational Science and Engineering