Правильное использование scipy.optimize.fmin_bfgs требуется по сравнению с кодом R

Question

Я использую все свои статистические данные в R и python для всех периферийных задач. Просто для удовольствия я попытался оптимизировать BFGS, чтобы сравнить его с обычным результатом LS - как в python, использующем scipy/numpy. Но результаты не совпадают. Я не вижу никаких ошибок. Я также прикрепляю эквивалентный код в R (который работает). Может ли кто-нибудь исправить мое использование scipy.optimize.fmin_bfgs, чтобы соответствовать результатам OLS или R?

import csv 
import numpy as np 
import scipy as sp 
from scipy import optimize 

class DataLine: 
    def __init__(self,row): 
     self.Y = row[0] 
     self.X = [1.0] + row[2:len(row)] 
     # 'Intercept','Food','Decor', 'Service', 'Price' and remove the name 
    def allDataLine(self): 
     return self.X + list(self.Y) # return operator.add(self.X,list(self.Y)) 
    def xData(self): 
     return np.array(self.X,dtype="float64") 
    def yData(self): 
     return np.array([self.Y],dtype="float64") 
def fnRSS(vBeta, vY, mX): 
    return np.sum((vY - np.dot(mX,vBeta))**2) 
if __name__ == "__main__": 
    urlSheatherData = "/Hans/workspace/optimsGLMs/MichelinNY.csv" 
    # downloaded from "http://www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv" 
    reader = csv.reader(open(urlSheatherData), delimiter=',', quotechar='"') 
    headerTuple = tuple(reader.next()) 
    dataLines = map(DataLine, reader) 
    Ys = map(DataLine.yData,dataLines) 
    Xs = map(DataLine.xData,dataLines) 
    # a check and an initial guess ... 
    vBeta = np.array([-1.5, 0.06, 0.04,-0.01, 0.002]).reshape(5,1) 
    print np.sum((Ys-np.dot(Xs,vBeta))**2) 
    print fnRSS(vBeta,Ys,Xs) 
    lsBetas = np.linalg.lstsq(Xs, Ys) 
    print lsBetas[1] 
    # prints the right numbers 
    print lsBetas[0] 
    optimizedBetas = sp.optimize.fmin_bfgs(fnRSS, x0=vBeta, args=(Ys,Xs)) 
    # completely off .. 
    print optimizedBetas 

Результатом оптимизации является:

Optimization terminated successfully. 
     Current function value: 6660.000006 
     Iterations: 276 
     Function evaluations: 448 
[ 4.51296549e-01 -5.64005114e-06 -3.36618459e-06 4.98821735e-06 
    9.62197362e-08] 

Но это действительно должно соответствовать результаты МНК достигнуты в lsBetas = np.linalg.lstsq(Xs, Ys):

[[-1.49209249] 
[ 0.05773374] 
[ 0.044193 ] 
[-0.01117662] 
[ 0.00179794]] 

Здесь R-код в случае полезно (он также имеет то преимущество, что он может читать непосредственно из URL-адреса):

urlSheatherData = "http://www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv" 
dfSheather = as.data.frame(read.csv(urlSheatherData, header = TRUE)) 
vY = as.matrix(dfSheather['InMichelin']) 
mX = as.matrix(dfSheather[c('Service','Decor', 'Food', 'Price')]) 
mX = cbind(1, mX) 
fnRSS = function(vBeta, vY, mX) { return(sum((vY - mX %*% vBeta)^2)) } 
vBeta0 = rep(0, ncol(mX)) 
optimLinReg = optim(vBeta0, fnRSS,mX = mX, vY = vY, method = 'BFGS', hessian=TRUE) 
print(optimLinReg$par) 

Answer 1

Во-первых, давайте создадим массивы из списка:

>>> Xs = np.vstack(Xs) 
>>> Ys = np.vStack(Ys) 

Затем fnRSS неправильно переведено, это аргумент, бета, передается транспонировать. Может быть установлено с

>>> def fnRSS(vBeta, vY, vX): 
...  return np.sum((vY.T - np.dot(vX, vBeta))**2) 

Окончательным результатом:

>>> sp.optimize.fmin_bfgs(fnRSS, x0=vBeta, args=(Ys,Xs)) 
Optimization terminated successfully. 
     Current function value: 26.323906 
     Iterations: 9 
     Function evaluations: 98 
     Gradient evaluations: 14 
array([-1.49208546, 0.05773327, 0.04419307, -0.01117645, 0.00179791]) 

Sidenote, рассмотреть возможность использования панд read_csv парсера или NumPy genfromtxt или recfromcsv для чтения CSV данных в массив вместо пользовательских написано парсеров. Нет проблем с чтением с URL:

>>> import pandas as pd 
>>> urlSheatherData = "http://www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv" 
>>> data = pd.read_csv(urlSheatherData) 
>>> data[['Service','Decor', 'Food', 'Price']].head() 
    Service Decor Food Price 
0  19  20 19  50 
1  16  17 17  43 
2  21  17 23  35 
3  16  23 19  52 
4  19  12 23  24 

[5 rows x 4 columns] 
>>> data['InMichelin'].head() 
0 0 
1 0 
2 0 
3 1 
4 0 
Name: InMichelin, dtype: int64