src/daComposant/daNumerics/ComputeStudent.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2010  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 numérique de calcul des variables de Student pour 2 vecteurs
  23     dépendants ou indépendants, avec variances supposées égales ou différentes
  24 """
  25 __author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Octobre 2008"
  26
  27 import sys ; sys.path.insert(0, "../daCore")
  28
  29 import numpy
  30 from scipy.stats import ttest_rel, ttest_ind, betai
  31 import logging
  32
  33 # ==============================================================================
  34 def DependantVectors(vector1 = None, vector2 = None, tolerance = 0.05 ):
  35     """
  36     Outil numérique de calcul de la variable de Student pour 2 vecteurs
  37     dépendants
  38     Ce calcul nécessite :
  39         - en input :
  40             - les deux vecteurs (comme liste, array ou matrix)
  41               d'échantillons dont on veut comparer la variance,
  42             - la tolérance
  43         - en output :
  44             - la p-value,
  45             - la valeur de la variable aléatoire,
  46             - la reponse au test pour une tolerance ainsi que
  47             - le message qui interprete la reponse du test.
  48     """
  49     if (vector1 is None) or (vector2 is None) :
  50         raise ValueError("Two vectors must be given to calculate the Student value")
  51     V1 = numpy.array(vector1)
  52     V2 = numpy.array(vector2)
  53     if (V1.size < 1) or (V2.size < 1):
  54         raise ValueError("The given vectors must not be empty")
  55     if V1.size != V2.size:
  56         raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
  57     #
  58     # Calcul de la p-value du Test de Student
  59     # --------------------------------------------------------------------
  60     [t, prob] = ttest_rel(V1, V2)
  61     areastudent = 100. * prob
  62     #
  63     logging.debug("DEPENDANTVECTORS t = %.3f, areastudent = %.3f"%(t, areastudent))
  64     #
  65     # Tests
  66     # --------------------------------------------------------------------
  67     message =  "DEPENDANTVECTORS Il y a %.2f %s de chance de se tromper en refusant l'hypothèse d'égalité des moyennes des 2 échantillons dépendants (si <%.2f %s on refuse effectivement l'égalité)"%(areastudent, "%", 100.*tolerance,"%")
  68     logging.debug(message)
  69     #
  70     if (areastudent < (100.*tolerance)) :
  71         answerTestStudent = False
  72     else:
  73         answerTestStudent = True
  74     #
  75     return  areastudent, t, answerTestStudent, message
  76
  77 # ==============================================================================
  78 def IndependantVectorsDifferentVariance(vector1 = None, vector2 = None, tolerance = 0.05 ):
  79     """
  80     Outil numerique de calcul de la variable de Student pour 2 vecteurs independants supposes de variances vraies differentes
  81     En input : la tolerance
  82     En output : la p-value, la valeur de la variable aleatoire, la reponse au test pour une tolerance ainsi que le message qui interprete la reponse du test.
  83     """
  84     if (vector1 is None) or (vector2 is None) :
  85         raise ValueError("Two vectors must be given to calculate the Student value")
  86     V1 = numpy.array(vector1)
  87     V2 = numpy.array(vector2)
  88     if (V1.size < 1) or (V2.size < 1):
  89         raise ValueError("The given vectors must not be empty")
  90     #
  91     # Calcul de la p-value du Test de Student
  92     # --------------------------------------------------------------------
  93     # t = (m1 - m2)/ sqrt[ (var1/n1 + var2/n2) ]
  94     # ou var est calcule comme var = somme (xi -xmena)**2 /(n-1)
  95     n1 = V1.size
  96     n2 = V2.size
  97     mean1 = V1.mean()
  98     mean2 = V2.mean()
  99     var1 = numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std() * numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std()
 100     var2 = numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std() * numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std()
 101     t = (mean1 - mean2)/ numpy.sqrt( var1/n1 + var2/n2 )
 102     df = ( (var1/n1 + var2/n2) * (var1/n1 + var2/n2) ) / ( (var1/n1)*(var1/n1)/(n1-1)  + (var2/n2)*(var2/n2)/(n2-1) )
 103     zerodivproblem = var1/n1 + var2/n2 == 0
 104     t = numpy.where(zerodivproblem, 1.0, t)           # replace NaN t-values with 1.0
 105     prob = betai(0.5*df,0.5,float(df)/(df+t*t))
 106     areastudent = 100. * prob
 107     #
 108     logging.debug("IndependantVectorsDifferentVariance t = %.3f, areastudent = %.3f"%(t, areastudent))
 109     #
 110     # Tests
 111     # --------------------------------------------------------------------
 112     message =  "IndependantVectorsDifferentVariance Il y a %.2f %s de chance de se tromper en refusant l'hypothèse d'égalité des moyennes des 2 échantillons indépendants supposés de variances différentes (si <%.2f %s on refuse effectivement l'égalité)"%(areastudent, "%", 100.* tolerance,"%")
 113     logging.debug(message)
 114     if (areastudent < (100.*tolerance)) :
 115         answerTestStudent = False
 116     else:
 117         answerTestStudent = True
 118     #
 119     return  areastudent, t, answerTestStudent, message
 120
 121 # ==============================================================================
 122 def IndependantVectorsEqualVariance(vector1 = None, vector2 = None, tolerance = 0.05 ):
 123     """
 124     Outil numerique de calcul de la variable de Student pour 2 vecteurs independants supposes de meme variance vraie
 125     En input : la tolerance
 126     En output : la p-value, la valeur de la variable aleatoire, la reponse au test pour une tolerance ainsi que le message qui interprete la reponse du test.
 127     """
 128     if (vector1 is None) or (vector2 is None) :
 129         raise ValueError("Two vectors must be given to calculate the Student value")
 130     V1 = numpy.array(vector1)
 131     V2 = numpy.array(vector2)
 132     if (V1.size < 1) or (V2.size < 1):
 133         raise ValueError("The given vectors must not be empty")
 134     #
 135     # Calcul de la p-value du Test de Student
 136     # --------------------------------------------------------------------
 137     # t = sqrt(n1+n2-2) * (m1 - m2)/ sqrt[ (1/n1 +1/n2) * ( (n1-1)var1 + (n2-1)var2 )]
 138     # ou var est calcule comme var = somme (xi -xmena)**2 /(n-1)
 139     [t, prob] = ttest_ind(V1, V2)
 140     areastudent = 100. * prob
 141     #
 142     logging.debug("IndependantVectorsEqualVariance t = %.3f, areastudent = %.3f"%(t, areastudent))
 143     # Tests
 144     # --------------------------------------------------------------------
 145     message =  "IndependantVectorsEqualVariance Il y a %.2f %s de chance de se tromper en refusant l'hypothèse d'égalité des moyennes des 2 échantillons indépendants supposés de même variance (si <%.2f %s on refuse effectivement l'égalité)"%(areastudent, "%", 100.* tolerance,"%")
 146     logging.debug(message)
 147     if (areastudent < (100.*tolerance)) :
 148         answerTestStudent = False
 149     else:
 150         answerTestStudent = True
 151
 152     return  areastudent, t, answerTestStudent, message
 153
 154 # ==============================================================================
 155 if __name__ == "__main__":
 156     print '\n AUTODIAGNOSTIC \n'
 157     # logging.getLogger().setLevel(logging.DEBUG)
 158
 159     print
 160     print " Test de Student pour des vecteurs dépendants"
 161     print " --------------------------------------------"
 162     # Tirage de l'echantillon
 163     V1 = numpy.matrix(([-1., 0., 4.])).T
 164     V2 = numpy.matrix(([-2., 0., 8.])).T
 165     V1 = V1.A1
 166     V2 = V2.A1
 167     #
 168     # Appel de l outil DependantVectors et initialisation des inputs
 169     [aire, Q, reponse, message] = DependantVectors(
 170         vector1 = V1,
 171         vector2 = V2,
 172         tolerance = 0.05)
 173     #
 174     # Verification par les calculs sans les routines de scipy.stats
 175     # (ref numerical recipes)
 176     n = V1.size
 177     df= n -1
 178     # Les routines de scipy.stats utilisent une variance calculee avec n-1 et non n comme dans std
 179     # t =  (m1 - m2)/ sqrt[(varx1 + varx2 - 2 cov(x1, x2))/n ]
 180     # ou var est calcule comme var = somme (xi -xmean)**2 /(n-1)
 181     var1 = numpy.sqrt(n)/numpy.sqrt(n-1)* V1.std() * numpy.sqrt(n)/numpy.sqrt(n-1) * V1.std()
 182     var2 = numpy.sqrt(n)/numpy.sqrt(n-1)* V2.std() * numpy.sqrt(n)/numpy.sqrt(n-1) * V2.std()
 183     m1 = V1.mean()
 184     m2 = V2.mean()
 185     cov = 0.
 186     for j in range(0, n) :
 187        cov = cov + (V1[j] - m1)*(V2[j] - m2)
 188     cov = cov /df
 189     sd = numpy.sqrt((var1 + var2 - 2. *cov) / n)
 190     tverif = (m1 -m2) /sd
 191     aireverif = 100. * betai(0.5*df,0.5,float(df)/(df+tverif*tverif))
 192     if  (aireverif - aire < 1.e-5)   :
 193        print " Le calcul est conforme à celui de l'algorithme du Numerical Recipes"
 194     else :
 195        raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes")
 196
 197     if  (numpy.abs(aire - 57.99159)< 1.e-5)  :
 198        print " Le calcul est JUSTE sur cet exemple."
 199     else :
 200        raise ValueError("Le calcul est FAUX sur cet exemple.")
 201
 202     print
 203     print " Test de Student pour des vecteurs independants supposés de même variance"
 204     print " ------------------------------------------------------------------------"
 205     # Tirage de l'echantillon
 206     V1 = numpy.matrix(([-1., 0., 4.])).T
 207     V2 = numpy.matrix(([-2., 0., 8.])).T
 208     V1 = V1.A1
 209     V2 = V2.A1
 210     #
 211     # Appel de l outil IndependantVectorsDifferentVariance et initialisation des inputs
 212     [aire, Q, reponse, message] = IndependantVectorsDifferentVariance(
 213         vector1 = V1,
 214         vector2 = V2,
 215         tolerance = 0.05)
 216     #
 217     if  (numpy.abs(aire - 78.91339)< 1.e-5)  :
 218        print " Le calcul est JUSTE sur cet exemple."
 219     else :
 220        raise ValueError("Le calcul est FAUX sur cet exemple.")
 221
 222     print
 223     print " Test de Student pour des vecteurs indépendants supposés de même variance"
 224     print " ------------------------------------------------------------------------"
 225     # Tirage de l'echantillon
 226     V1 = numpy.matrix(([-1., 0., 4.])).T
 227     V2 = numpy.matrix(([-2., 0., 8.])).T
 228     V1 = V1.A1
 229     V2 = V2.A1
 230     #
 231     # Appel de l outil IndependantVectorsEqualVariance et initialisation des inputs
 232     [aire, Q, reponse, message] = IndependantVectorsEqualVariance(
 233         vector1 = V1,
 234         vector2 = V2,
 235         tolerance = 0.05)
 236     #
 237     # Verification par les calculs sans les routines de scipy.stats (ref numerical recipes)
 238     n1 = V1.size
 239     n2 = V2.size
 240     df= n1 + n2 -2
 241     # Les routines de scipy.stats utilisent une variance calculee avec n-1 et non n comme dans std
 242     var1 = numpy.sqrt(n1)/numpy.sqrt(n1-1)* V1.std() * numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std()
 243     var2 = numpy.sqrt(n2)/numpy.sqrt(n2-1)* V2.std() * numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std()
 244     m1 = V1.mean()
 245     m2 = V2.mean()
 246     var = ((n1 -1.) *var1 + (n2 -1.) *var2 ) /df
 247     tverif = (m1 -m2) /numpy.sqrt(var*(1./n1 + 1./n2))
 248     aireverif = 100. * betai(0.5*df,0.5,float(df)/(df+tverif*tverif))
 249     #
 250     if  (aireverif - aire < 1.e-5)   :
 251        print " Le calcul est conforme à celui de l'algorithme du Numerical Recipes"
 252     else :
 253        raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes")
 254
 255     if  (numpy.abs(aire - 78.42572)< 1.e-5)  :
 256        print " Le calcul est JUSTE sur cet exemple."
 257     else :
 258        raise ValueError("Le calcul est FAUX sur cet exemple.")
 259
 260     print