doc/fr/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 par filtrage d'une variable constante")
  11 print("------------------------------------------------")
  12 print("  Observations bruitées acquises sur %i pas de temps"%(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.setObservationError   (ScalarSparseMatrix = 0.1**2)
  22 #
  23 case.setEvolutionModel     (Matrix             = [1.])
  24 case.setEvolutionError     (ScalarSparseMatrix = 1e-5)
  25 #
  26 case.setAlgorithmParameters(
  27     Algorithm="KalmanFilter",
  28     Parameters={
  29         "StoreSupplementaryCalculations":[
  30             "Analysis",
  31             "APosterioriCovariance",
  32             ],
  33         },
  34     )
  35 #
  36 # Boucle pour obtenir une analyse à l'arrivée de chaque observation
  37 #
  38 for i in range(1,len(Yobs)):
  39     case.setObservation(Vector = Yobs[i])
  40     case.execute( nextStep = True )
  41 #
  42 Xa = case.get("Analysis")
  43 Pa = case.get("APosterioriCovariance")
  44 #
  45 print("  État analysé à l'observation finale :", Xa[-1])
  46 print("")
  47 print("  Variance a posteriori finale :", 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='Mesures bruitées')
  61 plt.plot(Estimates,'r-',label='État estimé')
  62 plt.axhline(Xtrue,color='b',label='Valeur vraie')
  63 plt.legend()
  64 plt.title('Estimation de l\'état', fontweight='bold')
  65 plt.xlabel('Pas d\'observation')
  66 plt.ylabel('Tension')
  67 plt.savefig("simple_KalmanFilter2_state.png")
  68 #
  69 plt.figure()
  70 iobs = range(1,len(Observations))
  71 plt.plot(iobs,Variances[iobs],label='Variance d\'erreur a posteriori')
  72 plt.title('Estimation de la variance d\'erreur a posteriori', fontweight='bold')
  73 plt.xlabel('Pas d\'observation')
  74 plt.ylabel('$(Tension)^2$')
  75 plt.setp(plt.gca(),'ylim',[0,.01])
  76 plt.savefig("simple_KalmanFilter2_variance.png")