doc/en/scripts/simple_KalmanFilter1.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.setBackground         (Vector             = [0.])
  18 case.setBackgroundError    (ScalarSparseMatrix = 1.)
  19 #
  20 case.setObservationOperator(Matrix             = [1.])
  21 case.setObservation        (VectorSerie        = Yobs)
  22 case.setObservationError   (ScalarSparseMatrix = 0.1**2)
  23 #
  24 case.setEvolutionModel     (Matrix             = [1.])
  25 case.setEvolutionError     (ScalarSparseMatrix = 1e-5)
  26 #
  27 case.setAlgorithmParameters(
  28     Algorithm="KalmanFilter",
  29     Parameters={
  30         "StoreSupplementaryCalculations":[
  31             "Analysis",
  32             "APosterioriCovariance",
  33             ],
  34         },
  35     )
  36 case.setObserver(
  37     Info="  Analyzed state at current observation:",
  38     Template='ValuePrinter',
  39     Variable='Analysis',
  40     )
  41 #
  42 case.execute()
  43 Xa = case.get("Analysis")
  44 Pa = case.get("APosterioriCovariance")
  45 #
  46 print("")
  47 print("  Final a posteriori variance:", Pa[-1])
  48 print("")
  49 #
  50 #-------------------------------------------------------------------------------
  51 #
  52 Observations = array([yo[0]   for yo in Yobs])
  53 Estimates    = array([xa[0]   for xa in case.get("Analysis")])
  54 Variances    = array([pa[0,0] for pa in case.get("APosterioriCovariance")])
  55 #
  56 import matplotlib.pyplot as plt
  57 plt.rcParams['figure.figsize'] = (10, 4)
  58 #
  59 plt.figure()
  60 plt.plot(Observations,'kx',label='Noisy measurements')
  61 plt.plot(Estimates,'r-',label='Estimated state')
  62 plt.axhline(Xtrue,color='b',label='Truth value')
  63 plt.legend()
  64 plt.title('Estimate of the state', fontweight='bold')
  65 plt.xlabel('Observation step')
  66 plt.ylabel('Voltage')
  67 plt.savefig("simple_KalmanFilter1_state.png")
  68 #
  69 plt.figure()
  70 iobs = range(1,len(Observations))
  71 plt.plot(iobs,Variances[iobs],label='A posteriori error variance')
  72 plt.title('Estimate of the a posteriori error variance', fontweight='bold')
  73 plt.xlabel('Observation step')
  74 plt.ylabel('$(Voltage)^2$')
  75 plt.setp(plt.gca(),'ylim',[0,.01])
  76 plt.savefig("simple_KalmanFilter1_variance.png")