I am used to doing all my statistics in R and python for all the peripheral tasks. Just for fun I attempted a BFGS optimization to compare it to the ordinary LS result - both in python using scipy/numpy. But the results do not match. I do not see any errors. I also attach the equivalent code in R (which works). Can anyone correct my usage of scipy.optimize.fmin_bfgs to match the OLS or R results ?
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
The result of the optimization is:
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]
But it really should match the OLS results achieved in lsBetas = np.linalg.lstsq(Xs, Ys):
[[-1.49209249]
[ 0.05773374]
[ 0.044193 ]
[-0.01117662]
[ 0.00179794]]
Here is the R code in case it is useful (it also has the advantage of being able to read directly from the 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)
First, let's make arrays out of list:
>>> Xs = np.vstack(Xs)
>>> Ys = np.vStack(Ys)
Then, fnRSS is incorrectly translated, it's argument, beta, is passed transposed. Can be fixed with
>>> def fnRSS(vBeta, vY, vX):
... return np.sum((vY.T - np.dot(vX, vBeta))**2)
Final result:
>>> 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, consider using pandas read_csv parser or numpy genfromtxt or recfromcsv to read csv data into array instead of custom written parsers. No problems to read from url either:
>>> 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