Python as a data analysis platform

Despite the fact that I've been aware of Scikits Learn (sklearn) for some time, I never got the chance to really use Python for data analysis and, instead, I was a victim of my own inertia and limited myself to use R and especially Matlab.

I must say, in the beginning, Python looks awkward: it was inconceivable for me to use an invisible element (spaces or tabs) as a structural construction of a program (defining blocks), in a way much similar to Fortran, which I always considered weird (coming from the C world). This and the lack of the omnipresent, C-syntax end-of-line semicolon, prove to be a major boosting element when programming in Python. I must say that whatever lack in computer performance is overcome by the speed the programmer experiences when writing the software. This applies to general software, such as the App server that I am preparing, which is being written in Python using the Google App Engine, and I have to say that it just runs smoothly, no need for recompilations, clear syntax and one-line complex data-processing pieces of code.

Regarding data analysis, it is a little more complicated than Matlab's clear orientation towards numerical linear algebra (where everything is a Matrix). Good comparisons and reasons supporting my view are

• https://sites.google.com/site/pythonforscientists/python-vs-matlab
• http://www.stat.washington.edu/~hoytak/blog/whypython.html
• http://stevetjoa.com/305/

It was precisely the last blog the one that spurred me to give it a try.

Now, going to Machine Learning specifics, sklearn has everything you need for the majority of the work a machine learning practitioner will ever need.

Data preprocessors, including text vectorizers and TF IDF preprocessors

SVM implementations

Stochastic Gradient Descent algorithms for fast regression and classification

Random Forest and other ensemble methods for robust regression and classification

Clustering algorithms

Data dimensionality reduction algorithms such as LLE, ISOMAP and spectral embeddings

Results presentation, including mean squared error for regression and precision/recall tables for classification. It even computes the area under the ROC curve.

This, added to the clean, standardized and well-designed interface, which always has a .fit method for every object which performs the task of learning from samples, and then either a .transform method if the learning is unsupervised (such as LLE, ISOMAP, ICA, PCA, or the preprocessors, etc) or .predict if the learning is supervised (SVM, SGD, ensemble...). If enables a pipelining mechanism that allows us to build the whole pipeline from data reading to results output.

One of the lead programmers of the project, Andreas Müller has a very insightful blog. Check it out in the following URL

peekaboo-vision.blogspot.com.es

I decided to be more active on Kaggle. For the moment I scored 13th on the Leaderboard of the Amazon employee access competition that recently opened. Competing against Alexander Larko or any of the other high-standing data scientists chills my blood.

Last but not least, just to comment that future work seems to be bent on using the GPU to perform all the linear algebra. Check out
Gnumpy: http://www.cs.toronto.edu/~tijmen/gnumpy.html
Deep Belief Networks: http://deeplearning.net/tutorial/DBN.html
PyCUDA: http://documen.tician.de/pycuda/tutorial.html

Machine Learning video lectures (Beginners)

The first one is part of the PyData workshop held in March this year. This lecture is oriented to beginners, so if you already know the (very) basics, then it won't be of interest to you. It is nice, however, since you can see all the advanced capabilities of Scikits Learn while watching something you already know, such as linear, polynomial and Gaussian SVMs, logistic regression and naive bayes.

Improve the performance of your SVM

Though I am not very keen on differential geometry (others aren't either, but they claim to be researching on the field), I find it amusing to read a little bit of it when it is used along with kernel methods, and especially when you can improve the behavior of a SVM with it.

Amari and Wu are responsible for the following method: The idea is that, in order to increase class separability, we need to enlarge the spatial resolution around the boundary in the feature space. Take, for instance, the Riemannian distance along the manifold
$$
ds^2 = \sum_{i,j} g_{i,j} dx_i dx_j
$$
We need it to be large along the border of $f(\mathbf{x})=0$ and small between points of the same class. In practice, the boundary is not known, so we use the points the we know are closest to the boundary: the support vectors. A conformal transformation does the job
$$
\tilde{g}_{i,j}(\mathbf{x}) = \Omega (\mathbf{x}) g_{i,j} (\mathbf{x})
$$

This is very difficult to realize practically, so we consider a quasi-conformal transformation to induce the a similar map by directly modifying
$$
\tilde{K}(\mathbf{x_1},\mathbf{x_2}) = c(\mathbf{x_1}) c(\mathbf{x_2}) K(\mathbf{x_1},\mathbf{x_2})
$$
where $c(\mathbf{x})$ is a positive function, that can be built from the data as
$$
c(\mathbf{x}) = \sum_{i \in SV} h_i e^{\frac{\| \mathbf{x} - \mathbf{x}\|^2}{2\tau^2}}
$$
where $h_i$ is a parameter of the $i$-th support vector.

Thus, if you first train a SVM with a standard kernel, and then you compute $c(x)$ and make a new kernel with the previous expressions, your SVM will behave better.

The authors report higher classification accuracy and less support vectors than with standard kernels.

Check out the paper:
http://www.dcs.warwick.ac.uk/~feng/papers/Scaling%20the%20Kernel%20Function.pdf

L. Wasserman's & O. Bousquet's blogs

I was reading UCL CSML news and an entry of Larry Wasserman's blog popped up. Yes, this is the Wasserman's responsible for the statistical bible that can be used to exercise your mind by learning statistics and to exercise your muscles by lifting and putting the book on the desktop.

Also, I stumbled upon Olivier Bousquet's blog. It contains some thoughts about Machine Learning that will be of interest to practitioners.

I added a link to Wasserman's and Bousquet's blogs to the right.

Ohh hell: I just found out that the last entry of Bousquet's blog is 5 years old. Shame.

SVM plot decision function

In the paper I submitted, I have to deal with SVMs and I wanted to plot the decision function that my kernel made with a 2D dataset.

I stumbled upon a 2D problem that I am interested in visualizing (I already posted about it). The point is that I know now how to extract the decision function values for a grid so you are able to plot them. Again, I am surprised how difficult it is to find it on the Internet.

With e1071, you need the following code

im=predict(K.svm,  ... ,scale=F,decision.values=T)
im=matrix(attributes(im)$decision.values,nrow=100,byrow=F)
image(seq(0, 20, length.out=100), seq(0, 20, length.out=100), im,xlab="",ylab="")
points(y)
pdf.options(reset=T)

Notice that we are extracting the value of the decision function with the attribute decision.values from the prediction SVM object. This prediction SVM object, obviously, was created with a grid of points with a range larger than the input data.

Plot classification regions in an SVM

I've been really busy. I can't really claim to be an SVM expert since my postgraduate work up to now did not deal very much with machine learning, though since late I've been working with SVM's and kernels.

In fact I have a very advanced technique that really boost the SVM performance, but that will be left until the paper is published.

One of the test I ran dealt with a modified XOR data problem. The problem consist of four groups drawn from bivariate normal distributions. They are assign two classes such that the groups of the same class are always separated from each other by some group of another class.

y=mvrnorm(50,c(3,5),Sigma=diag(c(0.5,1.5)))
y=rbind(y,mvrnorm(50,c(15,13),Sigma=diag(c(1.5,0.5))))
y=rbind(y,mvrnorm(50,c(7,5),Sigma=diag(c(0.5,1.5))))
y=rbind(y,mvrnorm(50,c(15,17),Sigma=diag(c(1.5,0.5))))
labels=c(rep(1,100),rep(-1,100))

I put the figure here so that the problem is clear, the explanation of how to get it follows.

K.svm=svm(Phi, labels, type="C", kernel="linear",probability=T)
X=as.matrix(expand.grid(list(x = seq(0, 20, length.out=100), y = seq(0, 20, length.out=100))))
# compute the kernel on X here
im=predict(K.svm, PhiX,scale=F)
im=matrix(as.numeric(im),nrow=100,byrow=F)
image(seq(0, 20, length.out=100),seq(0, 20, length.out=100),im,xlab="",ylab="",col=c("#FFFCCCFF","#FFF000FF"))#heat.colors(2))
points(y)

We see that we first train the SVM with the kernel features as explained in the previous post.
Then we create a grid spanning all the points of the region we are interested in painting and evaluate the trained SVM it there. Then we recompose the grid of classified points into a 2D plane and plot it along with the original points.

Practical kernel trick

When using an SVM, a naive usage of the software may lead us to use the standard kernels that are offered to us and perform the analysis as well as we can.

However, one can choose to build his own kernels given some prior information about the problem. For example, we could have a dataset with values of instantiated variables for each datum, but we might also have correlations telling us when those data simultaneously appear. For example, in credit rating, one may have variables such as income for each of the customers. One might also have some kind of information linking customers, such as marital status, housing status, professional category and so on. Therefore, apart from the points in $\mathbb{R}^n$ (wages), we also have some kind of coexpression graph, where two individuals score higher with each other if their categorical variables coincide, so that we build a matrix $S_{i,j}=\mbox{#}(k,j)$, where $\mbox{#}(i,j)$ is the coincidence of the vales for all variables between individual $i$ and individual $j$. We then could make a custom kernels in this way:
$$
K=K_{\sigma} + \gamma S^*
$$
Where $K_{\sigma}$ is the classical Gaussian kernel build from the euclidean distances of the wages, with parameter $\sigma$, and $S^*$ is a projection of the previous $S$ onto the cone of positive semi-definite matrices. $\gamma$ is a tuning parameter. Now, how can we use it with our SVM software?

LibSVM can be fed with the previous kernel matrix $K$ if you work with the C/C++ or the Matlab interface. Just make sure to insert an index from one to the number of data as the first column. If you are working with R with the library e1071, you won't be able to do this.

It is well known to practitioners that the kernel-trick can help us in this case. Since the Mercer's kernel expansion is
$$
k(x,y)=\sum_i \lambda_i \phi_i(x) \phi_i(y) = \Phi(x)^T \Phi(y)
$$
Where $\lambda_i$ and $\phi_i$ are the eigenvalues and eigenfunctions of the integral operator whose positive definite kernel is $k$. $\Phi(x) = [\sqrt{\lambda_1} \phi_1(x), \sqrt{\lambda_2} \phi_2(x),\ldots]$ is, therefore, the "lifting" map from $\mathbb{R}^n$ to the Hilbert space of square-summable infinite sequences $l^2$. Since we are working with a finite approximation to this operator, the eigen-approximation is given by the eigen values and eigenvectors of the kernel matrices above. We will always be able to compute these. Therefore, we substitute $ \Phi(x)^T \Phi(y)$ by its finite approximation using the eigenvalues and eigenvectors of the matrix, so that using a standard linear kernel with the lifting map (the images under the approximated $\Phi$) we achieve the values in $K$.

The following example is in R. Assume we want to compute some credit score input, and the individuals' salaries are the ones given in the variable salaries. The experience with these customers is +1 if they were good customers, and -1 if they were bad customers, and stored in the variable y.

library(e1071)
salaries=c(20000,30000,55000,50000,30000,25000)
y=c(+1,+1,+1,-1,-1,-1)

 It is clear that the experience and the salaries are not very correlated. Imagine betting is a huge problem for some segment of the population, and we have the individuals' record of using betting facilities, in this matrix

betting=c(0,0,0,1,1,1)

S=betting%*%t(betting)
Seig=eigen(S)
Seig$vectors
values=Seig$values
values[values<0]=0
D=diag(values)
Ss=Ss=U%*%D%*%t(U)
s=0.000000005
Ks=as.matrix(exp(-s*dist(salaries)^2))
diag(Ks)<-1

The gaussian distances are:  
2 0.606530660                                               
3 0.002187491 0.043936934                                   
4 0.011108997 0.135335283 0.882496903                       
5 0.606530660 1.000000000 0.043936934 0.135335283           
6 0.882496903 0.882496903 0.011108997 0.043936934 0.88249690

We now see that the true correlations between these individuals are not very strong. Therefore, supplementing Ks with S is a good idea. Ks+gamma*S is:
            1          2           3          4          5          6
1 0.000000000 0.60653066 0.002187491 0.01110900 0.60653066 0.88249690
2 0.606530660 0.00000000 0.043936934 0.13533528 1.00000000 0.88249690
3 0.002187491 0.04393693 0.000000000 0.88249690 0.04393693 0.01110900
4 0.011108997 0.13533528 0.882496903 1.00000000 1.13533528 1.04393693
5 0.606530660 1.00000000 0.043936934 1.13533528 1.00000000 1.88249690
6 0.882496903 0.88249690 0.011108997 1.04393693 1.88249690 1.00000000

We see that there is some more correlation between the betting public. Following the R program

K=Ks+gamma*S
Keig=eigen(K)
phi=Keig$vectors%*%diag(sqrt(Keig$values))

Now we train the SVM with this embedding, or lifting map, $\Phi$. Observe that if $K=VEV^T=VE_q E_q V^T = \Phi^T \Phi$ then the embedding is $\Phi^T=VE_q$, where $E_q$ has the square root of the eigenvalues.

svm(phi, y, type="C", kernel="linear", cost=1)

Subsequent customers (after retraining) will be given a more fitting credit rating in this fashion than using only their salaries.

This trick is also the basis of the celebrated technique of random features (check out the paper).

Multiclass SVM with e1071

When dealing with multi-class classification using the package e1071 for R, which encapsulates LibSVM, one faces the problem of correctly predicting values, since the predict function doesn't seem to deal effectively with this case. In fact, testing the very example that comes in the svm help (?svm on the R command line), one sees the failing performance of the function (albeit working with a correctly fitted model).

data(iris)
attach(iris)

model <- svm(Species ~ ., data = iris)

x <- subset(iris, select = -Species)
y <- Species
model <- svm(x, y, probability = TRUE)

This model is correctly fitted. However

pred <- predict(model, x)

Does not show the correct values

> pred     1      2      3      4      5      6      7      8      9     10     11     12     13     14     15 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    16     17     18     19     20     21     22     23     24     25     26     27     28     29     30 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    31     32     33     34     35     36     37     38     39     40     41     42     43     44     45 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    46     47     48     49     50     51     52     53     54     55     56     57     58     59     60 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    61     62     63     64     65     66     67     68     69     70     71     72     73     74     75 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    76     77     78     79     80     81     82     83     84     85     86     87     88     89     90 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    91     92     93     94     95     96     97     98     99    100    101    102    103    104    105 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
   106    107    108    109    110    111    112    113    114    115    116    117    118    119    120 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
   121    122    123    124    125    126    127    128    129    130    131    132    133    134    135 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
   136    137    138    139    140    141    142    143    144    145    146    147    148    149    150 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
Levels: setosa versicolor virginica

We see that every prediction is setosa, even though the dataset is equally divided between the three classes (setosa, versicolor and virginica). It seems that something went wrong with the direct prediction, but one way to overcome this problem is to use the predicted probabilities, that seem to be well computed:

pred <- predict(model, x, decision.values = TRUE, probability = TRUE)

Observe how they look

> attr(pred, "probabilities")         setosa  versicolor   virginica
1   0.980325653 0.011291686 0.008382661
2   0.972927977 0.018061300 0.009010723
3   0.979044327 0.011921880 0.009033793
...
48  0.977140826 0.013050710 0.009808464
49  0.977831001 0.013359834 0.008809164
50  0.980099521 0.011501036 0.008399444
51  0.017740468 0.954734399 0.027525133
52  0.010394167 0.973918376 0.015687457
...
97  0.009263806 0.986123276 0.004612918
98  0.008503724 0.988168405 0.003327871
99  0.025068812 0.965415124 0.009516064
100 0.007514580 0.987584706 0.004900714
101 0.012482541 0.002502134 0.985015325
...
149 0.013669944 0.017618659 0.968711397
150 0.010205071 0.140882630 0.848912299

Now we see that the most probable class is indeed the ground truth and we can correctly classify with the following function

predsvm<-function(model,newdata)
{
  prob<-attr(predict(model, newdata, probability = TRUE),"probabilities")
  n<-dim(prob)[1]
  m<-dim(prob)[2]
 
  me<-which(prob==apply(prob,1,max))
  return(as.factor(model$labels[floor((me-1)/n)+1]))
}

One might also program the following function, that deals with the way the support vector coefficients are stored in the model object, in model$coefs and model$rho:

## Linear Kernel function
K <- function(i,j) crossprod(i,j)

predsvm <- function(object, newdata) {
  ## compute start-index
  start <- c(1, cumsum(object$nSV)+1)
  start <- start[-length(start)]

  ## compute kernel values
  kernel <- sapply (1:object$tot.nSV,
                    function (x) K(object$SV[x,], newdata))

  ## compute raw prediction for classifier (i,j)
  predone <- function (i,j) {
    ## ranges for class i and j:
    ri <- start[i] : (start[i] + object$nSV[i] - 1)
    rj <- start[j] : (start[j] + object$nSV[j] - 1)
    
    ## coefs for (i,j):
    coef1 <- object$coefs[ri, j-1]
    coef2 <- object$coefs[rj, i]

    ## return raw values:
    crossprod(coef1, kernel[ri]) + crossprod(coef2, kernel[rj])
  }

  ## compute votes for all classifiers
  votes <- rep(0,object$nclasses)
  c <- 0 # rho counter
  for (i in 1 : (object$nclasses - 1))
    for (j in (i + 1) : object$nclasses)
      if (predone(i,j) > object$rho[c <- c + 1])
        votes[i] <- votes[i] + 1
      else
        votes[j] <- votes[j] + 1

  ## return winner (index with max. votes)
  object$levels[which(votes %in% max(votes))[1]]
}