python - 与 R 代码相比，需要正确使用 scipy.optimize.fmin_bfgs

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) 中实现的 OLS 结果相匹配:

[[-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 被错误翻译，它的参数 beta 被转置。可以通过修复

>>> 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])

旁注，考虑使用 pandas 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