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.

Using matplotlib to plot the answer to ps1

Now that you know how to solve the first problem in Problem set 1, you have to graph the answer. I used matplotlib. The website for this plotting module can be found here.

Here's a copy of code to plot the data for ps 1:

#Plots x1 vs x2 for Problem 1 of the Stanford machine learning class

from pylab import *
import numpy as np
import numpy.linalg
from math import *

x1,x2=np.loadtxt('q1x.dat',unpack=True)
y=np.loadtxt('q1y.dat',unpack=True)

# we need to show when h=.5 as a line to separate the data. h=.5 when
# theta transpose x =0. so 0=theta0 +theta1*x1+theta2*x2
# solving for x2 give -theta0/theta2-theta1/theta2*x1

a=np.array([min(x1),max(x1),1,0.01])
b=-(-2.6205)/1.1719-0.7604/1.1719*a
plot(a,b)
# Use different colors and markers to plot x1 vs x2 if y=0 or y=1

for i in range(0,98):
    if y[i]==0:
        plot(x1[i],x2[i],'ro')
    if y[i]==1:
        plot(x1[i],x2[i],'bx')

xlabel('x1')
ylabel('x2')

show()

I used a for x1 and b as the calculated value for x2. Plotting a vs b gives that beautiful straight line across the plot.

One of my previous posts talked about expecting the unexpected in numpy. It happens again in this little short program. This line:
a=np.array([min(x1),max(x1),1,0.01])
was not what I though it was. I thought this gave me a list of values between the minimum of x1 and the maximum of x1. I was wrong. Here is what this actually gives:
  print a
[ 0.57079941  7.7054006   1.          0.01      ]

If you want a list of values, you need to use this code:

 a=np.linspace(min(z),max(z),25)

This gives 25 values between the minimum and maximum of your data set.

The rest of the code is very straight forward. Just set up a loop. If y=1, the data plots as a red circle. If y=0, the data plots as a blue x.

I know I have posted this graph before but it is so pretty, I'm going to post it again:

Code for the Newton's method problem

As promised, here is my code for Newton's method. Note that I am using Numpy but not pandas. A friend has lent me this book so I am currently working on learning what pandas can do for me.

The first thing you have to do is import numpy. I import as np so all of my numpy commands begin with np. For some reason, the linear algebra module does not load unless I specifically call it with an import command. I need this to invert the Hessian matrix.
import numpy as np
import numpy.linalg
from math import *


Here I read in the data. The size of x along axis=0 (rows) gives the number of training examples. In order to get the intercept term, I have to put a column of ones in front of the x data. In order to get the column of ones where I wanted it, I had to use the insert command. Then I checked the number of columns.
x=np.loadtxt('q1x.dat')
m=np.size(x,0) #number of training examples
x=np.insert(x,[0],[1],axis=1)
n=np.size(x,1) #number of columns


This piece of code sets up a vector of zeros for theta to start the iterations. I do the same thing in the code below for the Hessian matrix and gradient vector. Recall that everything in Python starts at zero so I have an n x 1 vector for the gradient instead of n+1 x 1 like in Matlab.
theta=np.zeros(n)

y=np.loadtxt('q1y.dat')

 Note the reshape command for theta when I am calculating h. For some reason, I couldn't get theta to transpose.
# set up for 10 iterations
for i in range(0,10):
    H=np.zeros((n,n))
    grad=np.zeros(n)
    print "Theta just before h is ", theta
    h=1/(1+e**(np.dot(x,theta.reshape(n,1))))
    ts=np.dot(x,theta.reshape(n,1))

This part of the code is pretty self explanatory and looks very much like the matlab code from the answer sheet. However, note the outer command in the calculation for H. Let's examine this a little more closely.

Here's an example. In this example, a is an example of the lprime vector with one variable (ignore the intercept) and 4 training examples. b is the same vector. Notice they are both 4x1, so you cannot use matrix multiplication.
 a=np.array([[1],[2],[3],[4]])
>>> print a
[[1]
 [2]
 [3]
 [4]]
>>> b=np.array([[1],[2],[3],[4]])
>>> print b
[[1]
 [2]
 [3]
 [4]]
>>> c=np.dot(a,b)

Traceback (most recent call last):
  File "<pyshell#23>", line 1, in <module>
    c=np.dot(a,b)
ValueError: objects are not aligned




However, look at what happens when you use the outer command:
>>> d=np.outer(a,b)
>>> print d
[[ 1  2  3  4]
 [ 2  4  6  8]
 [ 3  6  9 12]
 [ 4  8 12 16]]

You get a 4x4 matrix. Where did this come from? The outer product gives you

out[i, j] = a[i] * b[j]

This is what you need for Hessian matrix. It took me a really, really long time to find this. You're welcome.
So here is the rest of the code:
# This calculates the values for lprime vector which will be summed below
    for j in range(0,m):
     
        grad=grad+(x[j]*(y[j]-h[j]))

        H=H+h[j]*(1-h[j])*np.outer(x[j],x[j])
   
    Hinv=np.linalg.inv(H)
    theta=theta-np.dot(Hinv,grad)


The next post will show the code for the plot.