src/daComposant/daNumerics/mmqr.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2012 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
  22 # L'algorithme est base sur la publication : David R. Hunter, Kenneth Lange,
  23 # "Quantile Regression via an MM Algorithm", Journal of Computational and
  24 # Graphical Statistics, 9, 1, pp.60-77, 2000
  25
  26 import sys, math
  27 from numpy import sum, array, matrix, dot, linalg, asarray, asmatrix
  28
  29 # ==============================================================================
  30 def mmqr(
  31         func     = None,
  32         x0       = None,
  33         fprime   = None,
  34         quantile = 0.5,
  35         maxfun   = 15000,
  36         toler    = 1.e-06,
  37         y        = None,
  38         ):
  39     #
  40     # Recuperation des donnees et informations initiales
  41     # --------------------------------------------------
  42     variables = asmatrix(x0).A1
  43     mesures   = asmatrix(y).flatten().T
  44     increment = sys.float_info[0]
  45     p         = len(variables.flat)
  46     n         = len(mesures.flat)
  47     quantile  = float(quantile)
  48     #
  49     # Calcul des parametres du MM
  50     # ---------------------------
  51     tn      = float(toler) / n
  52     e0      = -tn / math.log(tn)
  53     epsilon = (e0-tn)/(1+math.log(e0))
  54     #
  55     # Calculs d'initialisation
  56     # ------------------------
  57     residus  = asmatrix( mesures - func( variables ) ).A1
  58     poids    = asarray( 1./(epsilon+abs(residus)) )
  59     veps     = 1. - 2. * quantile - residus * poids
  60     lastsurrogate = -sum(residus*veps) - (1.-2.*quantile)*sum(residus)
  61     iteration = 0
  62     #
  63     # Recherche iterative
  64     # -------------------
  65     while (increment > toler) and (iteration < maxfun) :
  66         iteration += 1
  67         #
  68         Derivees  = array(fprime(variables))
  69         Derivees  = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
  70         DeriveesT = array(matrix(Derivees).T)
  71         M         = - dot( DeriveesT , (array(matrix(p*[poids,]).T)*Derivees) )
  72         SM        =   dot( DeriveesT , veps ).T
  73         step      = linalg.lstsq( M, SM )[0]
  74         #
  75         variables = variables + step
  76         residus   = asmatrix( mesures - func(variables) ).A1
  77         surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  78         #
  79         while ( (surrogate > lastsurrogate) and ( max(list(abs(step))) > 1.e-16 ) ) :
  80             step      = step/2.
  81             variables = variables - step
  82             residus   = ( mesures-func(variables) ).A1
  83             surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  84         #
  85         increment     = lastsurrogate-surrogate
  86         poids         = 1./(epsilon+abs(residus))
  87         veps          = 1. - 2. * quantile - residus * poids
  88         lastsurrogate = -sum(residus * veps) - (1.-2.*quantile)*sum(residus)
  89     #
  90     # Mesure d'écart : q*Sum(residus)-sum(residus negatifs)
  91     # ----------------
  92     Ecart = quantile * sum(residus) - sum( residus[residus<0] )
  93     #
  94     return variables, Ecart, [n,p,iteration,increment,0]
  95
  96 # ==============================================================================
  97 if __name__ == "__main__":
  98     print '\n AUTODIAGNOSTIC \n'