1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2010 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
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
25 __author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Octobre 2008"
27 import sys ; sys.path.insert(0, "../daCore")
30 from scipy.stats import ttest_rel, ttest_ind, betai
33 # ==============================================================================
34 def DependantVectors(vector1 = None, vector2 = None, tolerance = 0.05 ):
36 Outil numérique de calcul de la variable de Student pour 2 vecteurs
40 - les deux vecteurs (comme liste, array ou matrix)
41 d'échantillons dont on veut comparer la variance,
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.
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")
58 # Calcul de la p-value du Test de Student
59 # --------------------------------------------------------------------
60 [t, prob] = ttest_rel(V1, V2)
61 areastudent = 100. * prob
63 logging.debug("DEPENDANTVECTORS t = %.3f, areastudent = %.3f"%(t, areastudent))
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)
70 if (areastudent < (100.*tolerance)) :
71 answerTestStudent = False
73 answerTestStudent = True
75 return areastudent, t, answerTestStudent, message
77 # ==============================================================================
78 def IndependantVectorsDifferentVariance(vector1 = None, vector2 = None, tolerance = 0.05 ):
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.
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")
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)
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
108 logging.debug("IndependantVectorsDifferentVariance t = %.3f, areastudent = %.3f"%(t, areastudent))
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
117 answerTestStudent = True
119 return areastudent, t, answerTestStudent, message
121 # ==============================================================================
122 def IndependantVectorsEqualVariance(vector1 = None, vector2 = None, tolerance = 0.05 ):
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.
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")
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
142 logging.debug("IndependantVectorsEqualVariance t = %.3f, areastudent = %.3f"%(t, areastudent))
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
150 answerTestStudent = True
152 return areastudent, t, answerTestStudent, message
154 # ==============================================================================
155 if __name__ == "__main__":
156 print '\n AUTODIAGNOSTIC \n'
157 # logging.getLogger().setLevel(logging.DEBUG)
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
168 # Appel de l outil DependantVectors et initialisation des inputs
169 [aire, Q, reponse, message] = DependantVectors(
174 # Verification par les calculs sans les routines de scipy.stats
175 # (ref numerical recipes)
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()
186 for j in range(0, n) :
187 cov = cov + (V1[j] - m1)*(V2[j] - m2)
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"
195 raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes")
197 if (numpy.abs(aire - 57.99159)< 1.e-5) :
198 print " Le calcul est JUSTE sur cet exemple."
200 raise ValueError("Le calcul est FAUX sur cet exemple.")
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
211 # Appel de l outil IndependantVectorsDifferentVariance et initialisation des inputs
212 [aire, Q, reponse, message] = IndependantVectorsDifferentVariance(
217 if (numpy.abs(aire - 78.91339)< 1.e-5) :
218 print " Le calcul est JUSTE sur cet exemple."
220 raise ValueError("Le calcul est FAUX sur cet exemple.")
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
231 # Appel de l outil IndependantVectorsEqualVariance et initialisation des inputs
232 [aire, Q, reponse, message] = IndependantVectorsEqualVariance(
237 # Verification par les calculs sans les routines de scipy.stats (ref numerical recipes)
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()
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))
250 if (aireverif - aire < 1.e-5) :
251 print " Le calcul est conforme à celui de l'algorithme du Numerical Recipes"
253 raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes")
255 if (numpy.abs(aire - 78.42572)< 1.e-5) :
256 print " Le calcul est JUSTE sur cet exemple."
258 raise ValueError("Le calcul est FAUX sur cet exemple.")