self._name = "3DVAR"
logging.debug("%s Initialisation"%self._name)
- def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None):
+ def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
"""
Calcul de l'estimateur 3D-VAR
"""
else:
logging.debug("%s Calcul de Hm(Xb)"%self._name)
HXb = Hm( Xb )
+ HXb = numpy.asmatrix(HXb).flatten().T
#
# Calcul du préconditionnement
# ----------------------------
logging.info("%s CostFunction Jb = %s"%(self._name, Jb))
logging.info("%s CostFunction Jo = %s"%(self._name, Jo))
logging.info("%s CostFunction J = %s"%(self._name, J))
+ self.StoredVariables["CurrentState"].store( _X.A1 )
self.StoredVariables["CostFunctionJb"].store( Jb )
self.StoredVariables["CostFunctionJo"].store( Jo )
self.StoredVariables["CostFunctionJ" ].store( J )
logging.debug("%s GradientOfCostFunction GradJb = %s"%(self._name, numpy.asmatrix( GradJb ).flatten()))
logging.debug("%s GradientOfCostFunction GradJo = %s"%(self._name, numpy.asmatrix( GradJo ).flatten()))
logging.debug("%s GradientOfCostFunction GradJ = %s"%(self._name, numpy.asmatrix( GradJ ).flatten()))
- # self.StoredVariables["GradientOfCostFunctionJb"].store( Jb )
- # self.StoredVariables["GradientOfCostFunctionJo"].store( Jo )
- # self.StoredVariables["GradientOfCostFunctionJ" ].store( J )
return GradJ.A1
#
# Point de démarrage de l'optimisation : Xini = Xb
#
# Paramètres de pilotage
# ----------------------
- if Par.has_key("Bounds") and (type(Par["Bounds"]) is type([]) or type(Par["Bounds"]) is type(())) and (len(Par["Bounds"]) > 0):
- Bounds = Par["Bounds"]
+ if Parameters.has_key("Bounds") and (type(Parameters["Bounds"]) is type([]) or type(Parameters["Bounds"]) is type(())) and (len(Parameters["Bounds"]) > 0):
+ Bounds = Parameters["Bounds"]
else:
Bounds = None
MinimizerList = ["LBFGSB","TNC", "CG", "BFGS"]
- if Par.has_key("Minimizer") and (Par["Minimizer"] in MinimizerList):
- Minimizer = str( Par["Minimizer"] )
+ if Parameters.has_key("Minimizer") and (Parameters["Minimizer"] in MinimizerList):
+ Minimizer = str( Parameters["Minimizer"] )
else:
Minimizer = "LBFGSB"
logging.debug("%s Minimiseur utilisé = %s"%(self._name, Minimizer))
- if Par.has_key("MaximumNumberOfSteps") and (Par["MaximumNumberOfSteps"] > -1):
- maxiter = int( Par["MaximumNumberOfSteps"] )
+ if Parameters.has_key("MaximumNumberOfSteps") and (Parameters["MaximumNumberOfSteps"] > -1):
+ maxiter = int( Parameters["MaximumNumberOfSteps"] )
else:
maxiter = 15000
logging.debug("%s Nombre maximal de pas d'optimisation = %s"%(self._name, maxiter))
bounds = Bounds,
maxfun = maxiter,
iprint = iprint,
+ factr = 1.,
)
logging.debug("%s %s Minimum = %s"%(self._name, Minimizer, Minimum))
logging.debug("%s %s Nb of F = %s"%(self._name, Minimizer, Informations['funcalls']))
from daCore import BasicObjects, PlatformInfo
m = PlatformInfo.SystemUsage()
+import numpy
+
# ==============================================================================
class ElementaryAlgorithm(BasicObjects.Algorithm):
def __init__(self):
self._name = "BLUE"
logging.debug("%s Initialisation"%self._name)
- def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None):
+ def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
"""
Calcul de l'estimateur BLUE (ou Kalman simple, ou Interpolation Optimale)
"""
else:
logging.debug("%s Calcul de Hm * Xb"%self._name)
HXb = Hm * Xb
-
+ HXb = numpy.asmatrix(HXb).flatten().T
+ #
# Calcul de la matrice de gain dans l'espace le plus petit
+ # --------------------------------------------------------
if Y.size <= Xb.size:
logging.debug("%s Calcul de K dans l'espace des observations"%self._name)
K = B * Ht * (Hm * B * Ht + R).I
self._name = "ENSEMBLEBLUE"
logging.debug("%s Initialisation"%self._name)
- def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None ):
+ def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None ):
"""
Calcul d'une estimation BLUE d'ensemble :
- génération d'un ensemble d'observations, de même taille que le
self._name = "KALMAN"
logging.debug("%s Initialisation"%self._name)
- def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None):
+ def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
"""
Calcul de l'estimateur de Kalman
"""
BasicObjects.Algorithm.__init__(self)
self._name = "LINEARLEASTSQUARES"
- def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None):
+ def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
"""
Calcul de l'estimateur au sens des moindres carres sans ebauche
"""
self.__M = {}
#
self.__X = Persistence.OneVector()
- self.__Parameters = {}
+ self.__Parameters = {}
self.__StoredDiagnostics = {}
+ self.__StoredInputs = {}
#
# Variables temporaires
- self.__algorithm = {}
- self.__algorithmFile = None
- self.__algorithmName = None
- self.__diagnosticFile = None
+ self.__algorithm = {}
+ self.__algorithmFile = None
+ self.__algorithmName = None
+ self.__diagnosticFile = None
#
# Récupère le chemin du répertoire parent et l'ajoute au path
# (Cela complète l'action de la classe PathManagement dans PlatformInfo,
asVector = None,
asPersistentVector = None,
Scheduler = None,
+ toBeStored = False,
):
"""
Permet de définir l'estimation a priori :
- asPersistentVector : entrée des données, comme un vecteur de type
persistent contruit avec la classe ad-hoc "Persistence"
- Scheduler est le contrôle temporel des données
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
if asVector is not None:
if type( asVector ) is type( numpy.matrix([]) ):
self.__Xb = asPersistentVector
else:
raise ValueError("Error: improperly defined background")
+ if toBeStored:
+ self.__StoredInputs["Background"] = self.__Xb
return 0
- def setBackgroundError(self, asCovariance=None):
+ def setBackgroundError(self,
+ asCovariance = None,
+ toBeStored = False,
+ ):
"""
Permet de définir la covariance des erreurs d'ébauche :
- asCovariance : entrée des données, comme une matrice compatible avec
le constructeur de numpy.matrix
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
self.__B = numpy.matrix( asCovariance, numpy.float )
+ if toBeStored:
+ self.__StoredInputs["BackgroundError"] = self.__B
return 0
# -----------------------------------------------------------
asVector = None,
asPersistentVector = None,
Scheduler = None,
+ toBeStored = False,
):
"""
Permet de définir les observations :
- asPersistentVector : entrée des données, comme un vecteur de type
persistent contruit avec la classe ad-hoc "Persistence"
- Scheduler est le contrôle temporel des données disponibles
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
if asVector is not None:
if type( asVector ) is type( numpy.matrix([]) ):
self.__Y = asPersistentVector
else:
raise ValueError("Error: improperly defined observations")
+ if toBeStored:
+ self.__StoredInputs["Observation"] = self.__Y
return 0
- def setObservationError(self, asCovariance=None):
+ def setObservationError(self,
+ asCovariance = None,
+ toBeStored = False,
+ ):
"""
Permet de définir la covariance des erreurs d'observations :
- asCovariance : entrée des données, comme une matrice compatible avec
le constructeur de numpy.matrix
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
self.__R = numpy.matrix( asCovariance, numpy.float )
+ if toBeStored:
+ self.__StoredInputs["ObservationError"] = self.__R
return 0
def setObservationOperator(self,
asFunction = {"Direct":None, "Tangent":None, "Adjoint":None},
asMatrix = None,
appliedToX = None,
+ toBeStored = False,
):
"""
Permet de définir un opérateur d'observation H. L'ordre de priorité des
X divers, l'opérateur par sa valeur appliquée à cet X particulier,
sous la forme d'un dictionnaire appliedToX[NAME] avec NAME un nom.
L'opérateur doit néanmoins déjà avoir été défini comme d'habitude.
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
if (type(asFunction) is type({})) and (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
else:
self.__H["AppliedToX"] = None
#
+ if toBeStored:
+ self.__StoredInputs["ObservationOperator"] = self.__H
return 0
# -----------------------------------------------------------
asFunction = {"Direct":None, "Tangent":None, "Adjoint":None},
asMatrix = None,
Scheduler = None,
+ toBeStored = False,
):
"""
Permet de définir un opérateur d'évolution M. L'ordre de priorité des
la matrice, et l'opérateur "Adjoint" à l'aide de la transposée. La
matrice fournie doit être sous une forme compatible avec le
constructeur de numpy.matrix.
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
if (type(asFunction) is type({})) and (asFunction["Tangent"] is not None) and (asFunction["Adjoint"] is not None):
if not asFunction.has_key("Direct") or (asFunction["Direct"] is None):
self.__M["Adjoint"] = Operator( fromMatrix = matrice.T )
else:
raise ValueError("Error: improperly defined evolution operator")
+ #
+ if toBeStored:
+ self.__StoredInputs["EvolutionModel"] = self.__M
return 0
- def setEvolutionError(self, asCovariance=None):
+ def setEvolutionError(self,
+ asCovariance = None,
+ toBeStored = False,
+ ):
"""
Permet de définir la covariance des erreurs de modèle :
- asCovariance : entrée des données, comme une matrice compatible avec
le constructeur de numpy.matrix
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
self.__Q = numpy.matrix( asCovariance, numpy.float )
+ if toBeStored:
+ self.__StoredInputs["EvolutionError"] = self.__Q
return 0
# -----------------------------------------------------------
- def setControls (self, asVector = None ):
+ def setControls (self,
+ asVector = None,
+ toBeStored = False,
+ ):
"""
Permet de définir la valeur initiale du vecteur X contenant toutes les
variables de contrôle, i.e. les paramètres ou l'état dont on veut
algorithme itératif/incrémental
- asVector : entrée des données, comme un vecteur compatible avec le
constructeur de numpy.matrix.
+ - toBeStored : booléen indiquant si la donnée d'entrée est sauvée pour
+ être rendue disponible au même titre que les variables de calcul
"""
if asVector is not None:
self.__X.store( asVector )
+ if toBeStored:
+ self.__StoredInputs["Controls"] = self.__X
return 0
# -----------------------------------------------------------
# Instancie un objet du type élémentaire du fichier
# -------------------------------------------------
self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
+ self.__StoredInputs["AlgorithmName"] = str(choice)
return 0
def setAlgorithmParameters(self, asDico=None):
Permet de définir les paramètres de l'algorithme, sous la forme d'un
dictionnaire.
"""
- self.__Parameters = dict( asDico )
+ if asDico is not None:
+ self.__Parameters = dict( asDico )
+ else:
+ self.__Parameters = {}
+ self.__StoredInputs["AlgorithmParameters"] = self.__Parameters
return 0
# -----------------------------------------------------------
#
# Instancie un objet du type élémentaire du fichier
# -------------------------------------------------
- if self.__StoredDiagnostics.has_key(name):
- raise ValueError("A diagnostic with the same name already exists")
+ if self.__StoredInputs.has_key(name):
+ raise ValueError("A default input with the same name \"%s\" already exists."%str(name))
+ elif self.__StoredDiagnostics.has_key(name):
+ raise ValueError("A diagnostic with the same name \"%s\" already exists."%str(name))
else:
self.__StoredDiagnostics[name] = self.__diagnosticFile.ElementaryDiagnostic(
name = name,
self.shape_validate()
#
self.__algorithm.run(
- Xb = self.__Xb,
- Y = self.__Y,
- H = self.__H,
- M = self.__M,
- R = self.__R,
- B = self.__B,
- Q = self.__Q,
- Par = self.__Parameters,
+ Xb = self.__Xb,
+ Y = self.__Y,
+ H = self.__H,
+ M = self.__M,
+ R = self.__R,
+ B = self.__B,
+ Q = self.__Q,
+ Parameters = self.__Parameters,
)
return 0
if key is not None:
if self.__algorithm.has_key(key):
return self.__algorithm.get( key )
+ elif self.__StoredInputs.has_key(key):
+ return self.__StoredInputs[key]
elif self.__StoredDiagnostics.has_key(key):
return self.__StoredDiagnostics[key]
else:
- raise ValueError("The requested key \"%s\" does not exists as a diagnostic or as a stored variable."%key)
+ raise ValueError("The requested key \"%s\" does not exists as an input, a diagnostic or a stored variable."%key)
else:
allvariables = self.__algorithm.get()
allvariables.update( self.__StoredDiagnostics )
+ allvariables.update( self.__StoredInputs )
return allvariables
- def get_available_variables(self, key=None):
+ def get_available_variables(self):
"""
Renvoie les variables potentiellement utilisables pour l'étude,
+ initialement stockées comme données d'entrées ou dans les algorithmes,
identifiés par les chaînes de caractères. L'algorithme doit avoir été
préalablement choisi sinon la méthode renvoie "None".
"""
-
- if type( self.__algorithm ) is type( {} ):
+ if len( self.__algorithm.keys()) == 0 and len( self.__StoredInputs.keys() ) == 0:
return None
- if key is not None:
- if self.__algorithm.has_key(key):
- return self.__algorithm.get( key )
- else:
- raise ValueError("The requested key \"%s\" does not exists as a stored variable."%key)
else:
- variables = self.__algorithm.get().keys()
+ variables = []
+ if len( self.__algorithm.keys()) > 0:
+ variables.extend( self.__algorithm.get().keys() )
+ if len( self.__StoredInputs.keys() ) > 0:
+ variables.extend( self.__StoredInputs.keys() )
variables.sort()
return variables
HookFunction = HookFunction,
HookParameters = HookParameters,
)
+
def prepare_to_pickle(self):
- self.__algorithmFile = None
- self.__diagnosticFile = None
- self.__H = {}
+ self.__algorithmFile = None
+ self.__diagnosticFile = None
+ self.__H = {}
# ==============================================================================
if __name__ == "__main__":
print "Innovation :", ADD.get("Innovation").valueserie(0)
print
- print "Algorithmes disponibles :", ADD.get_available_algorithms()
- # print " Chemin des algorithmes :", ADD.get_algorithms_main_path()
- print "Diagnostics disponibles :", ADD.get_available_diagnostics()
- # print " Chemin des diagnostics :", ADD.get_diagnostics_main_path()
+ print "Algorithmes disponibles.......................:", ADD.get_available_algorithms()
+ # print " Chemin des algorithmes.....................:", ADD.get_algorithms_main_path()
+ print "Diagnostics types disponibles.................:", ADD.get_available_diagnostics()
+ # print " Chemin des diagnostics.....................:", ADD.get_diagnostics_main_path()
+ print "Variables disponibles.........................:", ADD.get_available_variables()
print
ADD.setDiagnostic("RMS", "Ma RMS")
liste = ADD.get().keys()
liste.sort()
- print "Variables et diagnostics disponibles :", liste
- print
+ print "Variables et diagnostics nommés disponibles...:", liste
+ print
+ print "Exemple de mise en place d'un observeur :"
def obs(var=None,info=None):
print " ---> Mise en oeuvre de l'observer"
print " var =",var.valueserie(-1)
interne à l'objet, mais auquel on accède par la méthode "get".
Les variables prévues sont :
+ - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
+ - CostFunctionJb : partie ébauche ou background de la fonction-cout
+ - CostFunctionJo : partie observations de la fonction-cout
+ - GradientOfCostFunctionJ : gradient de la fonction-cout globale
+ - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
+ - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
+ - CurrentState : état courant lors d'itérations
- Analysis : l'analyse
- Innovation : l'innovation : d = Y - H Xb
- SigmaObs2 : correction optimale des erreurs d'observation
On peut rajouter des variables à stocker dans l'initialisation de
l'algorithme élémentaire qui va hériter de cette classe
"""
+ self._name = None
self.StoredVariables = {}
+ #
self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
- self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneScalar(name = "GradientOfCostFunctionJ")
- self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneScalar(name = "GradientOfCostFunctionJb")
- self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneScalar(name = "GradientOfCostFunctionJo")
+ self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
+ self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
+ self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
+ self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
self.StoredVariables["CovarianceAPosteriori"] = Persistence.OneMatrix(name = "CovarianceAPosteriori")
- self._name = None
def get(self, key=None):
"""
"""
return self.StoredVariables.has_key(key)
- def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Par=None):
+ def keys(self):
+ """
+ Renvoie la liste des clés de variables stockées.
+ """
+ return self.StoredVariables.keys()
+
+ def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
"""
Doit implémenter l'opération élémentaire de calcul d'assimilation sous
sa forme mathématique la plus naturelle possible.
"""
return len(self.__steps)
+ # ---------------------------------------------------------
+ # Méthodes d'accès de type dictionnaire
+ def keys(self):
+ return self.stepserie()
+
+ def values(self):
+ return self.valueserie()
+
+ def items(self):
+ pairs = []
+ for i in xrange(self.stepnumber()):
+ pairs.append( (self.stepserie(item=i), self.valueserie(item=i)) )
+ return pairs
+
# ---------------------------------------------------------
def mean(self):
"""
def __init__(self, name="", unit="", basetype = list):
Persistence.__init__(self, name, unit, basetype)
+def NoType( value ): return value
+
+class OneNoType(Persistence):
+ """
+ Classe définissant le stockage d'un objet sans modification (cast) de type.
+ Attention, selon le véritable type de l'objet stocké à chaque pas, les
+ opérations arithmétiques à base de numpy peuvent être invalides ou donner
+ des résultats inatendus. Cette classe n'est donc à utiliser qu'à bon escient
+ volontairement, et pas du tout par défaut.
+ """
+ def __init__(self, name="", unit="", basetype = NoType):
+ Persistence.__init__(self, name, unit, basetype)
+
+# ==============================================================================
+class CompositePersistence:
+ """
+ Structure de stockage permettant de rassembler plusieurs objets de
+ persistence.
+
+ Des objets par défaut sont prévus, et des objets supplémentaires peuvent
+ être ajoutés.
+ """
+ def __init__(self, name=""):
+ """
+ name : nom courant
+
+ La gestion interne des données est exclusivement basée sur les variables
+ initialisées ici (qui ne sont pas accessibles depuis l'extérieur des
+ objets comme des attributs) :
+ __StoredObjects : objets de type persistence collectés dans cet objet
+ """
+ self.__name = str(name)
+ #
+ self.__StoredObjects = {}
+ #
+ # Definition des objets par defaut
+ # --------------------------------
+ self.__StoredObjects["Background"] = OneVector("Background", basetype=numpy.array)
+ self.__StoredObjects["BackgroundError"] = OneMatrix("BackgroundError")
+ self.__StoredObjects["Observation"] = OneVector("Observation", basetype=numpy.array)
+ self.__StoredObjects["ObservationError"] = OneMatrix("ObservationError")
+ self.__StoredObjects["Analysis"] = OneVector("Analysis", basetype=numpy.array)
+ self.__StoredObjects["AnalysisError"] = OneMatrix("AnalysisError")
+ self.__StoredObjects["Innovation"] = OneVector("Innovation", basetype=numpy.array)
+ self.__StoredObjects["KalmanGainK"] = OneMatrix("KalmanGainK")
+ self.__StoredObjects["OperatorH"] = OneMatrix("OperatorH")
+ self.__StoredObjects["RmsOMA"] = OneScalar("RmsOMA")
+ self.__StoredObjects["RmsOMB"] = OneScalar("RmsOMB")
+ self.__StoredObjects["RmsBMA"] = OneScalar("RmsBMA")
+ #
+
+ def store(self, name=None, value=None, step=None):
+ """
+ Stockage d'une valeur "value" pour le "step" dans la variable "name".
+ """
+ if name is None: raise ValueError("Storable object name is required for storage.")
+ if name not in self.__StoredObjects.keys():
+ raise ValueError("No such name '%s' exists in storable objects."%name)
+ self.__StoredObjects[name].store( value=value, step=step )
+
+ def add_object(self, name=None, persistenceType=Persistence, basetype=numpy.array ):
+ """
+ Ajoute dans les objets stockables un nouvel objet défini par son nom, son
+ type de Persistence et son type de base à chaque pas.
+ """
+ if name is None: raise ValueError("Object name is required for adding an object.")
+ if name in self.__StoredObjects.keys():
+ raise ValueError("An object with the same name '%s' already exists in storable objects. Choose another one."%name)
+ self.__StoredObjects[name] = persistenceType( name=str(name), basetype=basetype )
+
+ def get_object(self, name=None ):
+ """
+ Renvoie l'objet de type Persistence qui porte le nom demandé.
+ """
+ if name is None: raise ValueError("Object name is required for retrieving an object.")
+ if name not in self.__StoredObjects.keys():
+ raise ValueError("No such name '%s' exists in stored objects."%name)
+ return self.__StoredObjects[name]
+
+ def set_object(self, name=None, objet=None ):
+ """
+ Affecte directement un 'objet' qui porte le nom 'name' demandé.
+ Attention, il n'est pas effectué de vérification sur le type, qui doit
+ comporter les méthodes habituelles de Persistence pour que cela
+ fonctionne.
+ """
+ if name is None: raise ValueError("Object name is required for setting an object.")
+ if name in self.__StoredObjects.keys():
+ raise ValueError("An object with the same name '%s' already exists in storable objects. Choose another one."%name)
+ self.__StoredObjects[name] = objet
+
+ def del_object(self, name=None ):
+ """
+ Supprime un objet de la liste des objets stockables.
+ """
+ if name is None: raise ValueError("Object name is required for retrieving an object.")
+ if name not in self.__StoredObjects.keys():
+ raise ValueError("No such name '%s' exists in stored objects."%name)
+ del self.__StoredObjects[name]
+
+ # ---------------------------------------------------------
+ # Méthodes d'accès de type dictionnaire
+ def __getitem__(self, name=None ):
+ return self.get_object( name )
+
+ def __setitem__(self, name=None, objet=None ):
+ self.set_object( name, objet )
+
+ def keys(self):
+ return self.get_stored_objects(hideVoidObjects = False)
+
+ def values(self):
+ return self.__StoredObjects.values()
+
+ def items(self):
+ return self.__StoredObjects.items()
+
+ # ---------------------------------------------------------
+ def get_stored_objects(self, hideVoidObjects = False):
+ objs = self.__StoredObjects.keys()
+ if hideVoidObjects:
+ usedObjs = []
+ for k in objs:
+ try:
+ if len(self.__StoredObjects[k]) > 0: usedObjs.append( k )
+ except:
+ pass
+ objs = usedObjs
+ objs.sort()
+ return objs
+
+ # ---------------------------------------------------------
+ def save_composite(self, filename=None, mode="pickle"):
+ """
+ Enregistre l'objet dans le fichier indiqué selon le "mode" demandé,
+ et renvoi le nom du fichier
+ """
+ import os
+ if filename is None:
+ filename = os.tempnam( os.getcwd(), 'dacp' ) + ".pkl"
+ else:
+ filename = os.path.abspath( filename )
+ #
+ import cPickle
+ if mode is "pickle":
+ output = open( filename, 'wb')
+ cPickle.dump(self, output)
+ output.close()
+ else:
+ raise ValueError("Save mode '%s' unknown. Choose another one."%mode)
+ #
+ return filename
+
+ def load_composite(self, filename=None, mode="pickle"):
+ """
+ Recharge un objet composite sauvé en fichier
+ """
+ import os
+ if filename is None:
+ raise ValueError("A file name if requested to load a composite.")
+ else:
+ filename = os.path.abspath( filename )
+ #
+ import cPickle
+ if mode is "pickle":
+ pkl_file = open(filename, 'rb')
+ output = cPickle.load(pkl_file)
+ for k in output.keys():
+ self[k] = output[k]
+ else:
+ raise ValueError("Load mode '%s' unknown. Choose another one."%mode)
+ #
+ return filename
+
# ==============================================================================
if __name__ == "__main__":
print '\n AUTODIAGNOSTIC \n'
del OBJET_DE_TEST
print
- print "======> Affichage d'objets stockés"
+ print "======> Utilisation des méthodes d'accès de type dictionnaire"
+ OBJET_DE_TEST = OneScalar("My int", unit="cm", basetype=int)
+ for i in range(5):
+ OBJET_DE_TEST.store( 7+i )
+ print "Taille \"len\" :", len(OBJET_DE_TEST)
+ print "Les pas de stockage :", OBJET_DE_TEST.keys()
+ print "Les valeurs :", OBJET_DE_TEST.values()
+ print "Les paires :", OBJET_DE_TEST.items()
+ del OBJET_DE_TEST
+ print
+
+ print "======> Persistence composite"
+ OBJET_DE_TEST = CompositePersistence("My CompositePersistence")
+ print "Objets stockables :", OBJET_DE_TEST.get_stored_objects()
+ print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True )
+ print "--> Stockage d'une valeur de Background"
+ OBJET_DE_TEST.store("Background",numpy.zeros(5))
+ print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True )
+ print "--> Ajout d'un objet nouveau par defaut, de type vecteur numpy par pas"
+ OBJET_DE_TEST.add_object("ValeursVectorielles")
+ OBJET_DE_TEST.store("ValeursVectorielles",numpy.zeros(5))
+ print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True )
+ print "--> Ajout d'un objet nouveau de type liste par pas"
+ OBJET_DE_TEST.add_object("ValeursList", persistenceType=OneList )
+ OBJET_DE_TEST.store("ValeursList",range(5))
+ print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True )
+ print "--> Ajout d'un objet nouveau, de type vecteur string par pas"
+ OBJET_DE_TEST.add_object("ValeursStr", persistenceType=Persistence, basetype=str )
+ OBJET_DE_TEST.store("ValeursStr","c020")
+ OBJET_DE_TEST.store("ValeursStr","c021")
+ print "Les valeurs :", OBJET_DE_TEST.get_object("ValeursStr").valueserie()
+ print "Acces comme dict :", OBJET_DE_TEST["ValeursStr"].stepserie()
+ print "Acces comme dict :", OBJET_DE_TEST["ValeursStr"].valueserie()
+ print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True )
+ print "--> Suppression d'un objet"
+ OBJET_DE_TEST.del_object("ValeursVectorielles")
+ print "Objets actifs :", OBJET_DE_TEST.get_stored_objects( hideVoidObjects = True )
+ print "--> Enregistrement de l'objet complet de Persistence composite"
+ OBJET_DE_TEST.save_composite("composite.pkl")
+ print
+
+ print "======> Affichage graphique d'objets stockés"
OBJET_DE_TEST = Persistence("My object", unit="", basetype=numpy.array)
D = OBJET_DE_TEST
vect1 = [1, 2, 1, 2, 1]
D.store(vect1)
D.store(vect2)
D.store(vect3)
- print "Affichage de l'ensemble du stockage sur une même image"
+ print "Affichage graphique de l'ensemble du stockage sur une même image"
D.stepplot(
title = "Tous les vecteurs",
filename="vecteurs.ps",
pause = False )
print "Stockage d'un quatrième vecteur de longueur différente"
D.store(vect4)
- print "Affichage séparé du dernier stockage"
+ print "Affichage graphique séparé du dernier stockage"
D.plot(
item = 3,
title = "Vecteurs",
del OBJET_DE_TEST
print
- print "======> Affichage dynamique d'objets"
+ print "======> Affichage graphique dynamique d'objets"
OBJET_DE_TEST = Persistence("My object", unit="", basetype=float)
D = OBJET_DE_TEST
D.plot(
# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
#
name = "Data Assimilation Package"
-version = "0.2.0"
-date = "lundi 23 septembre 2009, 11:11:11 (UTC+0200)"
+version = "0.3.0"
+date = "lundi 15 octobre 2010, 11:11:11 (UTC+0200)"
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Calcul du biais (i.e. la moyenne) à chaque pas. Ce diagnostic très simple
- est présent pour rappeller à l'utilisateur de l'assimilation qu'il faut
- qu'il vérifie le biais de ses erreurs en particulier.
-"""
-__author__ = "Sophie RICCI - Aout 2008"
-
-import numpy
-from daCore import BasicObjects, Persistence
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- BasicObjects.Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = float )
-
- def _formula(self, V):
- """
- Calcul du biais, qui est simplement la moyenne du vecteur
- """
- biais = V.mean()
- #
- return biais
-
- def calculate(self, vector = None, step = None):
- """
- Teste les arguments, active la formule de calcul et stocke le résultat
- """
- if vector is None:
- raise ValueError("One vector must be given to compute biais")
- V = numpy.array(vector)
- if V.size < 1:
- raise ValueError("The given vector must not be empty")
- #
- value = self._formula( V)
- #
- self.store( value = value, step = step )
-
-#===============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
- #
- # Instanciation de l'objet diagnostic
- # -----------------------------------
- D = ElementaryDiagnostic("Mon ComputeBiais")
- #
- # Tirage d un vecteur choisi
- # --------------------------
- x = numpy.matrix(([3., 4., 5.]))
- print " Le vecteur de type 'matrix' choisi est..:", x
- print " Le biais attendu de ce vecteur est......:", x.mean()
- #
- D.calculate( vector = x)
- print " Le biais obtenu de ce vecteur est.......:", D.valueserie(0)
- print
- #
- # Tirage d un vecteur choisi
- # --------------------------
- x = numpy.array(range(11))
- print " Le vecteur de type 'array' choisi est...:", x
- print " Le biais attendu de ce vecteur est......:", x.mean()
- #
- D.calculate( vector = x)
- print " Le biais obtenu de ce vecteur est.......:", D.valueserie(1)
- print
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Calcul de la fonction coût
-"""
-__author__ = "Sophie RICCI - Octobre 2008"
-
-import numpy
-from daCore import BasicObjects, Persistence
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- BasicObjects.Diagnostic.__init__(self, name)
- Persistence.OneScalar.__init__( self, name, unit, basetype = float)
-
- def _formula(self, X, HX, Xb, Y, R, B):
- """
- Calcul de la fonction cout
- """
-# Jb = 1./2. * (X - Xb).T * B.I * (X - Xb)
- Jb = 1./2. * numpy.dot((X - Xb) ,numpy.asarray(numpy.dot(B.I,(X - Xb)).A1))
- logging.info( "Partial cost function : Jb = %s"%Jb )
- #
-# Jo = 1./2. * (Y - HX).T * R.I * (Y - HX)
- Jo = 1./2. * numpy.dot((Y - HX) ,numpy.asarray(numpy.dot(R.I,(Y - HX)).A1))
- logging.info( "Partial cost function : Jo = %s"%Jo )
- #
- J = Jb + Jo
- logging.info( "Total cost function : J = Jo + Jb = %s"%J )
- return J
-
- def calculate(self, x = None, Hx = None, xb = None, yo = None, R = None, B = None , step = None):
- """
- Teste les arguments, active la formule de calcul et stocke le résultat
- """
- if (x is None) or (xb is None) or (yo is None) :
- raise ValueError("Vectors x, xb and yo must be given to compute J")
-# if (type(x) is not float) and (type(x) is not numpy.float64) :
-# if (x.size < 1) or (xb.size < 1) or (yo.size < 1):
-# raise ValueError("Vectors x, xb and yo must not be empty")
- if hasattr(numpy.matrix(x),'A1') :
- X = numpy.matrix(x).A1
- if hasattr(numpy.matrix(xb),'A1') :
- Xb = numpy.matrix(xb).A1
- if hasattr(numpy.matrix(yo),'A1') :
- Y = numpy.matrix(yo).A1
- B = numpy.matrix(B)
- R = numpy.matrix(R)
- if (Hx is None ) :
- raise ValueError("The given vector must be given")
-# if (Hx.size < 1) :
-# raise ValueError("The given vector must not be empty")
- HX = Hx.A1
- if (B is None ) or (R is None ):
- raise ValueError("The matrices B and R must be given")
-# if (B.size < 1) or (R.size < 1) :
-# raise ValueError("The matrices B and R must not be empty")
- #
- value = self._formula(X, HX, Xb, Y, R, B)
- #
- self.store( value = value, step = step )
-
-#===============================================================================
-if __name__ == "__main__":
- print "\nAUTOTEST\n"
- #
- D = ElementaryDiagnostic("Ma fonction cout")
- #
- # Vecteur de type array
- # ---------------------
- x = numpy.array([1., 2.])
- xb = numpy.array([2., 2.])
- yo = numpy.array([5., 6.])
- H = numpy.matrix(numpy.identity(2))
-# Hx = H*x
- Hx = numpy.dot(H,x)
- Hx = Hx.T
- B = numpy.matrix(numpy.identity(2))
- R = numpy.matrix(numpy.identity(2))
- #
- D.calculate( x = x, Hx = Hx, xb = xb, yo = yo, R = R, B = B)
- print "Le vecteur x choisi est...:", x
- print "L ebauche xb choisie est...:", xb
- print "Le vecteur d observation est...:", yo
- print "B = ", B
- print "R = ", R
- print "La fonction cout J vaut ...: %.2e"%D.valueserie(0)
- print "La fonction cout J vaut ...: ",D.valueserie(0)
-
- if (abs(D.valueserie(0) - 16.5) > 1.e-6) :
- raise ValueError("The computation of the cost function is NOT correct")
- else :
- print "The computation of the cost function is OK"
- print
- #
- # float simple
- # ------------
- x = 1.
- print type(x)
- xb = 2.
- yo = 5.
- H = numpy.matrix(numpy.identity(1))
- Hx = numpy.dot(H,x)
- Hx = Hx.T
- B = 1.
- R = 1.
- #
- D.calculate( x = x, Hx = Hx, xb = xb, yo = yo, R = R, B = B)
- print "Le vecteur x choisi est...:", x
- print "L ebauche xb choisie est...:", xb
- print "Le vecteur d observation est...:", yo
- print "B = ", B
- print "R = ", R
- print "La fonction cout J vaut ...: %.2e"%D.valueserie(1)
- if (abs(D.valueserie(1) - 8.5) > 1.e-6) :
- raise ValueError("The computation of the cost function is NOT correct")
- else :
- print "The computation of the cost function is OK"
- print
-
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Calcul de la variance d'un vecteur à chaque pas. Ce diagnostic très simple
- est présent pour rappeller à l'utilisateur de l'assimilation qu'il faut
- qu'il vérifie les variances de ses écarts en particulier.
-"""
-__author__ = "Jean-Philippe ARGAUD - Septembre 2008"
-
-import numpy
-from daCore import BasicObjects, Persistence
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- BasicObjects.Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = float)
-
- def _formula(self, V):
- """
- Calcul de la variance du vecteur en argument. Elle est faite avec une
- division par la taille du vecteur.
- """
- variance = V.var()
- #
- return variance
-
- def calculate(self, vector = None, step = None):
- """
- Teste les arguments, active la formule de calcul et stocke le résultat
- """
- if vector is None:
- raise ValueError("One vector must be given to compute biais")
- V = numpy.array(vector)
- if V.size < 1:
- raise ValueError("The given vector must not be empty")
- #
- value = self._formula( V)
- #
- self.store( value = value, step = step )
-
-#===============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
- #
- D = ElementaryDiagnostic("Ma variance")
- #
- # Vecteur de type matrix
- # ----------------------
- x = numpy.matrix(([3., 4., 5.]))
- print " Le vecteur de type 'matrix' choisi est..:", x
- print " Le moyenne de ce vecteur est............:", x.mean()
- print " La variance attendue de ce vecteur est..:", x.var()
- #
- D.calculate( vector = x)
- print " La variance obtenue de ce vecteur est...:", D.valueserie(0)
- print
- #
- # Vecteur de type array
- # ---------------------
- x = numpy.array(range(11))
- print " Le vecteur de type 'array' choisi est...:", x
- print " Le moyenne de ce vecteur est............:", x.mean()
- print " La variance attendue de ce vecteur est..:", x.var()
- #
- D.calculate( vector = x)
- print " La variance obtenue de ce vecteur est...:", D.valueserie(1)
- print
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic qui effectue le test du Khi2 pour juger de l'adéquation entre
- la distribution d'un échantillon et une distribution gaussienne dont la
- moyenne et l'écart-type sont calculés sur l'échantillon.
- En input : la tolerance(tolerance) et le nombre de classes(nbclasse)
- En output : Le resultat du diagnostic est une reponse booleenne au test :
- True si l adequation a une distribution gaussienne est valide
- au sens du test du Khi2,
- False dans le cas contraire.
-"""
-__author__ = "Sophie RICCI - Juillet 2008"
-
-import numpy
-from daCore import BasicObjects, Persistence
-from ComputeKhi2 import ComputeKhi2_Gauss
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- """
- """
- def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- BasicObjects.Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
- for key in ["tolerance", "dxclasse", "nbclasses"]:
- if not self.parameters.has_key(key):
- raise ValueError("A parameter named \"%s\" is required."%key)
-
- def formula(self, V):
- """
- Effectue le calcul de la p-value pour un vecteur et une distribution
- gaussienne et un nombre de classes donne en parametre du diagnostic.
- """
-
- [vectclasse, eftho, efobs, valeurKhi2, areaKhi2, message] = ComputeKhi2_Gauss(
- vectorV = V,
- dx = self.parameters["dxclasse"],
- nbclasses = self.parameters["nbclasses"],
- SuppressEmptyClasses = True)
-
- logging.info( message )
- logging.info( "(si <%.2f %s on refuse effectivement l'adéquation)"%(100.*self.parameters["tolerance"],"%") )
- logging.info("vecteur des classes=%s"%numpy.size(vectclasse) )
- logging.info("valeurKhi2=%s"%valeurKhi2)
- logging.info("areaKhi2=%s"%areaKhi2)
- logging.info("tolerance=%s"%self.parameters["tolerance"])
-
- if (areaKhi2 < (100.*self.parameters["tolerance"])) :
- answerKhisquareTest = False
- else:
- answerKhisquareTest = True
- logging.info( "La réponse au test est donc est %s"%answerKhisquareTest )
- return answerKhisquareTest
-
- def calculate(self, vector = None, step = None):
- """
- Active la formule de calcul
- """
- if vector is None:
- raise ValueError("One vector must be given to calculate the Khi2 test")
- V = numpy.array(vector)
- if V.size < 1:
- raise ValueError("The given vector must not be empty")
- #
- value = self.formula( V )
- #
- self.store( value = value, step = step)
-
-# ==============================================================================
-if __name__ == "__main__":
- print "\n AUTODIAGNOSTIC \n"
-
- print " Test d adequation du khi-2 a une gaussienne pour un vecteur x"
- print " connu de taille 1000, issu d'une distribution gaussienne normale"
- print " en fixant la largeur des classes"
- print
- #
- # Initialisation des inputs et appel du diagnostic
- # ------------------------------------------------
- tolerance = 0.05
- dxclasse = 0.1
- D = ElementaryDiagnostic("AdequationGaussKhi2", parameters = {
- "tolerance":tolerance,
- "dxclasse":dxclasse,
- "nbclasses":None,
- })
- #
- # Tirage de l'echantillon aleatoire
- # ---------------------------------
- numpy.random.seed(2490)
- x = numpy.random.normal(50.,1.5,1000)
- #
- # Calcul
- # ------
- D.calculate(x)
- #
- if D.valueserie(0) :
- print " L'adequation a une distribution gaussienne est valide."
- print
- else :
- raise ValueError("L'adéquation a une distribution gaussienne n'est pas valide.")
-
-
- print " Test d adequation du khi-2 a une gaussienne pour u:n vecteur x"
- print " connu de taille 1000, issu d'une distribution gaussienne normale"
- print " en fixant le nombre de classes"
- print
- #
- # Initialisation des inputs et appel du diagnostic
- # ------------------------------------------------
- tolerance = 0.05
- nbclasses = 70.
- D = ElementaryDiagnostic("AdequationGaussKhi2", parameters = {
- "tolerance":tolerance,
- "dxclasse":None,
- "nbclasses":nbclasses
- })
- #
- # Tirage de l'echantillon aleatoire
- # ---------------------------------
- numpy.random.seed(2490)
- x = numpy.random.normal(50.,1.5,1000)
- #
- # Calcul
- # ------
- D.calculate(x)
- #
- if D.valueserie(0) :
- print " L'adequation a une distribution gaussienne est valide."
- print
- else :
- raise ValueError("L'adequation a une distribution gaussienne n'est pas valide.")
-
-
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnotic de test sur la validité de l'hypothèse de linéarité de l'opérateur
- H entre xp et xm
-
- Pour calculer Hlin on utilise un schéma différences finies centrées 2
- points. On définit un dxparam tel que :
- xp = xb + dxparam
- et
- xm = xb - dxparam
- On calcule Hxp et Hxm pour obtenir Hlin. Hlin est utilise dans le Blue pour
- caler un paramêtre. La question importante est de choisir un dxparam pas
- trop grand.
-
- On veut vérifier ici que l'hypothèse de linéarite du modèle par rapport au
- paramêtre est valide sur l'intervalle du paramêtre [xm, xp]. Pour cela on
- s'assure que l'on peut retrouver la valeur Hxb par les développemenents de
- Taylor en xp et xm. Ainsi on calcule 2 estimations de Hxb, l'une à partir de
- Hxp (notee Hx1) et l'autre à partir de Hxm (notee Hx2), que l'on compare à
- la valeur calculée de Hxb. On s'intèresse ensuite a la distance entre Hxb et
- ses estimés Hx1 et Hx2. Si la distance est inférieure a un seuil de
- tolerance, l hypothese est valide.
-"""
-__author__ = "Sophie RICCI - Septembre 2008"
-
-import numpy
-from daCore import AssimilationStudy, BasicObjects, Persistence
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- BasicObjects.Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
- if not self.parameters.has_key("tolerance"):
- raise ValueError("A parameter named \"tolerance\" is required.")
-
- def formula(self, H, dxparam, Hxp, Hxm, Hx):
- """
- Test sur la validite de l hypothese de linearite de H entre xp et xm
- """
- dimension = numpy.size(Hx)
- #
- # Reconstruit les valeurs Hx1 et Hx2 de Hx a partir de Hxm et Hxp
- # ---------------------------------------------------------------
- Hx1 = Hxm + H.T * dxparam
- Hx2 = Hxp - H.T * dxparam
- #
- # Calcul de l'ecart entre Hx1 et Hx et entre Hx2 et Hx
- # ----------------------------------------------------
- ADD = AssimilationStudy.AssimilationStudy()
- ADD.setDiagnostic("RMS",
- name = "Calcul de la RMS entre Hx1 et Hx et entre Hx2 et Hx")
- RMS = ADD.get("Calcul de la RMS entre Hx1 et Hx et entre Hx2 et Hx")
- RMS.calculate(Hx1,Hx)
- std1 = RMS.valueserie(0)
- RMS.calculate(Hx2,Hx)
- std2 = RMS.valueserie(1)
- #
- # Normalisation des écarts par Hx pour comparer a un pourcentage
- # --------------------------------------------------------------
- RMS.calculate(Hx,Hx-Hx)
- std = RMS.valueserie(2)
- err1=std1/std
- err2=std2/std
- #
- # Comparaison
- # -----------
- if ( (err1 < self.parameters["tolerance"]) and (err2 < self.parameters["tolerance"]) ):
- reponse = True
- else:
- reponse = False
- return reponse
-
- def calculate(self, Hlin = None, deltaparam = None, Hxp = None, Hxm = None, Hx = None, step = None):
- """
- Arguments :
- - Hlin : Operateur d obsevation lineaire
- - deltaparam : pas sur le parametre param
- - Hxp : calcul en xp = xb + deltaparam
- - Hxm : calcul en xm = xb - deltaparam
- - Hx : calcul en x (generalement xb)
- """
- value = self.formula( Hlin, deltaparam, Hxp, Hxm, Hx )
- #
- self.store( value = value, step = step)
-
-#===============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
-
- print " Diagnotic de test sur la validité de l'hypothèse de linéarité de"
- print " l'opérateur H entre xp et xm."
- print
- #
- dimension = 3
- #
- # Définition des données
- # ----------------------
- Hx = numpy.array(([ 2., 4., 6.]))
- Hxp = numpy.array(([ 3., 5., 7.]))
- Hxm = numpy.array(([ 1., 3., 5.]))
- H = (Hxp - Hxm)/(2.)
- dxparam = 1.
- #
- # Instanciation de l'objet diagnostic
- # -----------------------------------
- D = ElementaryDiagnostic("Mon TestHlin", parameters = {"tolerance": 0.1})
- #
- # Calcul
- # ------
- D.calculate( Hlin = H, deltaparam = dxparam, Hxp = Hxp, Hxm = Hxm, Hx = Hx)
-
- # Validation du calcul
- # --------------------
- if not D.valueserie(0) :
- raise ValueError("La linearisation de H autour de x entre xm et xp est fausse pour ce cas test lineaire")
- else :
- print " La linéarisation de H autour de x entre xm et xp est valide pour ce cas-test linéaire."
- print
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic qui effectue le test du Khi2 pour juger de l'homogénéite entre
- les distributions de 2 vecteurs quelconques.
- - entrée : la tolerance (tolerance) et le nombre de classes (nbclasse),
- sous forme de paramètres dans le dictionnaire Par
- - sortie : le resultat du diagnostic est une reponse booleenne au test :
- True si l homogeneite est valide au sens du test du Khi2,
- False dans le cas contraire.
-"""
-__author__ = "Sophie RICCI - Juillet 2008"
-
-import numpy
-from daCore import BasicObjects, Persistence
-from ComputeKhi2 import ComputeKhi2_Homogen
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- def __init__(self, name="", unit="", basetype = None, parameters = {} ):
- BasicObjects.Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool )
- for key in ["tolerance", "dxclasse", "nbclasses"]:
- if not self.parameters.has_key(key):
- raise ValueError("A parameter named \"%s\" is required."%key)
-
- def _formula(self, V1, V2):
- """
- Effectue le calcul de la p-value pour deux vecteurs et un nombre de
- classes donne en parametre du diagnostic.
- """
- [classes, eftheo, efobs, valeurKhi2, areaKhi2, message] = ComputeKhi2_Homogen(
- vectorV1 = V1,
- vectorV2 = V2,
- dx = self.parameters["dxclasse"],
- nbclasses = self.parameters["nbclasses"],
- SuppressEmptyClasses = True)
- #
- logging.info( message )
- logging.info( "(si <%.2f %s on refuse effectivement l'homogeneite)"%(100.*self.parameters["tolerance"],"%") )
- #
- answerKhisquareTest = False
- if (areaKhi2 < (100.*self.parameters["tolerance"])) :
- answerKhisquareTest = False
- else:
- answerKhisquareTest = True
- #
- return answerKhisquareTest
-
- def calculate(self, vector1 = None, vector2 = None, step = None):
- """
- Active la formule de calcul
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Khi2 value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- if V1.size != V2.size:
- raise ValueError("The two given vectors must have the same size")
- #
- value = self._formula( V1, V2 )
- #
- self.store( value = value, step = step )
-
-# ==============================================================================
-if __name__ == "__main__":
- print "\n AUTODIAGNOSTIC \n"
-
- print " Test d'homogeneite du Khi-2 pour deux vecteurs de taille 10,"
- print " issus d'une distribution gaussienne normale"
- print
- #
- # Initialisation des inputs et appel du diagnostic
- # --------------------------------------------------------------------
- tolerance = 0.05
- dxclasse = 0.5
- D = ElementaryDiagnostic("HomogeneiteKhi2", parameters = {
- "tolerance":tolerance,
- "dxclasse":dxclasse,
- "nbclasses":None,
- })
- #
- # Tirage de l'echantillon aleatoire
- # --------------------------------------------------------------------
- numpy.random.seed(4000)
- x1 = numpy.random.normal(50.,1.5,10000)
- numpy.random.seed(2490)
- x2 = numpy.random.normal(50.,1.5,10000)
- #
- # Calcul
- # --------------------------------------------------------------------
- D.calculate(x1, x2)
- #
- print " La reponse du test est \"%s\" pour une tolerance de %.2e et une largeur de classe de %.2e "%(D.valueserie(0), tolerance, dxclasse)
- print
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic sur la reduction du biais lors de l'analyse
-"""
-__author__ = "Sophie RICCI - Aout 2008"
-
-import numpy
-from daCore import BasicObjects, Persistence
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- BasicObjects.Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool)
-
- def _formula(self, V1, V2):
- """
- Vérification de la reduction du biais entre OMB et OMA lors de l'analyse
- """
- biaisOMB = V1.mean()
- biaisOMA = V2.mean()
- #
- if biaisOMA > biaisOMB:
- reducebiais = False
- else :
- reducebiais = True
- #
- return reducebiais
-
- def calculate(self, vectorOMB = None, vectorOMA = None, step = None):
- """
- Teste les arguments, active la formule de calcul et stocke le résultat
- Arguments :
- - vectorOMB : vecteur d'écart entre les observations et l'ébauche
- - vectorOMA : vecteur d'écart entre les observations et l'analyse
- """
- if ( (vectorOMB is None) or (vectorOMA is None) ):
- raise ValueError("Two vectors must be given to test the reduction of the biais after analysis")
- V1 = numpy.array(vectorOMB)
- V2 = numpy.array(vectorOMA)
- if V1.size < 1 or V2.size < 1:
- raise ValueError("The given vectors must not be empty")
- if V1.size != V2.size:
- raise ValueError("The two given vectors must have the same size")
- #
- value = self._formula( V1, V2 )
- #
- self.store( value = value, step = step )
-
-#===============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
- #
- # Instanciation de l'objet diagnostic
- # -----------------------------------
- D = ElementaryDiagnostic("Mon ReduceBiais")
- #
- # Tirage des 2 vecteurs choisis
- # -------------------------------
- x1 = numpy.matrix(([3. , 4., 5. ]))
- x2 = numpy.matrix(([1.5, 2., 2.5]))
- print " L'écart entre les observations et l'ébauche est OMB :", x1
- print " La moyenne de OMB (i.e. le biais) est de............:", x1.mean()
- print " L'écart entre les observations et l'analyse est OMA :", x2
- print " La moyenne de OMA (i.e. le biais) est de............:", x2.mean()
- #
- D.calculate( vectorOMB = x1, vectorOMA = x2)
- if not D.valueserie(0) :
- print " Résultat : l'analyse NE RÉDUIT PAS le biais"
- else :
- print " Résultat : l'analyse RÉDUIT le biais"
- print
- #
- # Tirage des 2 vecteurs choisis
- # -------------------------------
- x1 = numpy.matrix(range(-5,6))
- x2 = numpy.array(range(11))
- print " L'écart entre les observations et l'ébauche est OMB :", x1
- print " La moyenne de OMB (i.e. le biais) est de............:", x1.mean()
- print " L'écart entre les observations et l'analyse est OMA :", x2
- print " La moyenne de OMA (i.e. le biais) est de............:", x2.mean()
- #
- D.calculate( vectorOMB = x1, vectorOMA = x2)
- if not D.valueserie(1) :
- print " Résultat : l'analyse NE RÉDUIT PAS le biais"
- else :
- print " Résultat : l'analyse RÉDUIT le biais"
- print
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Diagnostic sur les variances dans B et R par rapport à l'ébauche Xb et aux
- observations Y. On teste si on a les conditions :
- 1%*xb < sigma_b < 10%*xb
- et
- 1%*yo < sigma_o < 10%*yo
- Le diagnostic renvoie True si les deux conditions sont simultanément
- vérifiées, False dans les autres cas.
-"""
-__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Septembre 2008"
-
-import numpy
-from scipy.linalg import eig
-from daCore import BasicObjects, Persistence
-import logging
-
-# ==============================================================================
-class ElementaryDiagnostic(BasicObjects.Diagnostic,Persistence.OneScalar):
- def __init__(self, name = "", unit = "", basetype = None, parameters = {}):
- BasicObjects.Diagnostic.__init__(self, name, parameters)
- Persistence.OneScalar.__init__( self, name, unit, basetype = bool )
-
- def _formula(self, xb, B, yo, R):
- """
- Comparaison des variables et de leur variance relative
- """
- valpB = eig(B, left = False, right = False)
- valpR = eig(R, left = False, right = False)
- logging.info(" Si l on souhaite 1%s*xb < sigma_b < 10%s*xb, les valeurs propres de B doivent etre comprises dans l intervalle [%.3e,%.3e]"%("%","%",1.e-4*xb.mean()*xb.mean(),1.e-2*xb.mean()*xb.mean()))
- logging.info(" Si l on souhaite 1%s*yo < sigma_o < 10%s*yo, les valeurs propres de R doivent etre comprises dans l intervalle [%.3e,%.3e]"%("%","%",1.e-4*yo.mean()*yo.mean(),1.e-2*yo.mean()*yo.mean()))
- #
- limite_inf_valp = 1.e-4*xb.mean()*xb.mean()
- limite_sup_valp = 1.e-2*xb.mean()*xb.mean()
- variancexb = (valpB >= limite_inf_valp).all() and (valpB <= limite_sup_valp).all()
- logging.info(" La condition empirique sur la variance de Xb est....: %s"%variancexb)
- #
- limite_inf_valp = 1.e-4*yo.mean()*yo.mean()
- limite_sup_valp = 1.e-2*yo.mean()*yo.mean()
- varianceyo = (valpR >= limite_inf_valp).all() and (valpR <= limite_sup_valp).all()
- logging.info(" La condition empirique sur la variance de Y est.....: %s",varianceyo)
- #
- variance = variancexb and varianceyo
- logging.info(" La condition empirique sur la variance globale est..: %s"%variance)
- #
- return variance
-
- def calculate(self, Xb = None, B = None, Y = None, R = None, step = None):
- """
- Teste les arguments, active la formule de calcul et stocke le résultat
- Arguments :
- - Xb : valeur d'ébauche du paramêtre
- - B : matrice de covariances d'erreur d'ébauche
- - yo : vecteur d'observation
- - R : matrice de covariances d'erreur d'observation
- """
- if (Xb is None) or (B is None) or (Y is None) or (R is None):
- raise ValueError("You must specify Xb, B, Y, R")
- yo = numpy.array(Y)
- BB = numpy.matrix(B)
- xb = numpy.array(Xb)
- RR = numpy.matrix(R)
- if (RR.size < 1 ) or (BB.size < 1) :
- raise ValueError("The background and the observation covariance matrices must not be empty")
- if ( yo.size < 1 ) or ( xb.size < 1 ):
- raise ValueError("The Xb background and the Y observation vectors must not be empty")
- if xb.size*xb.size != BB.size:
- raise ValueError("Xb background vector and B covariance matrix sizes are not consistent")
- if yo.size*yo.size != RR.size:
- raise ValueError("Y observation vector and R covariance matrix sizes are not consistent")
- if yo.all() == 0. or xb.all() == 0. :
- raise ValueError("The diagnostic can not be applied to zero vectors")
- #
- value = self._formula( xb, BB, yo, RR)
- #
- self.store( value = value, step = step )
-
-#===============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
- #
- # Instanciation de l'objet diagnostic
- # -----------------------------------
- D = ElementaryDiagnostic("Mon OrdreVariance")
- #
- # Vecteur de type matrix
- # ----------------------
- xb = numpy.array([11000.])
- yo = numpy.array([1.e12 , 2.e12, 3.e12 ])
- B = 1.e06 * numpy.matrix(numpy.identity(1))
- R = 1.e22 * numpy.matrix(numpy.identity(3))
- #
- D.calculate( Xb = xb, B = B, Y = yo, R = R)
- print " L'ébauche est.......................................:",xb
- print " Les observations sont...............................:",yo
- print " La valeur moyenne des observations est..............: %.2e"%yo.mean()
- print " La valeur moyenne de l'ebauche est..................: %.2e"%xb.mean()
- print " La variance d'ébauche specifiée est.................: %.2e"%1.e6
- print " La variance d'observation spécifiée est.............: %.2e"%1.e22
- #
- if D.valueserie(0) :
- print " Les variances specifiées sont de l'ordre de 1% a 10% de l'ébauche et des observations"
- else :
- print " Les variances specifiées ne sont pas de l'ordre de 1% a 10% de l'ébauche et des observations"
- print
-
-
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Outil numerique de calcul de la variable Khi2
-
- On peut realiser deux types de test du Khi2 :
- - test d'adequation : comparer la distribution d'un echantillon a une
- distribution theorique,
- - test d'homogeneite : comparer les distributions de 2 vecteurs.
-
- Pour le test d'adequation, on travaille sur une gaussienne
- dont la moyenne et l'ecart type sont calcules sur
- l'echantillon, soit donnes.
-
- Ce fichier contient une classe "StatspourTests" de methodes qui realisent
- differentes etapes utiles aux calculs des tests du Khi2.
-
- Ce fichier contient de plus 3 methodes : ComputeKhi2_testGauss,
- ComputeKhi2_Gauss et ComputeKhi2_Homogen.
- - ComputeKhi2_testGauss : calcul la distance du Khi2 entre un vecteur
- aleatoire issu d un gaussienne et une distribution theorique gaussienne
- dont on specifie la moyenne et l ecart type
- - ComputeKhi2_Gauss : calcul la distance du Khi2 entre un vecteur donne et
- une distribution theorique gaussienne dont la moyenne et l ecart type sont
- calcules sur l echantillon
- - ComputeKhi2_Homogen : calcul la distance du Khi2 entre deux vecteurs donnes
-
- Ces methodes necessitent et fournissent :
- - en input :
- - le ou les vecteurs dont on etudie la distribution,
- - la distribution theorique et eventuellement la moyenne et ecart type,
- - la largeur des classes,
- - un booleen traduisant la suppression des classes vides
- - en output :
- - le vecteur des classes,
- - les pdf theorique et donnee,
- - la valeur du Khi2,
- - la p-value qui represent l'aire de la queue de la distribution du
- Khi2 et
- - le message qui interprete le test.
-"""
-__author__ = "Sophie RICCI - Mars 2010"
-
-import numpy
-from numpy import random
-from scipy import arange, asarray, stats
-from scipy.stats import histogram2, chisquare, chisqprob, norm
-from scipy.stats import distributions
-import scipy, sys, time, os
-import logging
-
-# ==============================================================================
-class StatspourTests :
- """
- Classe de methodes pour la preparation du test de Khi2
- """
- 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) :
-
-
- if (pdftest is None and obs is None) :
- raise ValueError('Donner soit une pdf de test soit un vecteur obs')
- if not obs is None :
- if pdftest is None :
- self.__obs=asarray(obs)
- if not pdftest is None :
- if obs is None :
- if len(pdftest) == 3:
- niter=eval(pdftest[2])
- obs=[eval(" ".join(pdftest[:2])) for z in range(niter)]
- self.__obs=asarray(obs)
- else :
- self.__obs=asarray(eval(" ".join(pdftest[:2])))
- if not (obsHomogen is None) :
- self.__obsHomogen = asarray(obsHomogen)
- self.__testHomogen = True
- else :
- self.__testHomogen = False
-
-
- self.__mean_exp = self.__obs.mean()
- self.__std_exp = self.__obs.std()
-
- if cdftheo is None : raiseValueError(" ... Definir le parametre cdftheo ...")
- if use_mean_std_exp :
- self.__cdf=cdftheo( self.__mean_exp, self.__std_exp).cdf
- else :
- self.__cdf=cdftheo( meantheo, stdtheo).cdf
-
- self.__min=min(self.__obs)
- self.__max=max(self.__obs)
- self.__N=len(self.__obs)
- self.__use_mean_std_exp=use_mean_std_exp
- self.__dxmin=dxmin
- self.__nbclasses = nbclasses
- if not (dxmin is None) and not (nbclasses is None) :
- raise ValueError("... Specifier soit le nombre de classes, soit la largeur des classes")
- if (dxmin is None) and (nbclasses is None) :
- raise ValueError("... Specifier soit le nombre de classes, soit la largeur des classes")
- if not (nbclasses is None) and (dxmin is None) :
- self.__dxmin = (self.__max - self.__min ) / float(self.__nbclasses)
- return None
-
- def MakeClasses(self) :
- """
- Classification en classes
- """
- self.__subdiv=arange(self.__min,self.__max+self.__dxmin,self.__dxmin)
- self.__modalites=len(self.__subdiv)
- return None
-
- def ComputeObs(self):
- """
- Calcul de la probabilite observee de chaque classe
- """
- self.__kobs=histogram2(self.__obs,self.__subdiv)[1:]
- return self.__kobs
-
- def ComputeObsHomogen(self):
- """
- Calcul de la probabilite observee pour le test homogeneite de chaque classe
- """
- self.__kobsHomogen=histogram2(self.__obsHomogen,self.__subdiv)[1:]
- return self.__kobsHomogen
-
- def ComputeTheo(self):
- """
- Calcul de la probabilite theorique de chaque classe
- """
- self.__ktheo=[self.__cdf(self.__subdiv[i+1])-self.__cdf(self.__subdiv[i]) for i in range(self.__modalites-1)]
- self.__ktheo=asarray(self.__ktheo)
- self.__ktheo=(sum(self.__kobs)/sum(self.__ktheo))*self.__ktheo
-
- def Computepdfs(self) :
-
- self.__subdiv=self.__subdiv[1:]
- self.__pdfobs=[self.__kobs[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)]
-
- if self.__testHomogen :
- self.__pdftheo=[self.__kobsHomogen[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)]
- else :
- self.__pdftheo=[self.__ktheo[i+1]/(self.__subdiv[i+1]-self.__subdiv[i]) for i in range(self.__modalites-2)]
-
- return self.__subdiv, self.__pdftheo, self.__pdfobs
-
- def Computeddl(self):
- """
- Calcul du nombre de degres de liberte
- """
- self.__ddl = self.__modalites - 1.
- if self.__use_mean_std_exp :
- self.__ddl = self.__ddl - 2.
- if (self.__ddl < 1.):
- raise ValueError("The ddl is 0, you must increase the number of classes nbclasse ")
- logging.debug("Nombre de degres de liberte=%s"%self.__ddl)
-
- def ComputeValue(self) :
- """
- Calcul de la variable Q qui suit une loi Khi-2
- """
- if self.__testHomogen :
- kobs,ktheo=self.__kobs.tolist(),self.__kobsHomogen.tolist()
- else :
- kobs,ktheo=self.__kobs.tolist(),self.__ktheo.tolist()
-
- # on supprime les classes theoriques qui ont moins d'un element (sinon la distance khi2 tendrait vers l'infini)
- ko,kt=[],[]
- self.__count0 = 0.
- for k,val in enumerate(ktheo):
- if val > 1.0:
- kt.append(val)
- ko.append(kobs[k])
- else :
- self.__count0 = self.__count0 +1.
- logging.debug("WARNING : nombre de classes vides supprimees (effectif theorique inferieur a 1.) pour le calcul de la valeur du Khi2 = %s"%self.__count0)
- ef1,ef2=asarray(ko),asarray(kt)
- count = 0.
- for el in ef1.tolist() :
- if el < 5. :
- count = count +1.
- for el in ef2.tolist() :
- if el < 5. :
- count = count +1.
- pourcent_nbclasse_effecinf = count /(2.*len(ef1.tolist())) *100.
- if (pourcent_nbclasse_effecinf > 20.) :
- logging.debug("WARNING : nombre de classes dont l effectif est inferieur a 5 elements %s"%pourcent_nbclasse_effecinf)
- k,p = chisquare(ef1, ef2)
- k2, p2 = [k],[p]
- for shift in range(1,6):
- k,p=chisquare(ef1[shift:],ef2[:-shift])
- k2.append(k)
- p2.append(p)
- k,p=chisquare(ef1[:-shift],ef2[shift:])
- k2.append(k)
- p2.append(p)
- logging.debug("Liste des valeurs du Khi2 = %s"%k2)
- self.__khi2=min(k2)
- self.__Q=self.__khi2
-
- logging.debug("Valeur du Khi2=%s"%self.__Q)
- return self.__Q
-
- def ComputeValue_ks(self) :
- """
- Calcul de la valeur ks pour le test de kolmogorov smirnov
- """
- if self.__testHomogen :
- kobs,ktheo=self.__kobs.tolist(),self.__kobsHomogen.tolist()
- else :
- kobs,ktheo=self.__kobs.tolist(),self.__ktheo.tolist()
-
- self.sortobs=asarray([scipy.sum(kobs[:i]) for i in range(1,1+(scipy.shape(kobs)[0]))])
- self.sortobs=self.sortobs/float(scipy.sum(kobs))
- self.sortthe=self.__cdf(self.__subdiv)
- DKS1=max(abs(self.sortthe-self.sortobs))
- DKS2=max(abs(self.sortthe[1:]-self.sortobs[0:-1]))
- self.__ks=max(DKS1,DKS2)
- logging.debug("Valeur du ks test=%s"%self.__ks)
- return self.__ks
-
- def ComputeArea(self):
- """
- Calcul de la p-value
- """
- self.__areakhi2 = 100 * chisqprob(self.__Q, self.__ddl)
- return self.__areakhi2
-
- def ComputeArea_ks(self):
- """
- Calcul de la p-value ks
- """
- self.__areaks = 100 * distributions.ksone.sf(self.__ks,self.__modalites)
- return self.__areaks
-
- def WriteMessage(self):
- """
- Interpretation du test
- """
- message = "Il y a %.2f%s de chance de se tromper en refusant l'adequation"%(self.__areakhi2,"%")
- return message
-
- def WriteMessage_ks(self):
- """
- Interpretation du test
- """
- message = "Il y a %.2f%s de chance de se tromper en refusant l'adequation"%(self.__areaks,"%")
- return message
-
- def WriteMessageHomogen(self):
- """
- Interpretation du test
- """
- message = "Il y a %.2f%s de chance de se tromper en refusant l'homogeneite"%(self.__areakhi2,"%")
- return message
-
-# ==============================================================================
-def ComputeKhi2_testGauss(
- meantheo = 0.,
- stdtheo = 1.,
- nech = 10,
- dx = 0.1,
- nbclasses = None,
- SuppressEmptyClasses = True,
- ):
- """
- Test du Khi2 d adequation entre tirage aleatoire dans gaussienne et une gaussienne theo
- """
- 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)
- essai.MakeClasses()
- essai.ComputeObs()
- essai.ComputeTheo()
- classes,eftheo, efobs = essai.Computepdfs()
- essai.Computeddl()
- valeurKhi2= essai.ComputeValue()
- areaKhi2 = essai.ComputeArea()
- message = essai.WriteMessage()
- logging.debug("message %s"%message)
- return classes, eftheo, efobs, valeurKhi2, areaKhi2, message
-
-def ComputeKhi2_Gauss(
- vectorV = None,
- dx = 0.1,
- SuppressEmptyClasses = True,
- nbclasses = None
- ):
- """
- Test du Khi2 d adequation entre un vecteur donne et une gaussienne theo de mean et std celles du vecteur
- """
- essai = StatspourTests( cdftheo=norm, pdftest = None, obs = vectorV, use_mean_std_exp=True,dxmin=dx, obsHomogen = None, nbclasses = nbclasses)
- essai.MakeClasses()
- essai.ComputeObs()
- essai.ComputeTheo()
- classes,eftheo, efobs = essai.Computepdfs()
- essai.Computeddl()
- valeurKhi2= essai.ComputeValue()
- areaKhi2 = essai.ComputeArea()
- message = essai.WriteMessage()
- logging.debug("message %s"%message)
- return classes, eftheo, efobs, valeurKhi2, areaKhi2, message
-
-def Compute_ks(
- vectorV = None,
- dx = 0.1,
- SuppressEmptyClasses = True,
- nbclasses = None
- ):
- """
- Test du KS adequation entre un vecteur donne et une gaussienne theo de mean et std celles du vecteur
- """
- essai = StatspourTests( cdftheo=norm, pdftest = None, obs = vectorV, use_mean_std_exp=True,dxmin=dx, obsHomogen = None, nbclasses = nbclasses)
- essai.MakeClasses()
- essai.ComputeObs()
- essai.ComputeTheo()
- classes,eftheo, efobs = essai.Computepdfs()
- essai.Computeddl()
- valeur_ks= essai.ComputeValue_ks()
- area_ks = essai.ComputeArea_ks()
- message = essai.WriteMessage_ks()
- logging.debug("message %s"%message)
- return classes, eftheo, efobs, valeur_ks, area_ks, message
-
-
-def ComputeKhi2_Homogen(
- vectorV1 = None,
- vectorV2 = None,
- dx = 0.1,
- SuppressEmptyClasses = True,
- nbclasses = None
- ):
- """
- Test du Khi2 d homogeniete entre 2 vecteurs
- """
- essai = StatspourTests( cdftheo=norm, pdftest = None, obs = vectorV1, use_mean_std_exp=True,dxmin=dx, obsHomogen = vectorV2, nbclasses = nbclasses)
- essai.MakeClasses()
- essai.ComputeObs()
- essai.ComputeObsHomogen()
- classes,eftheo, efobs = essai.Computepdfs()
- essai.Computeddl()
- valeurKhi2= essai.ComputeValue()
- areaKhi2 = essai.ComputeArea()
- message = essai.WriteMessageHomogen()
- logging.debug("message %s"%message)
- return classes, eftheo, efobs, valeurKhi2, areaKhi2, message
-
-# ==============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
- #
- numpy.random.seed(100)
-#-------------------------------------------------------------------------
- # Test de verification d adequation entre une gaussienne et un tirage gaussien
- print ''
- print 'Test de verification d adequation entre une gaussienne centree normale et un tirage gaussien'
- classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_testGauss(meantheo = 0., stdtheo = 1., nech = 1000., dx = 0.1, SuppressEmptyClasses = True, nbclasses = None)
- print ' valeurKhi2=',valeurKhi2
- print ' areaKhi2=',areaKhi2
- print ' ',message
-
- if (numpy.abs(areaKhi2 - 99.91)< 1.e-2) :
- print "The computation of the khisquare value is OK"
- else :
- raise ValueError("The computation of the khisquare value is WRONG")
-
- numpy.random.seed(2490)
-#-------------------------------------------------------------------------
- # Test de verification d adequation entre une gaussienne et un vecteur donne
- print ''
- print 'Test de verification d adequation entre une gaussienne et un vecteur donne'
- V = random.normal(50.,1.5,1000)
- classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = 0.1, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
- print ' valeurKhi2=',valeurKhi2
- print ' areaKhi2=',areaKhi2
- print ' ',message
-
- if (numpy.abs(areaKhi2 - 99.60)< 1.e-2) :
- print "The computation of the khisquare value is OK"
- else :
- raise ValueError("The computation of the khisquare value is WRONG")
-#-------------------------------------------------------------------------
- # Test de verification d adequation entre une gaussienne et un vecteur donne par KS
- print ''
- print 'Test de verification d adequation entre une gaussienne et un vecteur donne par KS'
-# V = random.normal(50.,1.5,1000)
- classes, eftheo, efobs, valeur_ks, area_ks, message = Compute_ks(dx = 0.1, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
- print ' valeur_ks=',valeur_ks
- print ' area_ks=',area_ks
- print ' ',message
-
- if (numpy.abs(area_ks - 82.17)< 1.e-2) :
- print "The computation of the ks value is OK"
- else :
- raise ValueError("The computation of the ks value is WRONG")
-
-#-------------------------------------------------------------------------
- # Test de d homogeneite entre 2 vecteurs donnes
- print ''
- print 'Test d homogeneite entre 2 vecteurs donnes'
- V1 = random.normal(50.,1.5,10000)
- numpy.random.seed(2490)
- V2 = random.normal(50.,1.5,10000)
- classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Homogen(dx = 0.5, vectorV1 = V1, vectorV2 = V2, SuppressEmptyClasses = True, nbclasses = None)
- print ' valeurKhi2=',valeurKhi2
- print ' areaKhi2=',areaKhi2
- print ' ',message
-
- if (numpy.abs(areaKhi2 - 99.98)< 1.e-2) :
- print "The computation of the khisquare value is OK"
- else :
- raise ValueError("The computation of the khisquare value is WRONG")
-#-------------------------------------------------------------------------
- # Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 10000
- print ''
- print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 10000'
-# file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech10000.gnu'
-# fid = open(file, "w")
-# lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' )
- numpy.random.seed(4000)
- V = random.normal(0., 1.,10000)
- aire = []
- for dx in arange(0.01, 1., 0.001) :
- classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
-# lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2)
- aire.append(areaKhi2)
- meanaire = numpy.asarray(aire)
-# fid.writelines(lines)
-
- 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"
- print
- if (numpy.abs( meanaire.mean() - 71.79)< 1.e-2) :
- print "The computation of the khisquare value is OK"
- else :
- raise ValueError("The computation of the khisquare value is WRONG")
-#-------------------------------------------------------------------------
- # Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 1000
- print ''
- print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 1000'
-# file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech1000.gnu'
-# fid = open(file, "w")
-# lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' )
- numpy.random.seed(4000)
- V = random.normal(0., 1.,1000)
- aire = []
- for dx in arange(0.05, 1., 0.001) :
- classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
-# lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2)
- aire.append(areaKhi2)
- meanaire = numpy.asarray(aire)
-# fid.writelines(lines)
-
- 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"
- print
- if (numpy.abs( meanaire.mean() - 90.60)< 1.e-2) :
- print "The computation of the khisquare value is OK"
- else :
- raise ValueError("The computation of the khisquare value is WRONG")
-#-------------------------------------------------------------------------
-
- # Test de verification d adequation entre une gaussienne et un tirage gaussien en faisant varier le nombre de classes, echantillon de taille 100
- print ''
- print 'Test de verification d adequation entre une gaussienne et un vecteur aleatoire gaussien de taille 100'
-# file = 'ComputeKhi2_adequationGauss_fctnbclasses_nech100.gnu'
-# fid = open(file, "w")
-# lines = '%s\n' % ('# dx , nbclasses, valeurKhi2, ProbKhi2' )
- numpy.random.seed(4000)
- V = random.normal(0., 1.,100)
- aire = []
- for dx in arange(0.1, 1., 0.01) :
- classes, eftheo, efobs, valeurKhi2, areaKhi2, message = ComputeKhi2_Gauss(dx = dx, vectorV = V, SuppressEmptyClasses = True, nbclasses = None)
-# lines += '%f %f %f %f\n' % (dx, numpy.size(classes), valeurKhi2, areaKhi2)
- aire.append(areaKhi2)
- meanaire = numpy.asarray(aire)
-# fid.writelines(lines)
-
- 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"
- print
+++ /dev/null
-#-*-coding:iso-8859-1-*-
-#
-# Copyright (C) 2008-2010 EDF R&D
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-#
-# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-#
-__doc__ = """
- Outil numérique de calcul des variables de Student pour 2 vecteurs
- dépendants ou indépendants, avec variances supposées égales ou différentes
-"""
-__author__ = "Sophie RICCI, Jean-Philippe ARGAUD - Octobre 2008"
-
-import sys ; sys.path.insert(0, "../daCore")
-
-import numpy
-from scipy.stats import ttest_rel, ttest_ind, betai
-import logging
-
-# ==============================================================================
-def DependantVectors(vector1 = None, vector2 = None, tolerance = 0.05 ):
- """
- Outil numérique de calcul de la variable de Student pour 2 vecteurs
- dépendants
- Ce calcul nécessite :
- - en input :
- - les deux vecteurs (comme liste, array ou matrix)
- d'échantillons dont on veut comparer la variance,
- - la tolérance
- - en output :
- - la p-value,
- - la valeur de la variable aléatoire,
- - la reponse au test pour une tolerance ainsi que
- - le message qui interprete la reponse du test.
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Student value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- if V1.size != V2.size:
- raise ValueError("The two given vectors must have the same size, or the vector types are incompatible")
- #
- # Calcul de la p-value du Test de Student
- # --------------------------------------------------------------------
- [t, prob] = ttest_rel(V1, V2)
- areastudent = 100. * prob
- #
- logging.debug("DEPENDANTVECTORS t = %.3f, areastudent = %.3f"%(t, areastudent))
- #
- # Tests
- # --------------------------------------------------------------------
- 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,"%")
- logging.debug(message)
- #
- if (areastudent < (100.*tolerance)) :
- answerTestStudent = False
- else:
- answerTestStudent = True
- #
- return areastudent, t, answerTestStudent, message
-
-# ==============================================================================
-def IndependantVectorsDifferentVariance(vector1 = None, vector2 = None, tolerance = 0.05 ):
- """
- Outil numerique de calcul de la variable de Student pour 2 vecteurs independants supposes de variances vraies differentes
- En input : la tolerance
- 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.
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Student value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- #
- # Calcul de la p-value du Test de Student
- # --------------------------------------------------------------------
- # t = (m1 - m2)/ sqrt[ (var1/n1 + var2/n2) ]
- # ou var est calcule comme var = somme (xi -xmena)**2 /(n-1)
- n1 = V1.size
- n2 = V2.size
- mean1 = V1.mean()
- mean2 = V2.mean()
- var1 = numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std() * numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std()
- var2 = numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std() * numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std()
- t = (mean1 - mean2)/ numpy.sqrt( var1/n1 + var2/n2 )
- df = ( (var1/n1 + var2/n2) * (var1/n1 + var2/n2) ) / ( (var1/n1)*(var1/n1)/(n1-1) + (var2/n2)*(var2/n2)/(n2-1) )
- zerodivproblem = var1/n1 + var2/n2 == 0
- t = numpy.where(zerodivproblem, 1.0, t) # replace NaN t-values with 1.0
- prob = betai(0.5*df,0.5,float(df)/(df+t*t))
- areastudent = 100. * prob
- #
- logging.debug("IndependantVectorsDifferentVariance t = %.3f, areastudent = %.3f"%(t, areastudent))
- #
- # Tests
- # --------------------------------------------------------------------
- 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,"%")
- logging.debug(message)
- if (areastudent < (100.*tolerance)) :
- answerTestStudent = False
- else:
- answerTestStudent = True
- #
- return areastudent, t, answerTestStudent, message
-
-# ==============================================================================
-def IndependantVectorsEqualVariance(vector1 = None, vector2 = None, tolerance = 0.05 ):
- """
- Outil numerique de calcul de la variable de Student pour 2 vecteurs independants supposes de meme variance vraie
- En input : la tolerance
- 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.
- """
- if (vector1 is None) or (vector2 is None) :
- raise ValueError("Two vectors must be given to calculate the Student value")
- V1 = numpy.array(vector1)
- V2 = numpy.array(vector2)
- if (V1.size < 1) or (V2.size < 1):
- raise ValueError("The given vectors must not be empty")
- #
- # Calcul de la p-value du Test de Student
- # --------------------------------------------------------------------
- # t = sqrt(n1+n2-2) * (m1 - m2)/ sqrt[ (1/n1 +1/n2) * ( (n1-1)var1 + (n2-1)var2 )]
- # ou var est calcule comme var = somme (xi -xmena)**2 /(n-1)
- [t, prob] = ttest_ind(V1, V2)
- areastudent = 100. * prob
- #
- logging.debug("IndependantVectorsEqualVariance t = %.3f, areastudent = %.3f"%(t, areastudent))
- # Tests
- # --------------------------------------------------------------------
- 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,"%")
- logging.debug(message)
- if (areastudent < (100.*tolerance)) :
- answerTestStudent = False
- else:
- answerTestStudent = True
-
- return areastudent, t, answerTestStudent, message
-
-# ==============================================================================
-if __name__ == "__main__":
- print '\n AUTODIAGNOSTIC \n'
- # logging.getLogger().setLevel(logging.DEBUG)
-
- print
- print " Test de Student pour des vecteurs dépendants"
- print " --------------------------------------------"
- # Tirage de l'echantillon
- V1 = numpy.matrix(([-1., 0., 4.])).T
- V2 = numpy.matrix(([-2., 0., 8.])).T
- V1 = V1.A1
- V2 = V2.A1
- #
- # Appel de l outil DependantVectors et initialisation des inputs
- [aire, Q, reponse, message] = DependantVectors(
- vector1 = V1,
- vector2 = V2,
- tolerance = 0.05)
- #
- # Verification par les calculs sans les routines de scipy.stats
- # (ref numerical recipes)
- n = V1.size
- df= n -1
- # Les routines de scipy.stats utilisent une variance calculee avec n-1 et non n comme dans std
- # t = (m1 - m2)/ sqrt[(varx1 + varx2 - 2 cov(x1, x2))/n ]
- # ou var est calcule comme var = somme (xi -xmean)**2 /(n-1)
- var1 = numpy.sqrt(n)/numpy.sqrt(n-1)* V1.std() * numpy.sqrt(n)/numpy.sqrt(n-1) * V1.std()
- var2 = numpy.sqrt(n)/numpy.sqrt(n-1)* V2.std() * numpy.sqrt(n)/numpy.sqrt(n-1) * V2.std()
- m1 = V1.mean()
- m2 = V2.mean()
- cov = 0.
- for j in range(0, n) :
- cov = cov + (V1[j] - m1)*(V2[j] - m2)
- cov = cov /df
- sd = numpy.sqrt((var1 + var2 - 2. *cov) / n)
- tverif = (m1 -m2) /sd
- aireverif = 100. * betai(0.5*df,0.5,float(df)/(df+tverif*tverif))
- if (aireverif - aire < 1.e-5) :
- print " Le calcul est conforme à celui de l'algorithme du Numerical Recipes"
- else :
- raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes")
-
- if (numpy.abs(aire - 57.99159)< 1.e-5) :
- print " Le calcul est JUSTE sur cet exemple."
- else :
- raise ValueError("Le calcul est FAUX sur cet exemple.")
-
- print
- print " Test de Student pour des vecteurs independants supposés de même variance"
- print " ------------------------------------------------------------------------"
- # Tirage de l'echantillon
- V1 = numpy.matrix(([-1., 0., 4.])).T
- V2 = numpy.matrix(([-2., 0., 8.])).T
- V1 = V1.A1
- V2 = V2.A1
- #
- # Appel de l outil IndependantVectorsDifferentVariance et initialisation des inputs
- [aire, Q, reponse, message] = IndependantVectorsDifferentVariance(
- vector1 = V1,
- vector2 = V2,
- tolerance = 0.05)
- #
- if (numpy.abs(aire - 78.91339)< 1.e-5) :
- print " Le calcul est JUSTE sur cet exemple."
- else :
- raise ValueError("Le calcul est FAUX sur cet exemple.")
-
- print
- print " Test de Student pour des vecteurs indépendants supposés de même variance"
- print " ------------------------------------------------------------------------"
- # Tirage de l'echantillon
- V1 = numpy.matrix(([-1., 0., 4.])).T
- V2 = numpy.matrix(([-2., 0., 8.])).T
- V1 = V1.A1
- V2 = V2.A1
- #
- # Appel de l outil IndependantVectorsEqualVariance et initialisation des inputs
- [aire, Q, reponse, message] = IndependantVectorsEqualVariance(
- vector1 = V1,
- vector2 = V2,
- tolerance = 0.05)
- #
- # Verification par les calculs sans les routines de scipy.stats (ref numerical recipes)
- n1 = V1.size
- n2 = V2.size
- df= n1 + n2 -2
- # Les routines de scipy.stats utilisent une variance calculee avec n-1 et non n comme dans std
- var1 = numpy.sqrt(n1)/numpy.sqrt(n1-1)* V1.std() * numpy.sqrt(n1)/numpy.sqrt(n1-1) * V1.std()
- var2 = numpy.sqrt(n2)/numpy.sqrt(n2-1)* V2.std() * numpy.sqrt(n2)/numpy.sqrt(n2-1) * V2.std()
- m1 = V1.mean()
- m2 = V2.mean()
- var = ((n1 -1.) *var1 + (n2 -1.) *var2 ) /df
- tverif = (m1 -m2) /numpy.sqrt(var*(1./n1 + 1./n2))
- aireverif = 100. * betai(0.5*df,0.5,float(df)/(df+tverif*tverif))
- #
- if (aireverif - aire < 1.e-5) :
- print " Le calcul est conforme à celui de l'algorithme du Numerical Recipes"
- else :
- raise ValueError("Le calcul n'est pas conforme à celui de l'algorithme Numerical Recipes")
-
- if (numpy.abs(aire - 78.42572)< 1.e-5) :
- print " Le calcul est JUSTE sur cet exemple."
- else :
- raise ValueError("Le calcul est FAUX sur cet exemple.")
-
- print
