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 sum, array, matrix, dot, linalg, asarray, asmatrix, ravel
  32
  33 # ==============================================================================
  34 def mmqr(
  35         func     = None,
  36         x0       = None,
  37         fprime   = None,
  38         quantile = 0.5,
  39         maxfun   = 15000,
  40         toler    = 1.e-06,
  41         y        = None,
  42         ):
  43     #
  44     # Recuperation des donnees et informations initiales
  45     # --------------------------------------------------
  46     variables = asmatrix(ravel( x0 ))
  47     mesures   = asmatrix(ravel( y )).T
  48     increment = sys.float_info[0]
  49     p         = len(variables.flat)
  50     n         = len(mesures.flat)
  51     quantile  = float(quantile)
  52     #
  53     # Calcul des parametres du MM
  54     # ---------------------------
  55     tn      = float(toler) / n
  56     e0      = -tn / math.log(tn)
  57     epsilon = (e0-tn)/(1+math.log(e0))
  58     #
  59     # Calculs d'initialisation
  60     # ------------------------
  61     residus  = asmatrix( mesures - func( variables ) ).A1
  62     poids    = asarray( 1./(epsilon+abs(residus)) )
  63     veps     = 1. - 2. * quantile - residus * poids
  64     lastsurrogate = -sum(residus*veps) - (1.-2.*quantile)*sum(residus)
  65     iteration = 0
  66     #
  67     # Recherche iterative
  68     # -------------------
  69     while (increment > toler) and (iteration < maxfun) :
  70         iteration += 1
  71         #
  72         Derivees  = array(fprime(variables))
  73         Derivees  = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
  74         DeriveesT = array(matrix(Derivees).T)
  75         M         = - dot( DeriveesT , (array(matrix(p*[poids,]).T)*Derivees) )
  76         SM        =   dot( DeriveesT , veps ).T
  77         step      = linalg.lstsq( M, SM )[0]
  78         #
  79         variables = variables + step
  80         residus   = asmatrix( mesures - func(variables) ).A1
  81         surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  82         #
  83         while ( (surrogate > lastsurrogate) and ( max(list(abs(step))) > 1.e-16 ) ) :
  84             step      = step/2.
  85             variables = variables - step
  86             residus   = ( mesures-func(variables) ).A1
  87             surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  88         #
  89         increment     = lastsurrogate-surrogate
  90         poids         = 1./(epsilon+abs(residus))
  91         veps          = 1. - 2. * quantile - residus * poids
  92         lastsurrogate = -sum(residus * veps) - (1.-2.*quantile)*sum(residus)
  93     #
  94     # Mesure d'écart : q*Sum(residus)-sum(residus negatifs)
  95     # ----------------
  96     Ecart = quantile * sum(residus) - sum( residus[residus<0] )
  97     #
  98     return variables, Ecart, [n,p,iteration,increment,0]
  99
 100 # ==============================================================================
 101 if __name__ == "__main__":
 102     print '\n AUTODIAGNOSTIC \n'