R vs Python for machine learning

In my last post, I talked about tuning an svm for the Kaggle competition. I submitted my tuned svm. My score on the leaderboard is .90350. Not only am I 198 on the leaderboard and sinking fast, but I didn't even reach the SVM Benchmark score. Additionally, the top person is at a score of .99031.

I figured that only an ensemble method would get me to a higher score and I started to experiment with these methods. I never managed to come up with an ensemble that even matched my original submission.  While I did this, I noticed some things about sci kit learn in Python that made me start to think about looking for other tools.

I decided to try R and Rapid Miner.  Rapid Miner has not been a successful experience. I can't seem to get passed the set up repository/import data stage. I have had much more success with R. Most of this is due to a wonderful set of videos by David Mease. If you are interested in learning R for data mining and machine learning, his videos are pure gold. There are 13 videos on Youtube. Not only does he show you how to use R, but he has all the example data sets online so that you can play along. He also does a wonderful job of explaining what benchmarks to use.

David uses a subset of a well known sonar data set. He uses 130 observations in the training set and 78 observations in the test set. There are 60 features. He goes over several methods with the same data set. I still have one more video, but so far he has covered decision trees, svm and k nearest neighbors. He uses k nearest neighbors with n=1 as a benchmark. This is the default in R. For this data set, it gives a missclassification rate of 21%. This is better that the decision tree misclassification rate which is about 30%. But the svm should be able to beat the untuned k nearest neighbors.

I used this same sonar data set to compare results in R and Python.

k nearest neighbors
Missclassification rate for R: 21%
Missclassification rate for Python: could not get this. I set the n_neighbors=1, but I got this error:

C:\Python27\lib\site-packages\sklearn\neighbors\classification.py:131: NeighborsWarning: kneighbors: neighbor k+1 and neighbor k have the same distance: results will be dependent on data order.

  neigh_dist, neigh_ind = self.kneighbors(X)

The default distance in k nearest neighbors is the Euclidian distance. The data should be scaled so that the variances of each variable are equal. R does this automatically. Python requires you to scale the data yourself.

Decision Tree

The following table shows the results I got:

Depth	R training accuracy	R test accuracy	Python training accuracy	Python test accuracy
1	.7769	.7179	.7769	.7179
2	.80	.7051	.8077	.7051
3	.8615	.6538	.8923	.6667
4	.8846	.6923	.9385	.7179
5	.8846	.6923	.9846	.7436
6	N/A		1.0	.7308

Note that the results are the same for a max depth of 1 and 2. As the max depth increases, it looks like sci kit learn gives the better results. However, the test accuracy stays fairly flat for both models while the Python model training accuracy increase to 1.0. It certainly looks like max depth 4 and 5 in Python have overfit the data. It would be nice to compare a picture of the two trees. The tree in R is quite easy to generate. Python requires some graphics modules that are fairly involved to use. At least, they were for me. I couldn't get either one to work. The R model won't fit max depth 6 because of overfitting issues.

Support Vector Machines

The first thing I did is run a default support vector machine in R and Python. Both programs use an rbf kernel as default.

R scales the data and uses cost=1 and gamma=1/number of features as default values. The untuned svm gives a missclassification error of 1.5% for the training data and  about 13% for the test data.

Python doesn't scale the data and neither did I. (Maybe this is not a fair comparison but it is an extra step in sci kit learn that isn't required in R.) The Python default values are C=1 (cost=1) and gamma=0. This untuned svm gives a missclassification error of about 30% for the training data and about 36% for the test data. 

I've already talked in a previous post about how the Python grid search crashes my computer. R has a procedure for tuning the svm, but it produces an error when I try to run it.

In addition to the questions I have about how sci kit learn models fit the data, there is the additional problem of categorical data. R usually recognizes categorical data. If it doesn't, you can set a variable to be categorical and R will know how to handle it. Python requires you to transform your own categorical data and it is a klugy process. There is a module called OneHotEncoder. But you can't run this module unless you transform all of your text data to numeric.

I still have a lot to learn about machine learning in R. But from I've seen so far, I think I'll stick to R when I want to run a machine learning algorithm.

Tuning C and gamma in the SVM model

In the previous post, the best SVM model for the Kaggle data is an rbf kernel. I'd like to find the best parameters for C and gamma. I did do a grid search with dictionary values for the parameters. I would have liked to have scikit learn search for the best parameters for me and I tried this code:

C_range = 10.0 ** np.arange(-2, 9)

gamma_range = 10.0 ** np.arange(-5, 4)

param_grid = dict(gamma=gamma_range, C=C_range)

cv = StratifiedKFold(cl, n_folds=3)

grid = GridSearchCV(SVC(), param_grid=param_grid, cv=cv)

This code should search for the best values for C in the range .01 to 1,000,000,000 and gamma in the range .00001 to 1000. Unfortunately, this just makes my computer crash and I never got an answer. 

But I found a reference here that says an exhaustive grid search is time consuming. It suggests to use a coarse search first and then a fine search when you are in the correct region.

So starting from C=10 and gamma = .01, I refined my search and here are the values that I got:

C           gamma           score

10            .01                .898

 9           .0095              .898

 8            .009               .901

 7            .0085             .901

 6            .0085             .901

5.5          .0085             .901

5.4          .0085             .901

5.3          .0086             .901

5.28        .0086             .901

You can see that started by using dictionary values for C with a step of 1 on each side of 10 and a step of .005 on each side of gamma. As the grid search stabilized, I narrowed the step on C to .1 and gamma to .0005. Gamma was very stable and I didn't change the step size. I arbitrarily stopped changing C when the step size reached 0.01.

When I ran these parameters using my 70/30 split on the data, I got a score of .9166. This is a 3.6% improvement on the previous score of .9133.

Reference: A Practical Guide to Support Vector Classification retrieved from http://www.csie.ntu.edu.tw/~cjlin/papers/guide/guide.pdf

More on the Data Science competition in Kaggle

In my last post, I talked about the Data Science competition in Kaggle. In that post, I ran an optimized SVM model with a gaussian kernel. In this post, I'll go a little further into depth about the data and models.

I characterized the data as "well structured". I have already mentioned that the data is continuous with no missing values. I used a combination of numpy and pandas to look for missing values, check the mean and standard deviations of each feature, produce histograms to look for data skewing and outliers and a correlation matrix to see if there were any features that had strong linear correlations. These are not specific statistical tests. But this process gave me a good feel for the data and whether I needed any preprocessing.

Once I determined that I had a good data set, I proceeded to modeling. Since there are no categorical features, I decided not to run any kind of decision tree analysis. Since the response value is a classifier, I started with logistic regression and a linear SVM. Each of these gave a score of .797.

At this point, I decided to try a grid search. Here's the description from the user's guide: GridSearchCV implements a “fit” method and a “predict” method like any classifier except that the parameters of the classifier used to predict is optimized by cross-validation.

Here's the code:

param_grid={'C':[.01,.1,1.0,10.0,100.0],'gamma':[.1,.01,.001,.0001],'kernel':['linear','rbf']}

svr=svm.SVC()

grid=grid_search.GridSearchCV(svr,param_grid)

grid.fit(x_train,y_train)

print "The best classifier is:", grid.best_estimator_

print "The best score is ", grid.best_score_

print "The best parameters are ", grid.best_params_

And here's the results:

The best classifier is: SVC(C=10.0, cache_size=200, class_weight=None, coef0=0.0, degree=3,
  gamma=0.01, kernel=rbf, max_iter=-1, probability=False, shrinking=True,
  tol=0.001, verbose=False)
The best score is  0.898426323319
The best parameters are  {'kernel': 'rbf', 'C': 10.0, 'gamma': 0.01}

Ironically, I had already come up with this optimized model just plugging in values. This is not a quick process. I can't give you the exact amount of time that this takes because I just go off and do something else while it is running. Note that the score is not quite as high as my model in the last post. I'm guessing this is because I split the data and used 70% for training and 30% for testing. I believe this model uses a cross validation which means it only used 70% of the data for cross validation.

I also ran a nearest neighbor model. Here's the code and the results:

from sklearn.neighbors import KNeighborsClassifier

neigh=KNeighborsClassifier()

neigh.fit(x_train,y_train)

y_pred3=neigh.predict(x_test)

neigh_score=neigh.score(x_test,y_test)

print "The score from K neighbors is", neigh_score

cm3=confusion_matrix(y_test,y_pred3)

print "This is the confusion matrix with for K neighbors",(cm3)

The score from K neighbors is 0.883333333333
This is the confusion matrix with for K neighbors [[133  22]
 [ 13 132]]

The score for the K neighbors classifier is almost as high as the optimized SVM with the rbf kernel.

I'd be very interested to hear what others are finding as they analyze this set.

Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.

Kaggle Data Science competition

Kaggle.com is sponsoring another learning competition for machine learning. This one specifically mentions using scikit-sklearn in Python. See the competition details here.

It is amazing how much more is available in scikits just since I have been writing this blog. Recently, I have switched to using Python(x,y) which is a distribution which includes everything you need for machine learning. And it's specifically for Windows!! See the information on this distribution here. You do have to be careful about the plug in though. Specifically, the latest version of scikit-sklearn is .13.1. The version that downloads with Python(x,y) is .12. You'll have to update it. Don't ask me how. I took lots of wrong turns, finally figured it out but probably can't reproduce it.

The data set from Kaggle is well structured. There are 40 features and 999 training examples. The feature data is all continuous and there are no missing values. I was able to write code that gives me the SVM standard score on the leaderboard: .913.

Someday I'll have time to figure out how to use github and I'll post my code there. For now, here's what I have:

import csv as csv

import numpy as np

import pandas as pd

import scipy as sp

import matplotlib.pyplot as plt

# Reading in training data for Kaggle sci kit competition

csv_file_object=csv.reader(open('C:/Users/numbersmom/Dropbox/kaggle sci kit competition/train.csv'))

header=csv_file_object.next()

records=[]

for row in csv_file_object:records.append(row)

records=np.array(records)

records=records.astype(np.float)

csv_file_object=csv.reader(open('C:/Users/numbersmom/Dropbox/kaggle sci kit competition/train_label.csv'))

header=csv_file_object.next()

cl=[]

for row in csv_file_object:cl.append(row)

cl=np.array(cl)

cl=cl.astype(np.int8)

cl=cl.reshape(999,)

tr_ex=np.size(cl)

#Need to use 70% of the data for training and 30% for testing

n_train=int(.7*tr_ex)

x_train,x_test=records[:n_train,:],records[n_train:,:]

y_train,y_test=cl[:n_train],cl[n_train:]

#SVM code

from sklearn import svm

from sklearn.svm import SVC

from sklearn.metrics import confusion_matrix

# I tried different models, but this one with c=10 and gamma=.01 gives

# gives the SVM benchmark score.

clf=svm.SVC(C=10.0,gamma=.01,kernel='rbf',probability=True)

clf.fit(x_train,y_train)

print clf.n_support_

y_pred1=clf.predict(x_test)

gau_score=clf.score(x_test,y_test)

print"This is the score for rbf model",gau_score

cm1=confusion_matrix(y_test,y_pred1)

print "This is the confusion matrix for rbf model",(cm1)

print "finished"

The confusion matrix looks like this: 

          pred 0         pred 1

act0    141             14

act1     12             133

There's lots of other stuff I can try to get that number higher. You can check out the helpful users guide to get more information.

SVM with Sage

It's been a long time since my last post, but I was very busy.

First, I had to install Virtual Box and figure out how that worked. Then I had to install Sage in the Virtual Box and figure out how that worked. Then I had to figure out how to run an SVM in Sage.

But I've done all that and I want to post the results and the code. Later on, I'll do a post on Virtual Box.

I used the data from Problem set #1 since it is a small set and easy to use. In order to use it in the SVM code, I had to do two things: first I had to combine the y data with the x data into one CSV file. Then I had to rearrange the data. The original data was sorted so that all the y=0 data comes first followed by the y=1 data. Since I wanted to use only 70% of the data for the analysis and 30% of the data to test accuracy, I had to resort the data so that the y=0 and y=1 data were interspersed evenly.

Below are the results and the graph. The first part of the solution is the output from the convex optimization done by the computer. The verbiage under "Optimal solution found" is my output from the analysis. I have three outputs because I used three different penalty constants: 0.1, 10, and 100.

     pcost       dcost       gap    pres   dres
 0: -7.0774e+00 -1.3421e+01  5e+02  2e+01  1e-14
 1: -1.2084e+00 -1.2358e+01  3e+01  8e-01  1e-14
 2: -5.3246e-01 -4.0847e+00  4e+00  3e-02  2e-15
 3: -6.3668e-01 -1.2676e+00  7e-01  5e-03  9e-16
 4: -7.7552e-01 -9.4394e-01  2e-01  1e-03  9e-16
 5: -8.1966e-01 -8.4799e-01  3e-02  6e-05  8e-16
 6: -8.2948e-01 -8.3229e-01  3e-03  5e-06  1e-15
 7: -8.3049e-01 -8.3053e-01  4e-05  6e-08  1e-15
 8: -8.3051e-01 -8.3051e-01  4e-07  6e-10  1e-15
Optimal solution found.
The percent of support vectors for C=.1 (note: this is underfitting) is
17 percent and there are  12 support vectors.
The weight matrix  is  [ 0.59580336  0.59525422]
The model accuracy is  0.896551724138

 

     pcost       dcost       gap    pres   dres
 0:  1.0143e+02 -3.0514e+05  7e+05  5e-01  7e-13
 1:  1.3516e+03 -5.6065e+04  9e+04  4e-02  6e-13
 2:  9.5288e+02 -1.1776e+04  2e+04  7e-03  4e-13
 3:  1.2074e+01 -3.0995e+03  3e+03  3e-15  3e-13
 4: -2.3429e+02 -4.7539e+02  2e+02  2e-15  2e-13
 5: -2.5664e+02 -4.1207e+02  2e+02  2e-15  3e-13
 6: -3.2509e+02 -4.1261e+02  9e+01  8e-15  4e-13
 7: -3.4242e+02 -3.4631e+02  4e+00  5e-15  3e-13
 8: -3.4389e+02 -3.4393e+02  4e-02  3e-15  4e-13
 9: -3.4390e+02 -3.4390e+02  4e-04  3e-15  5e-13
10: -3.4390e+02 -3.4390e+02  4e-06  4e-15  5e-13
Optimal solution found.
The percent of support vectors for C=100 (note: this is overfitting) is
5 percent and there are  4 support vectors.
The weight matrix  is  [ 2.02253604  1.20665256]
The model accuracy is  0.862068965517

 

     pcost       dcost       gap    pres   dres
 0: -6.2666e+01 -3.9389e+03  1e+04  7e-01  9e-14
 1: -2.0139e+01 -8.5498e+02  1e+03  5e-02  9e-14
 2: -1.2994e+01 -1.5479e+02  2e+02  6e-03  3e-14
 3: -2.2810e+01 -6.2679e+01  4e+01  1e-15  3e-14
 4: -2.8171e+01 -4.6504e+01  2e+01  8e-16  2e-14
 5: -3.3784e+01 -4.1343e+01  8e+00  1e-15  3e-14
 6: -3.2401e+01 -3.9318e+01  7e+00  5e-16  2e-14
 7: -3.4709e+01 -3.6954e+01  2e+00  1e-15  4e-14
 8: -3.5487e+01 -3.6151e+01  7e-01  6e-16  4e-14
 9: -3.5741e+01 -3.5894e+01  2e-01  1e-15  4e-14
10: -3.5787e+01 -3.5789e+01  2e-03  1e-15  4e-14
11: -3.5788e+01 -3.5788e+01  2e-05  2e-16  4e-14
Optimal solution found.
The percent of support vectors for C=10 ) is 7 percent and there are  5
support vectors.
The weight matrix  is  [ 0.95881111  0.84382082]
The model accuracy is  0.896551724138

 This data is so robust that it really almost doesn't matter which penalty constant that you use. 

Here is the code:

import csv as csv
import numpy as np
csv_file_object=csv.reader(open(DATA+'reps1.csv'))
header=csv_file_object.next()
records=[]
for row in csv_file_object:records.append(row)
records=np.array(records)
data=records.astype(np.float)
data[:,0][data[:,0]==0]=-1
#Linear model
m=np.size(data[:,0])
testm=int(.7*m)
datatr=data[0:testm,:]
datatest=data[testm+1:m,:]
from numpy import linalg
import cvxopt
import cvxopt.solvers
trx=datatr[:,1:]
trw=datatr[:,0]
y=trw.reshape(testm,1)
n=np.size(trx,1)
K=np.zeros((testm,testm))
for i in range(testm):
    for j in range(testm):
        K[i,j]=np.dot(trx[i],trx[j])
P=cvxopt.matrix(np.outer(trw,trw)*K)
q=cvxopt.matrix(np.ones(testm)*-1)
b=cvxopt.matrix(0.0)
A=cvxopt.matrix(np.array(y),(1,testm))
G1=cvxopt.matrix(np.diag(np.ones(testm)*-1))
G2=cvxopt.matrix(np.eye(testm))
G=cvxopt.matrix(np.vstack((G1,G2)))
h1=cvxopt.matrix(np.zeros(testm))
h2=cvxopt.matrix(np.ones(testm)*10)
h=cvxopt.matrix(np.vstack((h1,h2)))
solution=cvxopt.solvers.qp(P,q,G,h,A,b)
a=np.ravel(solution['x'])
sv=a[a>.00001]
p=np.size(sv)*100/testm
print 'The percent of support vectors for C=10 ) is',p,"percent and there are ", np.size(sv), "support vectors."
v=np.append(datatr,a.reshape(testm,1),axis=1)
n=np.size(v,axis=1)

z=v[v[:,n-1]>.00001]
ls=np.size(z[:,0])
w=np.zeros(np.size(trx, axis=1))
for t in range(ls):
    w=w+z[t,n-1]*z[t,0]*z[t,1:n-1]
print "The weight matrix  is ",w
yp=0
bp2=0
for i in range(ls):
    for j in range(ls):
        bp2=bp2+z[j,n-1]*z[j,0]*np.dot(z[i,1:n-1],z[j,1:n-1])
    yp=yp+z[i,0]
bf=(yp-bp2)/ls
ac=0
bc=0
cc=0
dc=0
l=np.size(datatest[:,0])
y_p=np.zeros(l)
for k in range(l):
    y_p[k]=np.dot(w,datatest[k,1:].reshape(np.size(trx,1),1))+bf
    if y_p[k]<0:
        s=-1
        if datatest[k,0]==s:
            ac=ac+1
        else:
            bc=bc+1
    else:
        s=1
        if datatest[k,0]==s:
            dc=dc+1
        else:
            dc=dc+1
acc=float(ac+dc)/float(l)
print "The model accuracy is ",acc
u1=v[v[:,0]==-1]
u2=v[v[:,0]==1]
from pylab import *
xvar=np.linspace(trx.min(),trx.max(),10)
clf()
x20=(w[0]*xvar)/(-1*w[1])+bf/(w[1]*-1)
x21=(w[0]*xvar+bf-1)/w[1]*-1
x2n1=(w[0]*xvar+bf+1)/w[1]*-1
plot(xvar,x20)
plot(xvar,x21)
plot(xvar,x2n1)
plot(u1[:,1],u1[:,2],'rx')
plot(u2[:,1],u2[:,2],'bx')
plot(z[:,1],z[:,2],'go')
savefig('sageplt.png')

It's not all good with Sage. I had some major problems with the CVXOPT module and some other annoying things that sent me back to IPython. I'll detail those in a later post.

Stuck

Just when you think you have it all figured out.....everything comes undone.

I had just finished all the lectures on Support Vector Machines and I thought I had a good handle on the concept (30000 ft view) and the details (street view). I opened up Problem set 2 only to find that all of the data files are in Matlab format. Not only that, but you have to import a library to solve for the Lagrange multipliers. What a gyp!! You mean we aren't even going to write the code? (I now know how naive that was.) I closed out the problem set with a vague thought of googling converting Matlab data files to Python data files.

Meanwhile, I continued to work on writing a program for Support Vector Machines. The concept of Support Vector Machines (what kind of name is that?) is really interesting. I'm going to attempt a general explanation with no mathematics.

When I did the Newton's method problem, I was looking for a line that cuts the data into two parts. If the data from PS1 represents survived and not survived, then the idea is to have all survived on one side of the line and not survived on the other side. And that's what happens with the data from PS1with about 89% of the data. In a real data set, if you had any outliers, you could add a fudge factor to the model that would give less weight to the outliers. I'm not kidding. That's what the statisticians do. They don't call it a fudge factor. They call it regularization. So it you have outliers and they are far from the line, you can weight them less than the points closer to the line.

Newton's method works well for a linear model with a small amount of features and data points. But when things get large and possible nonlinear, you need something else. 

Support Vector Machines use really complicated mathematics. The idea is to find a boundary that separates the data. Now, the boundary doesn't have to be a line. It can be an oval like in topography maps. Or it can be something more complicated. But if the data behaves well with your model, then the boundary should separate your data. For example, if you boundary is an oval, all survived data should be inside the oval and all not survived should be outside the oval. This boundary is defined by the data points closest to the boundary. For example, if you have a data point at the exact center of the oval, it is safely in the survived region.  I don't need to know anything else about that point. But the points at the outside edge of the inside area become very important. These points define where exactly the border is drawn. And since only these points define that boundary, they are the only points that need to be used in the model. They are call the support vector machines.

Here are some pictures of plots of support vector machines.

Of course, to solve a math problem, we first have to write an equation, then solve. It turns out that the equation that defines this problem is too hard to solve as is. Without getting into details and a lot of hand waving, the problem can be rewritten as something that can be solved with calculus and Lagrange multipliers. Then substituted back into the original equation to get the boundary. I can tell that your eyes have already glazed over. But let me tell you why this is important. It turns out that using Calculus and Lagrange multipliers allows us to turn the original problem into a convex problem that can be solved using convex optimization software. (Think about the parabola that you learned about in algebra: remember that we could find the maximum or minimum of this by using some formulas. This is just a more complex version of that.) The problem is that the only convex optimization software that I found for Python is CVXOPT and it doesn't work with, you guessed it, Windows. At least I can't get it to work. Here are the installation instructions. 

See this instruction:

tar -xvf blas.tgz

This is Linux and the command does not work in Vista. I can unzip this file, but I don't trust any of the instructions.

It turns out that there is a version of Python that already has this installed. It's called Sage and, you guessed it, it doesn't work on Windows. You have to set up a Virtual Box. I guess that is the next order of business.

Getting Ipython to work in the Windows environment

For those of your following along, you know that I have been using Python instead of Matlab or Octave for the Machine learning course.

It turns out that Python is a base for interactive computing. It's like buying the base model of a car: it comes with the standard equipment and not much else. It will get you from Point A to Point B. But it you want to do something fancy, you have to add on.

For interactive data analysis, there a modules which you can import into Python that make scientific computing easier. You can add in Numpy (numerical python) which gives you access to arrays. You can add in Pandas which gives you access to dataframes. Dataframes allow you to treat data as if it is in a spreadsheet. This makes is much easier to summarize the data. I'll do a separate post on dataframes later.

With each new module that you add in, there are new data structures and commands to learn. This makes it incredibly frustrating for a newbie like me.

So when a friend loaned me Wes McKinney's Python for Data Analysis book, I was thrilled. I figured I could just follow along and learn everything I need to know. Of course, life is never that easy as I found out when I got to Chapter 3. In Chapter 3, Mr. McKinney starts using IPython. In order to keep  using the book, I had to install this on my computer which uses Windows Vista. It turns out this is a big problem because all of the instructions for downloading IPython on your computer are written assuming you are using a Linux based system.

I have finally gotten IPython working on my computer, but it took a lot of research and finagling to do it. In order to help you, I'll try to walk you through the steps.

The completely unhelpful documentation for installation can be found at ipython.org.  Click on the link and read the documentation. The only thing I understood when I read that is that I need Python version 2.6 or higher already installed on my computer. I had already installed Python 2.7 so that I could use Numpy, Matplotlib and Pandas. But what are easy_install and pip? The documentation doesn't explain and there is no further information when you click on pypi.

I did find a blog (this is usually the best source for a newbie) that explains it all. Click on this link to get the instructions. Now that you have done all that, you are ready to use IPython and the interactive notebook. You can see a picture of it here.

Here's how I start up the notebook. It's not perfect, but it gets me where I want to go.

Click on the Windows icon circle.
Type cmd in the search box.
The window with the command prompt will open.
You must change the directory. Type cd c:\Python27\scripts
When the command line prompt comes back, type ipython notebook --pylab=inline
This opens up the notebook and allows you to get plots in the notebook and not a separate window. The only problem that I have is that it opens up the notebook in Explorer and it really doesn't work. I just copy the IP address into Firefox and it works for me.

I have just finished all the lectures for Support Vector Machines, so I will be working on the next problem set.