src/daComposant/daNumerics/ComputeKhi2.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2009  EDF R&D
   4 #
   5 #  This library is free software; you can redistribute it and/or
   6 #  modify it under the terms of the GNU Lesser General Public
   7 #  License as published by the Free Software Foundation; either
   8 #  version 2.1 of the License.
   9 #
  10 #  This library is distributed in the hope that it will be useful,
  11 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 #  Lesser General Public License for more details.
  14 #
  15 #  You should have received a copy of the GNU Lesser General Public
  16 #  License along with this library; if not, write to the Free Software
  17 #  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 #  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 __doc__ = """
  22     Outil numerique de calcul de la variable Khi2
  23
  24     On peut realiser deux types de test du Khi2 :
  25         - test d'adequation : comparer la distribution d'un echantillon a une
  26           distribution theorique,
  27         - test d'homogeneite : comparer les distributions de 2 vecteurs.
  28
  29     Pour le test d'adequation, on travaille sur une gaussienne
  30     dont la moyenne et l'ecart type sont calcules sur
  31     l'echantillon, soit donnes.
  32
  33     Ce fichier contient une classe "StatspourTests" de methodes qui realisent
  34     differentes etapes utiles aux calculs des tests du Khi2.
  35
  36     Ce fichier contient de plus 3  methodes : ComputeKhi2_testGauss,
  37     ComputeKhi2_Gauss et ComputeKhi2_Homogen.
  38     - ComputeKhi2_testGauss : calcul la distance du Khi2 entre un vecteur
  39       aleatoire issu d un gaussienne et une distribution theorique gaussienne
  40       dont on specifie la moyenne et l ecart type
  41     - ComputeKhi2_Gauss : calcul la distance du Khi2 entre un vecteur donne et
  42       une distribution theorique gaussienne dont la moyenne et l ecart type sont
  43       calcules sur l echantillon
  44     - ComputeKhi2_Homogen : calcul la distance du Khi2 entre deux vecteurs donnes
  45
  46     Ces methodes necessitent et fournissent :
  47         - en input :
  48             - le ou les vecteurs dont on etudie la distribution,
  49             - la distribution theorique et eventuellement la moyenne et ecart type,
  50             - la largeur des classes,
  51             - un booleen traduisant la suppression des classes vides
  52         - en output :
  53             - le vecteur des classes,
  54             - les pdf theorique et donnee,
  55             - la valeur du Khi2,
  56             - la p-value qui represent l'aire de la queue de la distribution du
  57               Khi2 et
  58             - le message qui interprete le test.
  59 """
  60 __author__ = "Sophie RICCI - Mars 2010"
  61
  62 import numpy
  63 from numpy import random
  64 from scipy import arange, asarray, stats
  65 from scipy.stats import histogram2, chisquare, chisqprob, norm
  66 import logging
  67
  68 # ==============================================================================
  69 class StatspourTests :
  70     """
  71     Classe de methodes pour la preparation du test de Khi2
  72     """
  73     def __init__(self, cdftheo=None, meantheo = None, stdtheo = None, pdftest=None,obs=None,use_mean_std_exp=True, dxmin=0.01, obsHomogen = None, nbclasses = None) :
  74
  75
  76         if (pdftest is None and obs is None) :
  77            raise ValueError('Donner soit une pdf de test soit un vecteur obs')
  78         if not obs is None :
  79             if pdftest is None :
  80               self.__obs=asarray(obs)
  81         if not pdftest is None :
  82             if obs is None :
  83                if len(pdftest) == 3:
  84                   niter=eval(pdftest[2])
  85                   obs=[eval(" ".join(pdftest[:2])) for z in range(niter)]
  86                   self.__obs=asarray(obs)
  87                else :
  88                   self.__obs=asarray(eval(" ".join(pdftest[:2])))
  89         if not (obsHomogen is None) :
  90           self.__obsHomogen = asarray(obsHomogen)
  91           self.__testHomogen =  True
  92         else :
  93           self.__testHomogen =  False
  94
  95
  96         self.__mean_exp = self.__obs.mean()
  97         self.__std_exp = self.__obs.std()
  98
  99         if cdftheo is None : raiseValueError(" ... Definir le parametre cdftheo ...")
 100         if  use_mean_std_exp :
 101           self.__cdf=cdftheo( self.__mean_exp, self.__std_exp).cdf
 102         else :
 103           self.__cdf=cdftheo( meantheo, stdtheo).cdf
 104
 105         self.__min=min(self.__obs)
 106         self.__max=max(self.__obs)
 107         self.__N=len(self.__obs)
 108         self.__use_mean_std_exp=use_mean_std_exp
 109         self.__dxmin=dxmin
 110         self.__nbclasses = nbclasses
 111         if not (dxmin is None) and  not (nbclasses is None) :
 112            raise ValueError("... Specifier soit le nombre de classes, soit la largeur des classes")
 113         if  (dxmin is None) and   (nbclasses is None) :
 114            raise ValueError("... Specifier soit le nombre de classes, soit la largeur des classes")
 115         if not (nbclasses is None) and (dxmin is None) :
 116           self.__dxmin = (self.__max - self.__min ) / float(self.__nbclasses)
 117         return None
 118
 119     def MakeClasses(self) :
 120         """
 121         Classification en classes
 122         """
 123         self.__subdiv=arange(self.__min,self.__max+self.__dxmin,self.__dxmin)
 124         self.__modalites=len(self.__subdiv)
 125         return None
 126
 127     def ComputeObs(self):
 128         """
 129         Calcul de la probabilite observee de chaque classe
 130         """
 131         self.__kobs=histogram2(self.__obs,self.__subdiv)[1:]
 132         return self.__kobs
 133
 134     def ComputeObsHomogen(self):
 135         """
 136         Calcul de la probabilite observee pour le test homogeneite de chaque classe
 137         """
 138         self.__kobsHomogen=histogram2(self.__obsHomogen,self.__subdiv)[1:]
 139         return self.__kobsHomogen
 140
 141     def ComputeTheo(self):
 142         """
 143         Calcul de la probabilite theorique de chaque classe
 144         """
 145         self.__ktheo=[self.__cdf(self.__subdiv[i+1])-self.__cdf(self.__subdiv[i]) for i in range(self.__modalites-1)]
 146         self.__ktheo=asarray(self.__ktheo)
 147         self.__ktheo=(sum(self.__kobs)/sum(self.__ktheo))*self.__ktheo
 148
 149     def Computepdfs(self) :
 150
 151         self.__subdiv=self.__subdiv[1:]
 152         self.__pdfobs=[self.__kobs[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)]
 153
 154         if self.__testHomogen :
 155             self.__pdftheo=[self.__kobsHomogen[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)]
 156         else :
 157             self.__pdftheo=[self.__ktheo[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)]
 158
 159         return self.__subdiv, self.__pdftheo, self.__pdfobs
 160
 161     def Computeddl(self):
 162         """
 163         Calcul du nombre de degres de liberte
 164         """
 165         self.__ddl = self.__modalites - 1.
 166         if self.__use_mean_std_exp :
 167             self.__ddl = self.__ddl - 2.
 168         if (self.__ddl < 1.):
 169             raise ValueError("The ddl is 0, you must increase the number of classes nbclasse ")
 170         logging.debug("Nombre de degres de liberte=%s"%self.__ddl)
 171
 172     def ComputeValue(self) :
 173         """
 174         Calcul de la variable Q qui suit une loi Khi-2
 175         """
 176         if self.__testHomogen :
 177           kobs,ktheo=self.__kobs.tolist(),self.__kobsHomogen.tolist()
 178         else :
 179           kobs,ktheo=self.__kobs.tolist(),self.__ktheo.tolist()
 180
 181         # on supprime les classes theoriques qui ont moins d'un element (sinon la distance khi2 tendrait vers l'infini)
 182         ko,kt=[],[]
 183         self.__count0 = 0.
 184         for k,val in enumerate(ktheo):
 185             if val > 1.0:
 186                 kt.append(val)
 187                 ko.append(kobs[k])
 188             else :
 189                 self.__count0 = self.__count0 +1.
 190         logging.debug("WARNING : nombre de classes vides supprimees (effectif theorique inferieur a 1.) pour le calcul de la valeur du Khi2 = %s"%self.__count0)
 191         ef1,ef2=asarray(ko),asarray(kt)
 192         count  = 0.
 193         for el in ef1.tolist() :
 194            if el < 5. :
 195               count = count +1.
 196         for el in ef2.tolist() :
 197            if el < 5. :
 198               count = count +1.
 199         pourcent_nbclasse_effecinf = count /(2.*len(ef1.tolist())) *100.
 200         if (pourcent_nbclasse_effecinf > 20.) :
 201            logging.debug("WARNING : nombre de classes dont l effectif est inferieur a 5 elements %s"%pourcent_nbclasse_effecinf)
 202         k,p = chisquare(ef1, ef2)
 203         k2, p2 = [k],[p]
 204         for shift in range(1,6):
 205             k,p=chisquare(ef1[shift:],ef2[:-shift])
 206             k2.append(k)
 207             p2.append(p)
 208             k,p=chisquare(ef1[:-shift],ef2[shift:])
 209             k2.append(k)
 210             p2.append(p)
 211         logging.debug("Liste des valeurs du Khi2 = %s"%k2)
 212         self.__khi2=min(k2)
 213         self.__Q=self.__khi2
 214
 215         logging.debug("Valeur du Khi2=%s"%self.__Q)
 216         return self.__Q
 217
 218     def ComputeArea(self):
 219         """
 220         Calcul de la p-value
 221         """
 222         self.__areakhi2 = 100 * chisqprob(self.__Q, self.__ddl)
 223         return self.__areakhi2
 224
 225     def WriteMessage(self):
 226         """
 227         Interpretation du test
 228         """
 229         message = "Il y a %.2f%s de chance de se tromper en refusant l'adequation"%(self.__areakhi2,"%")
 230         return message
 231
 232     def WriteMessageHomogen(self):
 233         """
 234         Interpretation du test
 235         """
 236         message = "Il y a %.2f%s de chance de se tromper en refusant l'homogeneite"%(self.__areakhi2,"%")
 237         return message
 238
 239 # ==============================================================================
 240 def ComputeKhi2_testGauss(
 241         meantheo = 0.,
 242         stdtheo = 1.,
 243         nech = 10,
 244         dx = 0.1,
 245         nbclasses = None,
 246         SuppressEmptyClasses = True,
 247         ):
 248     """
 249     Test du Khi2 d adequation entre tirage aleatoire dans gaussienne et une gaussienne theo
 250     """
 251     essai = StatspourTests( cdftheo=norm, meantheo = meantheo, stdtheo = stdtheo, pdftest = ["random.normal","(%.3f,%.2f,%d)"%(meantheo,stdtheo,nech)], obs = None, use_mean_std_exp=False,dxmin=dx, obsHomogen = None, nbclasses = nbclasses)
 252     essai.MakeClasses()
 253     essai.ComputeObs()
 254     essai.ComputeTheo()
 255     classes,eftheo, efobs = essai.Computepdfs()
 256     essai.Computeddl()
 257     valeurKhi2= essai.ComputeValue()
 258     areaKhi2 = essai.ComputeArea()
 259     message = essai.WriteMessage()
 260     logging.debug("message %s"%message)
 261     return classes, eftheo, efobs, valeurKhi2, areaKhi2, message
 262
 263 def ComputeKhi2_Gauss(
 264         vectorV = None,
 265         dx = 0.1,
 266         SuppressEmptyClasses = True,
 267         nbclasses = None
 268         ):
 269     """
 270     Test du Khi2 d adequation entre un vecteur donne et une gaussienne theo de mean et std celles du vecteur
 271     """
 272     essai = StatspourTests( cdftheo=norm, pdftest = None, obs = vectorV, use_mean_std_exp=True,dxmin=dx, obsHomogen = None, nbclasses = nbclasses)
 273     essai.MakeClasses()
 274     essai.ComputeObs()
 275     essai.ComputeTheo()
 276     classes,eftheo, efobs = essai.Computepdfs()
 277     essai.Computeddl()
 278     valeurKhi2= essai.ComputeValue()
 279     areaKhi2 = essai.ComputeArea()
 280     message = essai.WriteMessage()
 281     logging.debug("message %s"%message)
 282     return classes, eftheo, efobs, valeurKhi2, areaKhi2, message
 283
 284 def ComputeKhi2_Homogen(
 285         vectorV1 = None,
 286         vectorV2 = None,
 287         dx = 0.1,
 288         SuppressEmptyClasses = True,
 289         nbclasses = None
 290         ):
 291     """
 292     Test du Khi2 d homogeniete entre 2 vecteurs
 293     """
 294     essai = StatspourTests( cdftheo=norm, pdftest = None, obs = vectorV1, use_mean_std_exp=True,dxmin=dx, obsHomogen = vectorV2, nbclasses = nbclasses)
 295     essai.MakeClasses()
 296     essai.ComputeObs()
 297     essai.ComputeObsHomogen()
 298     classes,eftheo, efobs = essai.Computepdfs()
 299     essai.Computeddl()
 300     valeurKhi2= essai.ComputeValue()
 301     areaKhi2 = essai.ComputeArea()
 302     message = essai.WriteMessageHomogen()
 303     logging.debug("message %s"%message)
 304     return classes, eftheo, efobs, valeurKhi2, areaKhi2, message
 305
 306 # ==============================================================================
 307 if __name__ == "__main__":
 308     print '\n AUTODIAGNOSTIC \n'
 309     #
 310     numpy.random.seed(100)
 311
 312     # Test de verification d adequation entre une gaussienne et un tirage gaussien
 313     print ''
 314     print 'Test de verification d adequation entre une gaussienne centree normale et un tirage gaussien'
 315     classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_testGauss(meantheo = 0., stdtheo = 1., nech = 1000., dx = 0.1, SuppressEmptyClasses = True, nbclasses = None)
 316     print '  valeurKhi2=',valeurKhi2
 317     print '  areaKhi2=',areaKhi2
 318     print ' ',message
 319
 320     if (numpy.abs(areaKhi2 - 99.91)< 1.e-2) :
 321        print "The computation of the khisquare value is OK"
 322     else :
 323        raise ValueError("The computation of the khisquare value is WRONG")
 324
 325     numpy.random.seed(2490)
 326
 327     # Test de verification d adequation entre une gaussienne et un vecteur donne
 328     print ''
 329     print 'Test de verification d adequation entre une gaussienne et un vecteur donne'
 330     V = random.normal(50.,1.5,1000)
 331     classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = 0.1, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
 332     print '  valeurKhi2=',valeurKhi2
 333     print '  areaKhi2=',areaKhi2
 334     print ' ',message
 335
 336     if (numpy.abs(areaKhi2 - 99.60)< 1.e-2) :
 337        print "The computation of the khisquare value is OK"
 338     else :
 339        raise ValueError("The computation of the khisquare value is WRONG")
 340
 341     # Test de d homogeneite entre 2 vecteurs donnes
 342     print ''
 343     print 'Test d homogeneite entre 2 vecteurs donnes'
 344     V1 = random.normal(50.,1.5,10000)
 345     numpy.random.seed(2490)
 346     V2 = random.normal(50.,1.5,10000)
 347     classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Homogen(dx = 0.5, vectorV1 = V1, vectorV2 = V2, SuppressEmptyClasses = True, nbclasses = None)
 348     print '  valeurKhi2=',valeurKhi2
 349     print '  areaKhi2=',areaKhi2
 350     print ' ',message
 351
 352     if (numpy.abs(areaKhi2 - 99.98)< 1.e-2) :
 353        print "The computation of the khisquare value is OK"
 354     else :
 355        raise ValueError("The computation of the khisquare value is WRONG")
 356
 357     # Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 10000
 358     print ''
 359     print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 10000'
 360 #    file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech10000.gnu'
 361 #    fid = open(file, "w")
 362 #    lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' )
 363     numpy.random.seed(4000)
 364     V = random.normal(0., 1.,10000)
 365     aire = []
 366     for dx in arange(0.01, 1., 0.001) :
 367       classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
 368 #      lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2)
 369       aire.append(areaKhi2)
 370     meanaire = numpy.asarray(aire)
 371 #    fid.writelines(lines)
 372
 373     print  "  En moyenne, il y a ", meanaire.mean(),"% de chance de se tromper en refusant l adequation a une loi gaussienne  pour un echantillon de taille 10000"
 374     print
 375     if (numpy.abs( meanaire.mean() - 71.79)< 1.e-2) :
 376        print "The computation of the khisquare value is OK"
 377     else :
 378        raise ValueError("The computation of the khisquare value is WRONG")
 379
 380     # Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 1000
 381     print ''
 382     print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 1000'
 383 #    file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech1000.gnu'
 384 #    fid = open(file, "w")
 385 #    lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' )
 386     numpy.random.seed(4000)
 387     V = random.normal(0., 1.,1000)
 388     aire = []
 389     for dx in arange(0.05, 1., 0.001) :
 390       classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
 391 #      lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2)
 392       aire.append(areaKhi2)
 393     meanaire = numpy.asarray(aire)
 394 #    fid.writelines(lines)
 395
 396     print  "  En moyenne, il y a ", meanaire.mean(),"% de chance de se tromper en refusant l adequation a une loi gaussienne  pour un echantillon de taille 1000"
 397     print
 398     if (numpy.abs( meanaire.mean() - 90.60)< 1.e-2) :
 399        print "The computation of the khisquare value is OK"
 400     else :
 401        raise ValueError("The computation of the khisquare value is WRONG")
 402
 403    # Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 100
 404     print ''
 405     print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 100'
 406 #    file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech100.gnu'
 407 #    fid = open(file, "w")
 408 #    lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' )
 409     numpy.random.seed(4000)
 410     V = random.normal(0., 1.,100)
 411     aire = []
 412     for dx in arange(0.1, 1., 0.01) :
 413       classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
 414 #      lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2)
 415       aire.append(areaKhi2)
 416     meanaire = numpy.asarray(aire)
 417 #    fid.writelines(lines)
 418
 419     print  "  En moyenne, il y a ", meanaire.mean(),"% de chance de se tromper en refusant l adequation a une loi gaussienne  pour un echantillon de taille 100"
 420     print