Bayesian Naive Bayes for Classification with the Dirichlet Distribution

I have a classification task and was reading up on various approaches. In the specific case where all inputs are categorical, one can use “Bayesian Naïve Bayes” using the Dirichlet distribution. 

Poking through the freely available text by Barber, I found a rather detailed discussion in chapters 9 and 10, as well as example matlab code for the book, so took it upon myself to port it to R as a learning exercise.

I was not immediately familiar with the Dirichlet distribution, but in this case it appeals to the intuitive counting approach to discrete event probabilities.

In a nutshell we use the training data to learn the posterior distribution, which turns out to be counts of how often a given event occurs, grouped by class, feature and feature state.

Prediction is a case of counting events in the test vector. The more this count differs from the per-class trained counts, the lower the probability the current candidate class is a match.

Anyway, there are three files. The first is a straightforward port of Barber’s code, but this wasn’t very R-like, and in particular only seemed to handle input features with the same number of states.

I developed my own version that expects everything to be represented as factors. It is all a bit rough and ready but appears to work and there is a test/example script up here. As a bigger test I ran it on a sample  car evaluation data set from here, the confusion matrix is as follows:

testY   acc good unacc vgood

  acc    83    3    29     0

  good   16    5     0     0

  unacc  17    0   346     0

  vgood  13    0     0     6

That’s it for now. Comments/feedback appreciated. You can find me on twitter here

Links to files:

Barber Port

R amenable implementation

Example Usage

Everything in one directory (with data) here

Book Review: Applied Predictive Modeling by Max Kuhn and Kjell Johnson

This is a gem of a book.

From the introduction:

We intend this work to be a practitioner’s guide to the predictive modeling process and a place where one can come to learn about the approach and to gain intuition about the many commonly used and modern, powerful models.

…it was our goal to be as hands-on as possible, enabling the readers to reproduce the results within reasonable precision as well as being able to naturally extend the predictive modeling approach to their own data.

The book is structured into four main sections. First is General Strategies, which provides an introduction and discusses things like pre-processing and tuning. 

The next two sections cover regression and classification, each with chapters on linear and non-linear methods, as well as tree and rule based methods, with one to two chapters on practical issues such as measuring performance. 

The final section covers feature selection, predictor importance and a discussion around model performance. 

Features

There are a few things I really like:

It is not an academic or mathematical treatise; the emphasis is on practice, discussing the issues that commonly arise and how they can be approached. Plenty of references are provided for those wanting to dig deeper.

Every example has its data set and code available so one can work through the examples as presented. In most cases they are real world datasets and there is great discussion of the real world issues that arise, what should be considered and the various tradeoffs that can be made.

Discussion and code are separate. Aside from the excellent content, this is probably what I appreciate the most. Each chapter presents its content, with charts where appropriate, while the actual walk through of the code and raw output is in a separate section of the chapter. 

This makes it much easier to focus on the material being presented. It is always difficult to present source code along with discussion. This is not a book about programming per-se, it is about using existing tools to make intelligent and reasoned decisions about the task at hand. It makes a lot of sense to have the code presented separately.

Also, as far as I have read, each chart is at most only one page away from the text discussing it. This is a small thing but I feel there has been serious consideration about the presentation of the material and it has been done very well. [Update: This is not strictly true but is mostly the case throughout the book]

It is not a book about caret, the package of author Max Kuhn. To be honest I would be pretty happy even if it were about caret, which certainly does get some use in the code, but it is relatively package agnostic. 

Final thoughts

This is a great book, providing both the trees and the forest so to speak. I am unaware of any other book with similar content, and I wish I had something like this when I was first getting interested in machine learning. 

There are books that are very introductory, books that cover the details of the algorithms, and books that provide rigorous coverage of the theory, but these are not really accessible to those without a serious amount of mathematics. There are a few equations presented where appropriate, but it is certainly not the focus of the book.  

There are no real shortcomings, though if there were ever a second edition, coverage of time series methods and deep learning would be welcome. I appreciate they are both book worthy topics by themselves, and the latter is still very much a moving target. 

In summary: Great content, well written and well presented. This book would be my top recommendation to anyone looking to get started or working with predictive modeling. Well worth checking out. 

You can read more reviews on amazon here: Applied Predictive Modeling

You can find me on twitter and G+

The case for data snooping

When we are backtesting automated trading systems, accidental data snooping or look forward errors are an easy mistake to make. The nature of the error in this context is making our predictions using the data we are trying to predict. Typically, it comes from a mistake with our calculations of time offsets somewhere.

However, it can be a useful tool. If we give our system perfect forward knowledge:

1) We establish an upper bound for performance.

2) We can get a quick read if something is worth pursuing further, and

3) It can help highlight other coding errors.

The first two are pretty closely related. If our wonderful model is built using the values it is trying to predict, and still performs no better than random guessing, it’s probably not worth the effort trying to salvage it.

The flip side is when it performs well, that will be as good as it will ever get.

There are two main ways it can help identifying errors. Firstly, if our subsequent testing on non-snooped data provides comparable performance, we probably have another look ahead bug lurking somewhere.

Secondly, things like having amazing accuracy yet still performing poorly is another sign of a bug lurking somewhere.

Example

I wanted to compare SVM models when trained with actual prices vs a series of log returns, using the rolling model code I put up earlier. As a baseline, I also added in a 200 day simple moving average model.

(S) Indicates snooped data

A few things strike me about this.

For the SMA system, peeking ahead by a day only provides a small increase in accuracy. Given the longer-term nature of the 200 day SMA this is probably to be expected.

For the SVM trained systems, the results are somewhat contradictory.

For the look forward models, training on price data had much lower accuracy than the log returns, and the log return model performed much better. Note that both could have achieved 100% accuracy by predicting its first column of training data.

However, when not snooping, the models trained on closing prices did much better than those trained on returns. I’m not 100% sure there isn’t still some bug lurking somewhere, but hey if the code was off it would’ve shown up in the forward tested results no?

Feel free to take a look and have a play around with the code, which is up here.

Tracking down errors in R

It's that moment we all know and love, somewhere in our code something has gone wrong. We think we have done everything right, but instead of expected glory we find only terse red text lain below our lintel. 

This can be very frustrating, and trouble shooting these issues can often be very time consuming. 

All is not lost. There are a few bits of R that can greatly help finding out what exactly has gone wrong and where, which in turn should suggest a reasonable course of action.

First we will look at some simple methods we can use to track down issues, namely the warn option and traceback, and then we will look at stepping through functions with debug. 

First Steps

I'm going to use an example that has stuck with me from when I first started using R, using neural nets for classification with the iris data.

Let's take a look at the error:

library(nnet)

X <- iris[,1:4]

Y <- iris[,5]

mod <- nnet(X, Y, size=2)

# weights: 19

Error in nnet.default(X, Y, size = 2) : 

  NA/NaN/Inf in foreign function call (arg 2)

In addition: Warning message:

In nnet.default(X, Y, size = 2) : NAs introduced by coercion

Urgh. Really? What the hell does that mean?

We can use the built in traceback() function to see where this error occured

> traceback()

2: nnet.default(X, Y, size = 2)

1: nnet(X, Y, size = 2)

> 

We can see our call to nnet(), which in turn has called nnet.default() and which is where our error has come from. 

In the error output, we can see there was also a warning "NAs introduced by coercion". As we weren't expecting any warnings, let's track down that, as errors tend to compound. 

Warnings

To find out where that message was coming from, we will use options(warn = 2) which will turn warning messages into errors. We can do this by setting the warn option to a specific level, in this case 2. 

The default is warn = 0, which means warnings will be stored until the top level function returns. We could use warn = 1 which wil print the warning as it is encountered, but in this case we want to stop straight away, so will set it to 2.

options(warn = 2)

Let's try again:

> mod <- nnet(X, Y, size=2)

# weights: 19

Error in nnet.default(X, Y, size = 2) : 

  (converted from warning) NAs introduced by coercion

>

Hmm, still coming from nnet.default, let's see if traceback() is offering any new information.

> traceback()

6: doWithOneRestart(return(expr), restart)

5: withOneRestart(expr, restarts[[1L]])

4: withRestarts({

       .Internal(.signalCondition(simpleWarning(msg, call), msg, 

           call))

       .Internal(.dfltWarn(msg, call))

   }, muffleWarning = function() NULL)

3: .signalSimpleWarning("NAs introduced by coercion", quote(nnet.default(X, 

       Y, size = 2)))

2: nnet.default(X, Y, size = 2)

1: nnet(X, Y, size = 2)

> 

We see a whole bunch of extra stuff in the traceback(), but from 3 onwards it appears to be sideffects of having set warn = 2. We do however see our warning has come from nnet.default again, so we will dig into that to see if we can find out what is going on.

Getting dirty with debug()

To do this, we can use the debug function. We will turn on debugging for nnet, which will let us step through the code line by line as it is executed.

> debug(nnet)

> mod <- nnet(X, Y, size=2)

debugging in: nnet(X, Y, size = 2)

debug: UseMethod("nnet")

Browse[2]> 

The Browse> prompt tells us we are in the debugger. The debug: UseMethod("nnet") tells us the next line of code to be executed is UseMethod("nnet"). We could enter 'n' here to continue to the next line, however a convenient default is just hitting enter (i.e. an empty line). 

Browse[2]> 

debugging in: nnet.default(X, Y, size = 2)

debug: {

    net <- NULL

...

    class(net) <- "nnet"

    net

}

Browse[3]> 

Here R has printed out the R source for the function we just entered. We can see we are still at our Browse> prompt, so will continue on by hitting enter again and again:

Browse[3]> 

debug: net <- NULL

Browse[3]> 

debug: x <- as.matrix(x)

Browse[3]> 

debug: y <- as.matrix(y)

Browse[3]> 

debug: if (any(is.na(x))) stop("missing values in 'x'")

Browse[3]> 

debug: NULL

Browse[3]> 

Well this is neat. Hitting enter at the prompt, R shows us each line that is about to be executed. We will continue on hitting enter until we see our error message.

Browse[3]> 

… 

Browse[3]>

debug: if (length(weights) != ntr || any(weights < 0)) stop("invalid weights vector")

Browse[3]> 

debug: NULL

Browse[3]>

One thing of interest, we can see a conditional if statement is about to be run. When the conditional is evaluated as false, meaning the conditioning code won't be executed, we will see the debug: NULL printed out. 

Browse[3]> 

debug: Z <- as.double(cbind(x, y))

Browse[3]> 

Error in nnet.default(X, Y, size = 2) : 

  (converted from warning) NAs introduced by coercion

> 

Well, there is our warning message. Unfortunately we have lost the Browse> prompt, meaning we are no longer inside the function being debugged, but back at the main prompt. 

This is a side effect of our rather aggressive options("warn") setting. Let's tone it down a bit and set it to 1 so the warnings will be printed as they occur, then jump back into debugging.

> options(warn=1)

> mod <- nnet(X, Y, size=2)

debugging in: nnet(X, Y, size = 2)

debug: UseMethod("nnet")

Browse[2]> 

debugging in: nnet.default(X, Y, size = 2)

debug: {

    net <- NULL

…

Browse[3]> 

…

Browse[3]> 

debug: Z <- as.double(cbind(x, y))

Browse[3]> 

Warning in nnet.default(X, Y, size = 2) : NAs introduced by coercion

debug: storage.mode(weights) <- "double"

Browse[3]> 

After some time, we get our warning message, and we are still in our debugged function. The warning is coming from the line

Z <- as.double(cbind(x, y))

What could be the problem here? Something is going wrong when nnet is converting x and y to doubles. Let's take a look at them to see if there is anything going on.

We can do this by using get("variable"), where variable is the quoted name of the variable. First let's take a look at other two variables in the function to see how it works:

Browse[3]> get("nout")

[1] 1

Browse[3]> get("ntr")

[1] 150

Browse[3]> 

Looking back at the debug output of each line, we can see they were set as some of the dimension info for X and Y, nout has a value of 1 and ntr has a value of 150. 

Let's take a look at x and y now:

Browse[3]> get("x")

       Sepal.Length Sepal.Width Petal.Length Petal.Width

  [1,] 5.1 3.5 1.4 0.2

  [2,] 4.9 3.0 1.4 0.2

  [3,] 4.7 3.2 1.3 0.2

  [4,] 4.6 3.1 1.5 0.2

…

Browse[3]>

This is our input X data.

Browse[3]> get("y")

       [,1]        

  [1,] "setosa"    

  [2,] "setosa"    

  [3,] "setosa"    

  [4,] "setosa"    

...

Browse[3]>

And here is our Y data, just as we passed in.

Recall the warning was triggered by this line:

Z <- as.double(cbind(x, y))

Which is converting x and y to the numeric double type. 

We can see that x is numeric data, while y is character data. What happens if you convert character data to numeric? It doesn't seem to make sense, but let's try:

Browse[3]> as.numeric(get("y"))

Warning: NAs introduced by coercion

  [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

 [37] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

 [73] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[109] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[145] NA NA NA NA NA NA

Browse[3]> 

A-ha!

This is just the warning we see, and thinking about it, we can understand why converting strings to numeric data is probably not going to be particularly meaningful.

Let's get out of the debugger and think about what might be going on. Hit 'c' to continue execution 

Browse[3]> c

Error in nnet.default(X, Y, size = 2) : 

  NA/NaN/Inf in foreign function call (arg 2)

>

We are now back at the main prompt. 

Now what??

What is going on? It seems reasonable that we should be able to pass a factor for classification, in fact we are pretty sure that's what we saw being used in the examples in the nnet package documentation.

A careful reading of help(nnet) reveals some details. In particular:

If the response in formula is a factor, an appropriate classification network is constructed ... If the response is not a factor, it is passed on unchanged to nnet.default.

It is possible to pass a factor in Y, but we must use the formula syntax. Looking at the examples, we see the matrix syntax in use as well, however it is transforming the y values using class.ind(). 

Right about now I contemplate calling the police, as I'm pretty sure someone has snuck in and changed the docs while I wasn't looking. There is no way I would make such a simple mistake ...  

Anyway, let's turn off debugging, and see if we can get this working:

> undebug(nnet)

> mod <- nnet(X, class.ind(Y), size=2)

> mod <- nnet(Y~., data=cbind(X, Y), size=2)

Both methods run without error (or warning). Success at last.

Outro

We've seen a few ways we can dig into R and track down where things are going wrong. First is using options(warn=2) to make the R convert warnings to errors, and using traceback() to find out in which function the issue is arising. 

Often, this may be enough to get things back on track, especially if the function causing trouble is small. For more complex issues, we can use debug() which will let us step through the function line by line, inspecting variables and internal state as needed.

It should be said this will only be of use when the source is available to R, i.e. the function in question has been implemented in R and not compiled C/C++/Fortran or whatever that is imported as a shared lib. Deubgging that is more involved and there is some reasonable documentation on this available here.

That's all for now, happy hunting and thanks for stopping by.

Multidimensional Scaling and Company Similarity

Background and idea

Often we are looking at a particular sector, and want to get a quick overview of a group of companies relative to one another. I thought I might apply Multidimensional Scaling (MDS) to various financial ratios and see if it gave us anything useful.

The premise is that companies in similar industries should all have a degree of sameness, so MDS might be useful to highlight the companies that stand out from the crowd, perhaps in some literal sense ...

Method

I mostly use the data functions from quantmod to retrieve the financial statements from Google Finance. As always with free data, the quality is variable, but good enough for our purpose today. We need to do a bit of dancing to get the market price at the time the results were released, and this uses data from Yahoo Finance. It was a little bit more work to implement, but worth it so we can include P/E in the comparison.

I looked at two groups of companies, tech stocks and financials/banks.

For the tech stocks I used ROE, EPS, P/E, Operating Margin, Current Ratio, Gearing, Asset Turnover and Debt Ratio. For the financials, I used ROE, EPS, P/E, Gearing and Debt Ratio, mainly because the data available did not have the line items required to calculate the other ratios. 

The data from Google gives the last four periods, with the most recent coming first. It also gives Annual and Quarterly data and the charts below use the annual results. Annual Period 1 means the most recent results. Due to the scaling function, the actual scales on the graphs are not particularly meaningful, so I took them out. 

Charts

These are the charts for the most recent results (so end of year 2011). Overall, I am quite pleased with the results. We can see how most of the companies cluster together, while a few seem to be quite different. This shows at a glance the companies that might be worthy of further investigation. 

Tech Stocks

Financials

Outro

Code is up here MDS Company Similarity with R, it should hopefully be documented enough for others to mess around with. Any questions, comments or suggestions are very much appreciated as always.

As an aside, this is the first R program I wrote devoid of any for loops. I finally feel I am coming to grips with the language.

Another cut at market randomness

I have some background in computer security and one day found myself tasked with assessing the quality of randomness for session id tokens generated by popular web frameworks (namely Java and .NET). As it turns out, NIST have developed a series of tests for just this purpose detailed here.

As a non-believer in the absolute randomness of markets, I thought I might take a look at the series of returns from SPX/GSPC to see how they held up.

About the tests

The NIST tests are designed to be run on bit streams from a random number generator, so some liberties had to be taken. I am mainly interested in the predictability of up days vs down days, so I first generated a series of returns using quantmod then encoded them using a 1 for an up day and -1 for a down day and run the tests on the encoded stream. The NIST tests use -1 in place of a 0.

The main test is called the monobit test, to quote from the document:

The focus of the test is the proportion of zeroes and ones for the entire sequence.  The purpose of this test is to determine whether the number of ones and zeros in a sequence are approximately the same as would be expected for a truly random sequence.  The test assesses the closeness of the fraction of ones to ½, that is, the number of ones and zeroes in a sequence should be about the same.  

The tests are done using a significance level of 0.01, so come with a good degree in confidence assuming the underlying method is sound. 

One caveat is the length of the runs compared and how it relates to the distributions used to model the results. For the monobit test, the suggested input size is 100, requiring 101 days of data to determine 100 up or down days. If we were looking at 32 bit integers, 100 bits would only be 3 "full" random numbers, so arguably we would want to look at shorter time periods (e.g. 3-5 days of data). Given the difficulties around distributions which require a large n, I thought I would vary the significance level instead of a lower n, as our requirements are not as stringent as those for cryptographic random numbers.   

Results

At a basic level, this series does appear to be random, at least the vast majority of the time with n = 100 and alpha = 0.01. My confirmation bias was very upset with this. 

However, if we plot the proportion of runs deemed random vs the significance level, we see the proportion rising as one might expect. One thing that remains unexplained is why this appears to rise in steps rather than something more linear, though I expect this to be a side affect of either methodology or the normalisation done by the tests. I also took a look at the weekly data, which tends to a greater proportion of non-random runs quicker than daily data.

I am interested in the applications of machine learning to financial markets. Close to close returns that we have been looking at here are not the only information we have available, nor are they what I trade on a personal level. Also this is only one price series, and one could argue that in practise it is not actually a tradable series. 

Close to close returns are very useful in lots of applications, but if we are trying to build some predictive model we might need to look for more predictable pastures. Machine learning algorithms are great, but can only do so much. Finding some better potential inputs is what I will take a look at next.

The test code is available on github here: R Monobit test. Would be very interested to hear if anyone else takes a look.

I also took a visual and binomial look at randomness of the series in this post: A visual look at market randomness.

Oh and in case you were wondering, the web session tokens all turned out to be very strong.