doc/en/scripts/simple_KalmanFilter2.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 from numpy import array, random
   4 random.seed(1234567)
   5 Xtrue = -0.37727
   6 Yobs = []
   7 for i in range(51):
   8     Yobs.append([random.normal(Xtrue, 0.1, size=(1,)),])
   9 #
  10 print("Estimation of a constant variable by filtering")
  11 print("----------------------------------------------")
  12 print("  Noisy measurements acquired on %i time steps"%(len(Yobs)-1,))
  13 print("")
  14 from adao import adaoBuilder
  15 case = adaoBuilder.New('')
  16 #
  17 case.setObservationOperator(Matrix             = [1.])
  18 case.setObservationError   (ScalarSparseMatrix = 0.1**2)
  19 #
  20 case.setEvolutionModel     (Matrix             = [1.])
  21 case.setEvolutionError     (ScalarSparseMatrix = 1e-5)
  22 #
  23 case.setAlgorithmParameters(
  24     Algorithm="KalmanFilter",
  25     Parameters={
  26         "StoreSupplementaryCalculations":[
  27             "Analysis",
  28             "APosterioriCovariance",
  29             ],
  30         },
  31     )
  32 case.setObserver(
  33     Info="  Analyzed state at current observation:",
  34     Template='ValuePrinter',
  35     Variable='Analysis',
  36     )
  37 #
  38 XaStep, VaStep = 0., 1.
  39 for i in range(1,len(Yobs)):
  40     case.setBackground         (Vector             = "%s"%float(XaStep))
  41     case.setBackgroundError    (ScalarSparseMatrix = "%s"%float(VaStep))
  42     case.setObservation        (Vector             = Yobs[i])
  43     case.execute( nextStep = True )
  44     XaStep = case.get("Analysis")[-1]
  45     VaStep = case.get("APosterioriCovariance")[-1]
  46 #
  47 Xa = case.get("Analysis")
  48 Pa = case.get("APosterioriCovariance")
  49 #
  50 print("")
  51 print("  Final a posteriori variance:",Pa[-1])
  52 print("")
  53 #
  54 #-------------------------------------------------------------------------------
  55 #
  56 Observations = array([yo[0]   for yo in Yobs])
  57 Estimates    = array([xa[0]   for xa in case.get("Analysis")])
  58 Variances    = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
  59 #
  60 import matplotlib.pyplot as plt
  61 plt.rcParams['figure.figsize'] = (10, 4)
  62 #
  63 plt.figure()
  64 plt.plot(Observations,'kx',label='Noisy measurements')
  65 plt.plot(Estimates,'r-',label='Estimated state')
  66 plt.axhline(Xtrue,color='b',label='Truth value')
  67 plt.legend()
  68 plt.title('Estimate of the state', fontweight='bold')
  69 plt.xlabel('Observation step')
  70 plt.ylabel('Voltage')
  71 plt.savefig("simple_KalmanFilter2_state.png")
  72 #
  73 plt.figure()
  74 iobs = range(1,len(Observations))
  75 plt.plot(iobs,Variances[iobs],label='A posteriori error variance')
  76 plt.title('Estimate of the a posteriori error variance', fontweight='bold')
  77 plt.xlabel('Observation step')
  78 plt.ylabel('$(Voltage)^2$')
  79 plt.setp(plt.gca(),'ylim',[0,.01])
  80 plt.savefig("simple_KalmanFilter2_variance.png")