Voici un exemple d'estimation de régression logistique utilisant l ' algorithme d'optimisation BFGS [L-BFGS] à mémoire limitée. J'utiliserai la fonction optimx de la bibliothèque optimx de R et de scipy.optimize.fmin_l_bfgs_b en Python.
Python
L'exemple que j'utilise est tiré de Sheather (2009, p. 264). Le code Python suivant montre l'estimation de la régression logistique à l'aide de l'algorithme BFGS:
# load required librariesimport numpy as npimport scipy as spimport scipy.optimizeimport pandas as pdimport os # hyperlink to data locationurlSheatherData = "http : //www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv "# lire les données dans un tableau NumPyarrSheatherData = np.asarray (pd.read_csv (urlSheatherData)) # découper les données en obtenir la variable dépendantevY = arrSheatherData [:, 0] .astype ('float64') # découper les données pour obtenir la matrice des variables prédictives mX = np.asarray (arrSheatherData [:, 2:]). astype ('float64') # ajouter une interception aux variables prédictivesintercept = np.ones (mX.shape [0]). reshape (mX.shape [0], 1) mX = np.concatenate ((intercept, mX), axis = 1) # le nombre de variables et d'observationsiK = mX.shape [1] iN = mX.shape [0] # transformation logistiquedef logit (mX, vBeta): return ((np.exp (np.dot (mX, vBeta)) / (1.0 + np .exp (np.dot (mX, vBeta))))) # paramétrage stable de la fonction de coûtdef logLikelihoodLogitStable (vBeta, mX, vY): return (- (np.sum (vY * (np.dot (mX, vBeta) - np.log ((1.0 + np.exp (np.dot (mX, vBeta) ))))) + (1-vY) * (- np.log ((1.0 + np.exp (np.dot (mX, vBeta)))))))) # score functiondef likelihoodScore (vBeta, mX, vY ): retour (np.dot (mX.T, (logit (mX, vBeta) - vY))) # ========================= ============================================ # Optimiser pour obtenir le MLE en utilisant l'optimiseur BFGS (dérivés numériques) # ========================================== ===========================
optimLogitBFGS = sp.optimize.minimize (logLikelihoodLogitStable, x0 = np.array ([10, 0.5, 0.1, -0.3, 0.1]), args = (mX, vY), method = 'BFGS', options = {'gtol' : 1e-3, 'disp': True}) print (optimLogitBFGS) # imprimer les résultats de l'optimisation
Et cela peut facilement être adapté au scipy.optimize.fmin_l_bfgs_b fonction:
# ================================= ==================================== # optimiser pour obtenir le MLE en utilisant l'optimiseur L-BFGS (dérivés analytiques ) # ================================================= ==================== optimLogitLBFGS = sp.optimize.fmin_l_bfgs_b (logLikelihoodLogitStable, x0 = np.array ([10, 0.5, 0.1, -0.3, 0.1]), args = (mX, vY), fprime = likelihoodScore, pgtol = 1e-3, disp = True) print (optimLogitLBFGS) # affiche le résultat lts de l'optimisation
R
L'utilisation de l'optimiseur L-BFGS-B dans R est tout aussi simple. D'abord la version avec l'algorithme BFGS:
library (optimx) # read in the dataurlSheatherData = "http://www.stat.tamu.edu/~sheather/book/docs/datasets/ MichelinNY.csv "dfSheatherData = as.data.frame (read.csv (urlSheatherData, header = T)) # créer les matrices de conceptionvY = as.matrix (dfSheatherData ['InMichelin']) mX = as.matrix (dfSheatherData [c ( 'Service', 'Decor', 'Food', 'Price')]) # ajouter une interception aux variables prédictives mX = cbind (rep (1, nrow (mX)), mX) # le nombre de variables et d'observations iK = ncol (mX) iN = nrow (mX) # définir la transformation logistiquelogit = function (mX, vBeta) {return (exp (mX% *% vBeta) / (1+ exp (mX% *% vBeta)))} # paramétrage stable de la fonction log-vraisemblance # Remarque: Le négatif du log-vraisemblance est renvoyé, puisque nous allons # / minimiser / la fonction.logLikelihoodLogitStable = function (vBeta, mX, vY) {return (-sum (
vY * (mX% *% vBeta - log (1 + exp (mX% *% vBeta))) + (1-vY) * (- log (1 + exp (mX% *% vBeta)))) # somme) # return} # score functionlikelihoodScore = function (vBeta, mX, vY) {return (t (mX)% *% (logit (mX, vBeta) - vY))} # jeu initial de paramètres vBeta0 = c (10, -0.1, -0,3, 0,001, 0,01) # paramètres de départ arbitraires # ====================================== =============================== # Optimiser pour obtenir le MLE en utilisant l'optimiseur BFGS (dérivés numériques) # ===== =================================================== ============= optimLogitBFGS = optim (vBeta0, logLikelihoodLogitStable, mX = mX, vY = vY, method = 'BFGS', hessian = TRUE) optimLogitBFGS # obtenir les résultats de l'optimisation
puis la version avec le L-BFGS-B du package optimx :
# ======== =================================================== ========== # Optimiser pour obtenir le MLE en utilisant l'optimiseur L-BFGS (dérivés analytiques) # ======================= ======================= ======================= optimLogitLBFGS = optimx (vBeta0, logLikelihoodLogitStable, method = 'L-BFGS-B', gr = likelihoodScore, mX = mX, vY = vY, hessian = TRUE) résumé (optimLogitLBFGS)