src/daComposant/daNumerics/mmqr.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 # Copyright (C) 2008-2013 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 __doc__ = """
  24     Implémentation informatique de l'algorithme MMQR, basée sur la publication :
  25     David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
  26     Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
  27 """
  28 __author__ = "Jean-Philippe ARGAUD"
  29
  30 import sys, math
  31 from numpy import array, matrix, asarray, asmatrix
  32 from numpy import sum, dot, linalg, ravel, max, min, hstack, argmin, argmax
  33
  34 # ==============================================================================
  35 def mmqr(
  36         func     = None,
  37         x0       = None,
  38         fprime   = None,
  39         bounds   = None,
  40         quantile = 0.5,
  41         maxfun   = 15000,
  42         toler    = 1.e-06,
  43         y        = None,
  44         ):
  45     #
  46     # Recuperation des donnees et informations initiales
  47     # --------------------------------------------------
  48     variables = asmatrix(ravel( x0 ))
  49     mesures   = asmatrix(ravel( y )).T
  50     increment = sys.float_info[0]
  51     p         = len(variables.flat)
  52     n         = len(mesures.flat)
  53     quantile  = float(quantile)
  54     #
  55     # Calcul des parametres du MM
  56     # ---------------------------
  57     tn      = float(toler) / n
  58     e0      = -tn / math.log(tn)
  59     epsilon = (e0-tn)/(1+math.log(e0))
  60     #
  61     # Calculs d'initialisation
  62     # ------------------------
  63     residus  = ravel( mesures - func( variables ) )
  64     poids    = asarray( 1./(epsilon+abs(residus)) )
  65     veps     = 1. - 2. * quantile - residus * poids
  66     lastsurrogate = -sum(residus*veps) - (1.-2.*quantile)*sum(residus)
  67     iteration = 0
  68     #
  69     # Recherche iterative
  70     # -------------------
  71     while (increment > toler) and (iteration < maxfun) :
  72         iteration += 1
  73         #
  74         Derivees  = array(fprime(variables))
  75         Derivees  = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
  76         DeriveesT = array(matrix(Derivees).T)
  77         M         =   dot( DeriveesT , (array(matrix(p*[poids,]).T)*Derivees) )
  78         SM        =   dot( DeriveesT , veps ).T
  79         step      = - linalg.lstsq( M, SM )[0]
  80         #
  81         variables = variables + step
  82         if bounds is not None:
  83             while( (variables < ravel(asmatrix(bounds)[:,0])).any() or (variables > ravel(asmatrix(bounds)[:,1])).any() ):
  84                 step      = step/2.
  85                 variables = variables - step
  86         residus   = ravel( mesures - func(variables) )
  87         surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  88         #
  89         while ( (surrogate > lastsurrogate) and ( max(list(abs(step))) > 1.e-16 ) ) :
  90             step      = step/2.
  91             variables = variables - step
  92             residus   = ravel( mesures-func(variables) )
  93             surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  94         #
  95         increment     = lastsurrogate-surrogate
  96         poids         = 1./(epsilon+abs(residus))
  97         veps          = 1. - 2. * quantile - residus * poids
  98         lastsurrogate = -sum(residus * veps) - (1.-2.*quantile)*sum(residus)
  99     #
 100     # Mesure d'écart : q*Sum(residus)-sum(residus negatifs)
 101     # ----------------
 102     Ecart = quantile * sum(residus) - sum( residus[residus<0] )
 103     #
 104     return variables, Ecart, [n,p,iteration,increment,0]
 105
 106 # ==============================================================================
 107 if __name__ == "__main__":
 108     print '\n AUTODIAGNOSTIC \n'