1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2019 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
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.
28 __author__ = "Jean-Philippe ARGAUD"
31 from numpy import array, matrix, asarray, asmatrix
32 from numpy import sum, dot, linalg, ravel, max, min, hstack, argmin, argmax
34 # ==============================================================================
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)
53 quantile = float(quantile)
55 # Calcul des parametres du MM
56 # ---------------------------
58 e0 = -tn / math.log(tn)
59 epsilon = (e0-tn)/(1+math.log(e0))
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)
71 while (increment > toler) and (iteration < maxfun) :
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, rcond=-1 )[0]
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() ):
85 variables = variables - step
86 residus = ravel( mesures - func(variables) )
87 surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
89 while ( (surrogate > lastsurrogate) and ( max(list(abs(step))) > 1.e-16 ) ) :
91 variables = variables - step
92 residus = ravel( mesures-func(variables) )
93 surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus)
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)
100 # Mesure d'écart : q*Sum(residus)-sum(residus negatifs)
102 Ecart = quantile * sum(residus) - sum( residus[residus<0] )
104 return variables, Ecart, [n,p,iteration,increment,0]
106 # ==============================================================================
107 if __name__ == "__main__":
108 print('\n AUTODIAGNOSTIC \n')