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(asmatrix(y).A1).T
  44     increment = sys.float_info[0]
  45     p         = len(variables.flat)
  46     n         = len(mesures.flat)
  47     #
  48     # Calcul des parametres du MM
  49     # ---------------------------
  50     tn      = float(toler) / n
  51     e0      = -tn / math.log(tn)
  52     epsilon = (e0-tn)/(1+math.log(e0))
  53     #
  54     # Calculs d'initialisation
  55     # ------------------------
  56     residus  = asmatrix( mesures - func( variables ) ).A1
  57     poids    = 1./(epsilon+abs(residus))
  58     veps     = 1. - 2. * quantile - residus * poids
  59     lastsurrogate = -sum(residus*veps) - (1.-2.*quantile)*sum(residus)
  60     iteration = 0
  61     #
  62     # Recherche iterative
  63     # -------------------
  64     while (increment > toler) and (iteration < maxfun) :
  65         iteration += 1
  66         #
  67         Derivees  = fprime(variables)
  68         DeriveesT = matrix(Derivees).T
  69         M         = - dot( DeriveesT , (array(matrix(p*[poids]).T)*array(Derivees)) )
  70         SM        =   dot( DeriveesT , veps ).T
  71         step      = linalg.lstsq( M, SM )[0].A1
  72         #
  73         variables = asarray(variables) + asarray(step)
  74         residus   = ( mesures - func(variables) ).A1
  75         surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  76         #
  77         while ( (surrogate > lastsurrogate) and ( max(list(abs(step))) > 1.e-16 ) ) :
  78             step      = step/2.
  79             variables = variables - step
  80             residus   = ( mesures-func(variables) ).A1
  81             surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
  82         #
  83         increment     = lastsurrogate-surrogate
  84         poids         = 1./(epsilon+abs(residus))
  85         veps          = 1. - 2. * quantile - residus * poids
  86         lastsurrogate = -sum(residus * veps) - (1.-2.*quantile)*sum(residus)
  87     #
  88     # Mesure d'écart : q*Sum(residus)-sum(residus negatifs)
  89     # ----------------
  90     Ecart = quantile * sum(residus) - sum( residus[residus<0] )
  91     #
  92     return variables, Ecart, [n,p,iteration,increment,0]
  93
  94 # ==============================================================================
  95 if __name__ == "__main__":
  96     print '\n AUTODIAGNOSTIC \n'