Mechanistician

Methods for comparing rankings of search engine results

2011-09-29T15:00:00.000-07:00

Interesting visualization from a paper on comparing search engine results:

InMaps - visualizing your LinkedIn network

2011-08-05T20:31:00.000-07:00

LinkedIn's InMaps does a pretty good job at detecting communities within your connections. It requires that you have at least 50 connections to use it (otherwise there is not enough info to make it interesting), but that threshold might still be too small, depending on your career trajectory. For me, it got really good once I had 75+ connections:

Search, Network, Analytics @ LinkedIn

2011-07-02T15:37:00.000-07:00

Tons of interesting projects from LinkedIn: SNA.

For those (still?) not familiar with the whole NoSQL movement, here is a nice intro, from a data modeling perspective: No Relation - The Mixed Blessings of Non-Relational Databases.

I'm particularly fond of the approach taken in HBase, based on the BigTable paper from Google.

As if that wasn't enough, here's something else to sink another couple of hours into: Challenges in Building Large-Scale Information Retrieval Systems.

Mining of massive datasets

2011-06-17T15:52:00.000-07:00

Good resource for mining large datasets: book and class. MinHash, a technique used in the google news personalization system is described in detail in the book.

TED talk on issues with online personalization

2011-05-15T11:36:00.000-07:00

IBM's Watson

2011-03-11T11:45:00.001-08:00

To build your own watson, see this article.

Scala

2011-01-25T20:43:00.000-08:00

I've been playing with Scala lately, and I have to admit it is growing on me. Below is a simple recursive interpreter written in it for a small Turing-complete language. Check out how it keeps state as it executes (the 'state0' and 'update' methods, with 'update' returning the head function of a linked list of functions as new variables are assigned to, and 'state0' being the tail function, which returns 0 as a default for all unassigned variables):


abstract class Aexp
case class IntAexp(i: Int) extends Aexp
case class VarAexp(v: String) extends Aexp
case class SumAexp(a1: Aexp, a2: Aexp) extends Aexp
case class MulAexp(a1: Aexp, a2: Aexp) extends Aexp

abstract class Bexp
case class BoolBexp(bool: Boolean) extends Bexp
case class LessBexp(a1: Aexp, a2: Aexp) extends Bexp
case class AndBexp(b1: Bexp, b2: Bexp) extends Bexp

abstract class Cmd
case class SkipCmd() extends Cmd
case class AssignCmd(v: String, a: Aexp) extends Cmd
case class SeqCmd(c1: Cmd, c2: Cmd) extends Cmd
case class IfCmd(b: Bexp, c1: Cmd, c2: Cmd) extends Cmd
case class WhileCmd(b: Bexp, c: Cmd) extends Cmd

object eval {
  type State = (String) => Int

  def state0(v: String): Int = 0

  def update(s: State, v: String, n: Int): State = { 
   v2 => { if (v.equals(v2)) n else s(v2) }
  }

  def A(a: Aexp, s: State): Int = {
    a match {
      case IntAexp(i) => i
      case VarAexp(v) => s(v)
      case SumAexp(a1, a2) => A(a1, s) + A(a2, s)
      case MulAexp(a1, a2) => A(a1, s) * A(a2, s)
    }
  }

  def B(b: Bexp, s: State): Boolean = {
    b match {
      case BoolBexp(bool) => bool
      case LessBexp(a1, a2) => A(a1, s) < A(a2, s)
      case AndBexp(b1, b2) => B(b1, s) && B(b2, s)
    }
  }

  def bigstep(c: Cmd, s: State): State = {
    c match {
      case AssignCmd(v, a) => update(s, v, A(a, s))
      case SkipCmd() => s
      case SeqCmd(c1, c2) => bigstep(c2, bigstep(c1, s))
      case IfCmd(b, c1, c2) => if (B(b, s)) bigstep(c1, s) else bigstep(c2, s)
      case WhileCmd(b, c) => if (B(b, s)) bigstep(WhileCmd(b, c), bigstep(c, s)) else s
    }
  }

  def main(args: Array[String]): Unit = {
    var c: Cmd = WhileCmd(LessBexp(VarAexp("x"), IntAexp(10)),
      AssignCmd("x", SumAexp(VarAexp("x"), IntAexp(1))))

    System.out.println("Big-step evaluating: " + c)
    val res: State = bigstep(c, state0)
    System.out.println("Result: " + res("x") + "\n\n")
  }
}

Talk on Ensemble methods

2011-01-10T14:57:00.000-08:00

Upcoming interesting talk at the local Bay Area Artificial Intelligence Meetup Group: Diversity, Complexity, and Regularization in Ensemble Models.

Gamma database machine

2010-12-13T13:46:00.001-08:00

While the Titanium DBMS may not have gone anywhere, here's an example of an older system that was certainly influential in today's hot technologies: the Gamma database machine. Contrast Gamma with a combination of Hadoop map-reduce, Pig, and HBase. Another interesting system along those lines is HadoopDB.

Blast from the past: Titanium DBMS

2010-12-13T13:22:00.001-08:00

Came across an old article from Dr. Dobb's while looking for something completely unrelated. It talks about some of the challenges faced by the developers of the system back then known as the NTC Ship Manager, then SafeNet, and now known as NS5, an Enterprise Resource Planning (ERP) system for the maritime shipping industry which I worked on some years ago (NS5, today, uses Java and MySQL). Here is an article about Titanium form 1993, claiming it outpeforms "traditional relational databases" by 5 to 1:

CrowdConf 2010

2010-10-25T18:24:00.001-07:00

Videos of the conference are up.

Adding a new operator to PostgreSQL

2010-10-24T08:53:00.000-07:00

So, I'm adding a new operator to PostgreSQL as part of a project for my database systems implementation class in order to enhance pg with the functionality described in the paper titled A Scalable, Predictable Join Operator for Highly Concurrent Data Warehouses.

Fortunately pg uses Flex/Bison to manage parsing SQL, greatly simplifying the task of making modifications to the grammar.

One of the steps for doing so includes adding the new keyword to the file pgsql/src/include/parser/kwlist.h, an excerpt of which is shown here:


PG_KEYWORD("isnull", ISNULL, TYPE_FUNC_NAME_KEYWORD)
PG_KEYWORD("isolation", ISOLATION, UNRESERVED_KEYWORD)
PG_KEYWORD("join", JOIN, TYPE_FUNC_NAME_KEYWORD)
PG_KEYWORD("key", KEY, UNRESERVED_KEYWORD)

The new operator described in the paper is aptly named CJOIN, so naturally, I defined my keyword right below the existing JOIN keyword:


PG_KEYWORD("isnull", ISNULL, TYPE_FUNC_NAME_KEYWORD)
PG_KEYWORD("isolation", ISOLATION, UNRESERVED_KEYWORD)
PG_KEYWORD("join", JOIN, TYPE_FUNC_NAME_KEYWORD)
PG_KEYWORD("cjoin", CJOIN, TYPE_FUNC_NAME_KEYWORD)
PG_KEYWORD("key", KEY, UNRESERVED_KEYWORD)

Only to find out my changes weren't being recognized (i.e. I would issue a cjoin query from the client and the server would still complain of a syntax error). Why? Where in the near 1 million lines of code had I missed something?

Well, it turns out the fact that the list in kwlist.h was in alphabetical order was not by coincidence or due to a desire to be neat, it actually affected execution, as explained in the comments at the beginning of that file:


/*
* List of keyword (name, token-value, category) entries.
*
* !!WARNING!!: This list must be sorted by ASCII name, because binary
*   search is used to locate entries.
*/

D'oh!

How to start a social movement

2010-10-15T17:47:00.000-07:00

Interesting video on the dynamics of starting a social movement:

Crowdsourcing

2010-09-30T21:26:00.000-07:00

This Monday (10/4) is crowdconf, the "world's first conference on the future of distributed work."

Here are some links for papers/research on crowdsourcing:

MapReduce algorithm design book

2010-08-05T22:10:00.000-07:00

Jimmy Lin's book is a great resource for people writing map-reduce algorithms (the final pre-production manuscript is available as a free download).

Architecture for Highly Concurrent Server Apps

2010-07-26T02:45:00.000-07:00

Came across this paper today, SEDA: An Architecture for Well-Conditioned, Scalable Internet Services. From the author's website:

"SEDA is an acronym for staged event-driven architecture, and decomposes a complex, event-driven application into a set of stages connected by queues. This design avoids the high overhead associated with thread-based concurrency models, and decouples event and thread scheduling from application logic. By performing admission control on each event queue, the service can be well-conditioned to load, preventing resources from being overcommitted when demand exceeds service capacity. SEDA employs dynamic control to automatically tune runtime parameters (such as the scheduling parameters of each stage), as well as to manage load, for example, by performing adaptive load shedding. Decomposing services into a set of stages also enables modularity and code reuse, as well as the development of debugging tools for complex event-driven applications."

Even though SEDA hasn't been worked on since 2002 (you can download the code from the website), it presents some very interesting ideas. For completeness, here's one paper (Why Events Are A Bad Idea (for high-concurrency servers)) that critiques SEDA's performance, while this other one (Comparing the Performance of Web Server Architectures) presents findings in favor of a SEDA-like architecture.

In the Java world, the event-based model for building highly scalable server applications is preferred over the simpler thread-per-connection model, and this article (Architecture of a highly scalable NIO-based server) provides a good intro. Netty, an open-source project built around the concepts presented in the article and in SEDA, provides "an asynchronous event-driven network application framework and tools for rapid development of maintainable high performance & high scalability protocol servers & clients."

Update 1: This page also has a lot of good info on the topic.
Update 2: This post also presents some interesting findings.

Large-Scale Machine Learning with Mahout - Part 3

2010-07-07T21:50:00.001-07:00

RecommenderJob essentially coordinates a series of map-reduce tasks to accomplish the goal of producing recommendations for the specified users from the specified training data. It performs the computation in a batch fashion and has no support for incremental processing. For example, considering the GroupLens dataset mentioned in the previous post, if suddenly you received another 100k more ratings, RecommenderJob would have to process the accumulated 200k ratings in one shot to compute new recommendations based on the latest training data. Despite this and other limitations, it illustrates Hadoop's map-reduce paradigm well.

RecommenderJob contains 4 map-reduce pairs, in the following order:

ItemIDIndexMapper and ItemIDIndexReducer
ToItemPrefsMapper and ToUserVectorReducer
UserVectorToCoocurrenceMapper and UserVectorToCoocurrenceReducer
RecommenderMapper and IdentityReducer

You can think of the data flowing from mapper to reducer, and from map-reduce pair to map-reduce pair.

The first map-reduce pair is configured as follows:


    JobConf itemIDIndexConf = AbstractJob.prepareJobConf(inputPath, itemIDIndexPath, jarFile,
      TextInputFormat.class, ItemIDIndexMapper.class, IntWritable.class, LongWritable.class,
      ItemIDIndexReducer.class, IntWritable.class, LongWritable.class, MapFileOutputFormat.class);

Where inputPath is a directory containing one or more files, each line in a file being a record of interest. inputPath denotes the input is for the mapper. ItemIDIndexPath is the path to the directory where the output of the reducer will be saved. TexInputFormat is a class that handles a specific file format and will be used to process the files in inputPath. ItemIDIndexMapper is the map class. IntWritable and LongWritable are the output key and value types for the mapper. ItemIDIndexReducer is the reducer class. The second IntWritable and LongWritable classes represent the format of the key and value output types of the reducer. MapFileOutputFormat is the class the handles the format that the output of the reducer will be in.

Once this first job is configured, it is run by the JobClient. The goal of the first map-reduce pair in this job is to handle the possibility that the spectrum of itemID's falls outside the range of java's int. The reason we would like to map itemID's to ints is that data structures in the remaining tasks contained in the job (arrays) use ints as the indexing method. So, the job will hash all (long) userIDs to ints (where there will inevitably be collisions, which in this case are unhandled) and save that map in the output file, which not coincidentally is of type MapFileOutputFormat. The possibility that the range of int would not be sufficient to uniquely identify all itemID's is something that would be handled in a more robust recommender.

Note that the hash function, idToIndex, has been corrected in Mahout's trunk so as to not return negative values:


static int idToIndex(long itemID) {
  return 0x7FFFFFFF & ((int) itemID ^ (int) (itemID >>> 32));
}

ItemIDIndexMapper's map method is called once per line in the input file. Each time, the key will be the line number and the value will be the line. This is so because that is what TexInputFormat, which we specified, does. In turn, the output will be another set of key-value pairs, the key being the hash value of the itemID and the value being the itemID.

What happens next, the magic that Hadoop performs between the mapper and the reducer and which is commonly known as the "shuffle," is complicated. Tom White's book has a great explanation. Another good and concise explanation can be found on page 3 of this paper, which discusses modifications to Hadoop so as to enable data to be pipelined between tasks and between jobs.

Once it comes time to run ItemIDIndexReducer's reduce method, Hadoop has gathered all values with the same key into a list. The reduce method will then choose the smallest of those values and output it to be the corresponding value (itemID) for the key (hashed itemID). This gets saved into a file-based map-like data structure, MapFile, in a directory specified by ItemIDIndexPath.

The second map-reduce task is configured as follows:


    JobConf toUserVectorConf = AbstractJob.prepareJobConf(inputPath, userVectorPath, jarFile,
      TextInputFormat.class, ToItemPrefsMapper.class, LongWritable.class, ItemPrefWritable.class,
      ToUserVectorReducer.class, LongWritable.class, VectorWritable.class, SequenceFileOutputFormat.class);

Its purpose is to is to create the user preference vector mentioned in the previous post. It uses the same input as the previous map-reduce pair and the reducer outputs a file composed of a sequence of key-vector instances. This SequenceFile is distinct from the MapFile used in the previous task only slightly. The MapFile can be thought of as SequenceFile that has an index that allows random access lookups by key. The SequenceFile, on the other hand, as its name suggests, is meant to be accessed sequentially.

Notice that ToUserVectorReducer hashes the itemID using the same method as ItemIDIndexMapper. This is because the user preference vector will be stored in a RandomAccessSparseVector, which is a HashMap-like structure, where keys are ints and values are doubles. I have looked to see if there is a RandomAccessSparseVector in Mahout that allows for long keys, but one sure would be handy.

As we follow along in ToUserVectorReducer's reduce method, we see that the user preference vector is limited to a configurable number of preferences, and by default that limit is set to 20. This surfaces another area that can be improved in this current implementation, which even if it allowed for itemID's to be longs, it does not account for a user preference vector that does not fit in RAM. I do not mean these observations as criticism of Mahout, which already in its early stages has a tremendous amount of useful functionality; they are simply meant to get your mind thinking about some of the issues you could solve should you venture into building your own recommendation system or into contributing to Mahout itself. This size limitation of the user preference vector is restricted to non-zero values, since the RandomAccessSparseVector is initialized to a cardinality of Integer.MAX_VALUE. The cardinality is the 'effective size' of the RandomAccessSparseVector, meaning that it represents the number of allowable indices into the vector. RandomAccessSparseVector is useful for storing sparse vectors, as it only keeps track of non-default (zero) values stored in it.

The method that it uses to limit the size of the user preference vector won't exactly cut its length at the limit, but it does a decent job at picking a user's top preferred items. On the other hand, this kind of scheme would not necessarily to well work if your problem domain denotes preferences as binary values, such as the one mentioned in the Google news personalization system described in this paper (and you can also watch a Google talk video about it here).

The third map-reduce task is configured as follows:


    JobConf toCooccurrenceConf = AbstractJob.prepareJobConf(userVectorPath, cooccurrencePath, jarFile,
      SequenceFileInputFormat.class, UserVectorToCooccurrenceMapper.class, IntWritable.class,
      IntWritable.class, UserVectorToCooccurrenceReducer.class, IntWritable.class, VectorWritable.class,
      MapFileOutputFormat.class);

Its purpose is to build the other major data structure RecommenderJob computes: the coocurrence matrix. This job will read the user preference vector file computed in the previous job and save the computed coocurrence matrix to the cooccurrencePath directory in the form of a MapFile.

UserVectorToCooccurrenceMapper's map method is called once for each user preference vector, and it will output every permutation of 2 itemID's for that user (excluding pairs of same itemID's), in other words, the items that co-occur in a user's preference list. Note that it will output seemingly repeated pairs, for example, it will output [Key:itemID=123, value:itemID=456] and [Key:itemID=456, value:itemID=123]. This is needed because we are gonna build the rows (or columns, since the matrix is symmetrical) of the coocurrence matrix and we need to access each row (or column) for every itemID. The Hadoop will then collect all values for a given key, and feed the key and list of values into the reducer.

UserVectorToCooccurrenceReducer's reduce method will build one RandomAccessSparseVector for every itemID (hashed itemID, that is); the rows (or columns) of the matrix. it simply iterates over the items in cooccurrenceRow and keeps track of how many times it sees a given item. It also has a small optimization where it zeros out any coocurrence of 1. Finally each cooccurrenceRow will be saved into a MapFile, keyed by itemID, so it can be quickly accessed in random order.

The fourth, and last, map-reduce task is configured as follows:


    JobConf recommenderConf = AbstractJob.prepareJobConf(userVectorPath, outputPath, jarFile,
      SequenceFileInputFormat.class, RecommenderMapper.class, LongWritable.class,
      RecommendedItemsWritable.class, IdentityReducer.class, LongWritable.class,
      RecommendedItemsWritable.class, TextOutputFormat.class);
    recommenderConf.set(RecommenderMapper.COOCCURRENCE_PATH, cooccurrencePath);
    recommenderConf.set(RecommenderMapper.ITEMID_INDEX_PATH, itemIDIndexPath);
    recommenderConf.setInt(RecommenderMapper.RECOMMENDATIONS_PER_USER, recommendationsPerUser);
    recommenderConf.set(RecommenderMapper.USERS_FILE, usersFile);
    recommenderConf.setClass("mapred.output.compression.codec", GzipCodec.class, CompressionCodec.class);

The additional job configuration settings simply make available information such as paths of previously computed output files that this job will need and other settings such as how many recommendations are to made to each user. It is also specified that the output of this job is to be compressed.

The reducer specified in this job is the IdentityReducer, which simply outputs its input (notice it is deprecated, as there is a newer API that simplifies things a bit - the concepts, however, remain the same in the newer API):


@Deprecated
public class IdentityReducer<K, V>
    extends MapReduceBase implements Reducer<K, V, K, V> {

  /** Writes all keys and values directly to output. */
  public void reduce(K key, Iterator<V> values,
                     OutputCollector<K, V> output, Reporter reporter)
    throws IOException {
    while (values.hasNext()) {
      output.collect(key, values.next());
    }
  }
}

RecommenderMapper is a bit more involved than the tasks we've seen so far. In addition to overriding the map method, it also overrides configure and close. The calling of these methods is orchestrated by org.apache.hadoop.mapred.MapRunner, which will call configure upon instantiating RecommenderMapper, and will call close once it has called its map method for every key-value pair assigned to it.

RecommenderMapper will perform the matrix-vector multiplication of the coocurrence matrix with a user's preference vector to compute the recommendation vector mentioned in the previous post. It will then select the highest valued indices (as many as were specified in RECOMMENDATIONS_PER_USER, which defaults to 10) of the recommendation vector (filtering them so as to not select something the user has already seen) and them map those indices back to the original itemIDs. The recommendation are encapsulated on a RecommendedItemsWritable object, which together with the LongWritable which encapsulates the userID, are accumulated in the output file.

Notable observations include how the matrix-vector multiplication is done (keep a running sum of the result of multiplying the ith column of the matrix by the ith value of the vector) and the mechanics of how the columns of the coocurrence matrix are loaded and cached (which is complex enough to deserve its own post).

Another thing that should jump at you, is the file that is specified as input into the last task, which goes to the RecommenderMapper. It is the file containing the user preference vectors, and thus if you were specifying multiple mappers, it would be based on that information that input splits would be performed. Whereas ideally, if you had many users for who to make recommendations, you would want to distribute based on that criteria instead. Similarly, when you have a small number of users, you would just want to perform one matrix - vector multiplication for that user and save yourself the trouble of calling the RecommenderMapper's map method for users to whom you are not interested in making recommendations. This is handled in the current implementation by returning from the map method if the user preference vector does not belong to a user we are intrested in:


  @Override
  public void map(LongWritable userID,
                  VectorWritable vectorWritable,
                  OutputCollector<LongWritable,RecommendedItemsWritable> output,
                  Reporter reporter) throws IOException {
    
    if ((usersToRecommendFor != null) && !usersToRecommendFor.contains(userID.get())) {
      return;
    }
    ...
  }

Large-Scale Machine Learning with Mahout - Part 2

2010-07-06T21:29:00.001-07:00

Following up on the previous post, today we examine a map-reduce-based job that does a bit more than count how many times each word occurs in a document. You'll wanna already be familiar with Hadoop and have it installed and running on your computer, and for that you can use the references described previously.

RecommenderJob runs a collaborative filtering algorithm that is fairly easy to distribute compared with others (such as the predominantly online SVD algorithms popularized in the recent Netflix challenge, see here and here). You can check Mahout In Action for a more detailed description, but briefly, 2 main data structures are computed:

The coocurrence matrix
The user preference vector

The coocurrence matrix is a symmetric matrix of items by items, where the intersection of an item row and an item column holds the count of how many times that item pair was present in the set of rated items of a given user. The user preference vector is simply a vector of all the items, each user having its own preference vector. The only non-zero items in a user's preference vector are those which the user has rated, and the rating is the value for that item in the vector.

To compute recommendations for a given user, the coocurrence matrix is multiplied on the right by the user preference vector, and the result is another vector with the same dimensions as the user preference vector, let's call it the "recommendation vector." What you'll wanna do then is recommend to the user the highest-valued items in their recommendation vector (presumably excluding those items the user has already rated).

One simple dataset to use is the 100k GroupLens dataset, which can be downloaded from here, and which has the following tab-separated format: user id | item id | rating | timestamp. RecommenderJob expects a comma-separated file and you can either write a script to convert it to the format expected by the RecommenderJob or edit the the RecommenderJob to process that format (you could also create a map-reduce task with an IdentityReducer to do the preprocessing).

For starters, the job can be run from the IDE as a junit unit test, similar to the one shown below:


package dr.item_rec;

import org.apache.hadoop.fs.FileSystem;
import org.apache.hadoop.fs.Path;
import org.apache.hadoop.mapred.JobConf;
import org.junit.Test;
import org.junit.Assert;

public class RecommenderJob_Tester {
    @Test
    public void MainTest() throws Exception {
        JobConf conf = new JobConf();
        conf.set("fs.default.name", "file:///");
        conf.set("mapred.job.tracker", "local");

        Path input = new Path("/Users/mario/Desktop/fastcf/data/rec_job");
        Path output = new Path("/Users/mario/Desktop/fastcf/data/rec_job_out");
        Path usersFile = new Path("/Users/mario/Desktop/fastcf/data/recs/users");

        FileSystem fs = FileSystem.getLocal(conf);
        // delete output dirs, including those defined in AbstractJob,
        // RecommenderJob's parent class
        fs.delete(output, true);
        fs.delete(new Path("/Users/mario/Desktop/fastcf/distributed_recommender/temp/userVectors"), true);
        fs.delete(new Path("/Users/mario/Desktop/fastcf/distributed_recommender/temp/itemIDIndex"), true);
        fs.delete(new Path("/Users/mario/Desktop/fastcf/distributed_recommender/temp/cooccurrence"), true);

        RecommenderJob driver = new RecommenderJob();
        driver.setConf(conf);

        int exitCode = driver.run(new String[]{
                "--input",
                input.toString(),
                "--output",
                output.toString(),
                "--usersFile",
                usersFile.toString()
        });
        
        Assert.assertEquals(0, exitCode);
    }
}

The source for the project can be downloaded from here. A couple of modifications to the project since the last post include (these details would only be relevant to those not terribly familiar with the Java tools used here):

Adding junit as a dependency to the pom.xml file.
Under module settings - libraries, look for he mahout-core:0.3 library and attach its source by going to the directory where you downloaded Mahout to (the same can be done with the Hadoop jar). This is convenient for being able to step into Mahout and Hadoop code from within you project while debugging within the IDE.

Other details include the input and output files, which can be deduced by reading the source for RecommenderJob (/Users/mario/Desktop/fastcf/data/recs/users can be a file as simple as having a single userID). We'll continue this next time.

Large-Scale Machine Learning with Mahout - Part 1

2010-07-04T18:36:00.000-07:00

Mining very large data sets is within everyone's reach nowadays, with open source libraries like Hadoop and Mahout.

Hadoop needs no introduction since most anyone interested in cloud computing has an idea of what it is. Even at such an early release number (0.2x) it is already heavily used in big name places like Facebook and Yahoo, just to mention a couple. For completeness sake, some of the seminal papers on which the Hadoop ecosystem is based include:

What is not as well-known is the Mahout project, which aims to build scalable machine learning libraries (and which leverages Hadoop but does not restrict itself to it).

Mahout was also inspired by a research paper:

Map-Reduce for Machine Learning on Multicore

That paper makes the observation that machine learning algorithms that fit the Statistical Query Model (those which compute sufficient statistics over the data) fit well into the map-reduce paradigm. For example, batch gradient descent, which computes the gradient contribution from each training instance and sums them all together to get the gradient, can be parallelized by splitting up the training data in several chunks, computing the local chunk contributions to the gradient in parallel, and then adding the chunk contributions for the overall gradient. An interesting result in that paper is that they achieved linear computational speed up with an increasing number of cores.

Even if machine learning is not particularly your cup-of-tea, examining the Mahout code can be a good exercise in learning how to leverage the power of Hadoop, and that is the point of this post and the next.

For starters, you'll wanna have access to the best guide to Hadoop currently out there, Tom White's book:

Hadoop: The Definitive Guide

In addition, there is a book about Mahout in the making:

Mahout in Action

And having that one within reach is also a good idea. This post and the next are somewhat gonna mirror part of chapter 6 from that book, "Distributing Recommendation Computations," but I'm aiming at going into a little more detail to get things setup, particularly for those out there who don't have a lot of experience with the Java ecosystem.

We'll be working with with tagged releases from each project, rather than the trunk. This is just so this discussion does not become obsolete in a couple of days given that the rapidly evolving nature of these projects.

I'm currently working on my Mac laptop, where I use RapidSVN to manage my subversion repositories. RapidSVN can be downloaded here. Once you have it, add bookmarks for the URL's of the existing Hadoop and Mahout repositories as shown in the screenshot below. Next, you will want to check out a working copy of the release tagged "Mahout-0.3" to a local directory of your choice.

After that, you will wanna do the same with Hadoop's release tagged "release-0.20.2" to its own directory.

Now that you have the sources locally, we are gonna open them up with an IDE. I'm currently using IntelliJ IDEA, which has a community edition version that can be downloaded for free here.

To open Mahout, create a new IDEA project and select "Import Project from external model" and the model will be a "Maven" model. Select the directory you downloaded the Mahout sources to as the project's root directory and with the remaining default settings untouched, click "Next" until the end to finish IDEA's project creation wizard.

At this point you'll wanna build the project. Mahout uses Maven, which is a "software project management and comprehension tool." In simple terms, Maven does in Mahout what Ant and Ivy do in Hadoop. Ant is the well-known tool for automating software build processes in the Java world, whereas Ivy is a dependency management tool. Maven, in turn, merges those capabilities (and more) into one tool. You'll wanna download Maven, for which you can follow these instructions, and then configure IDEA to use Maven. To configure IDEA, go to File - Settings - Maven and enter in the "Maven Home Directory" field the path to where you just installed Maven.

You'll wanna read up on the Maven lifecycle, but for now, just to make sure everything is alright, open the Maven Projects window in IDEA, and run the Clean followed by the Install tasks under the Apache Lucene Mahout node, as shown in the screenshot below.

This will download a bunch of dependencies, compile the code, and deploy it to you local Maven repository. My local repository happens to be under "/Users/mario/.m2" and you can find all the libraries needed by your Maven projects under there.

The whole thing may take a few minutes, but once it's done, you're all set. Building using IDEA's Build - Make Project menu command would be a pain (if you try it, you'll likely get several errors), since you would have to add a lot of dependencies, directories for code that gets generated during the build, etc. It's best to just let the Maven build handle it. Later on, we will be creating a project that references the Mahout library (and thereby all of Mahout's dependencies, including Hadoop), and so there you will be able to use the IDE. For now, let's open up Hadoop with IDEA.

Create again a new IDEA project, but this time select "create Java project from existing sources," and enter the directory where you downloaded the Hadoop sources to. IDEA will present you with all the source directories it has found; leave them all selected and click next. Next, merge all the libraries IDEA has found into one and call it "lib" and click next. Then, merge all the modules IDEA has as suggested into one and call it "Hadoop." Finally, open up the Ant Build tab in IDEA and run the "compile" target as shown below (other interesting targets include "compile-core-test" and "generate-test-records"). Again, building this project using the IDE can be a pain, so let Ant handle it. The Ant task includes Ivy tasks, which download dependencies just as Maven did. In Hadoop's case though, everything gets downloaded to within the Hadoop directory.

After the builds (in both Hadoop and Mahout), additional code is generated. If you wanna pickup the additional directories as modules, you can just recreate the project following the same steps. Alternatively, prior to opening the project in the IDEA, you can run the Maven and Ant tasks from the command line.

Once you have Hadoop built, you'll wanna install it locally and make sure you can run jobs in stand-alone and pseudo-distributed modes. Tom White's book is an excellent resource for setting this up, and you can just copy the contents of your Hadoop folder to another directory, say "hadoop-installation" and follow the instructions in Appendix A of the book to install hadoop from the hadoop-installation directory.

Once you are confident in your local Hadoop installation, create a new IDEA project, where you will copy some files from Mahout and examine one job in detail, the RecommenderJob, found in the org.apache.mahout.cf.taste.hadoop.item package.

The easiest way to do so is to create a new IDEA project from scratch and select Maven Module as the project type. Leave all the defaults alone and click next to finish the wizard. Next, you'll need to do 4 things:

Edit the pom.xml file to include a dependency on Mahout-0.3.
Add slf4j-jcl-1.5.8.jar as a library to your project (this library can be found in your local Maven repository). Add a new Library under Module Settings and attach the folder you copied the Jar to as a Jar Directory.
Copy all the files relevant to the RecommenderJob from Mahout to your project (the ones shown in the screenshot).
Under Project settings, set the project language level to "6 - @Override in interfaces"

Your project should look somewhat like the screenshot below. You should now be able to build using IDEA and using the Maven lifecycle tasks.

Next time we'll discuss RecommenderJob in detail.

Bayesian statistics

2010-06-09T18:24:00.000-07:00

I just finished a 2-course sequence on Bayesian statistics (ams206 and ams207) and doing so has certainly clarified many things I hadn't picked up on while trying to learn the topic by myself.

The Bayesian paradigm provides a very simple, elegant, and coherent foundation for building statistical models. I think of these models as being analogous to architectural blueprints for a building or UML diagrams for a software system.

In the same way that a software system may be implemented in Java or C++, using this or that framework, we may also implement our statistical model in a variety of ways. From the full Bayesian model, we may concede in "correctness" for the sake of computation complexity (simulation vs optimization), we may try this or that prior, or make this or that independence assumption; however, the full model remains explicit and our compromises are made clear.

For anybody looking for introductory resources on Bayesian stats as I once was, I would recommend Andrew Gelman's book, Bayesian Data Analysis, which includes answers to selected exercises (useful for self-study). Jim Albert's book, Bayesian Computation in R, will also come in handy from an application perspective, although I think the use of the LearnBayes package obscures things a bit on a first encounter (it would have been preferable to have seen the MCMC coded from scratch with lots of comments explaining the 'why' of things). Finally, once you're done with Albert's book, I would recommend Bayesian Modeling Using WinBUGS, and a good strategy would be to try to model and solve the examples in that book using R and then compare your answers with their WinBUGS implementation.

Bayesian Methods for Large Datasets

2010-05-12T18:43:00.000-07:00

Interesting talk on applying full bayesian models to large datasets (with the netflix challenge as the example application).

Here is the paper that goes with the talk.

A sign of the times

2010-05-09T07:36:00.000-07:00

Stanford has a new course offering in their statistics department this summer: Algorithmic Trading.

To my knowledge, this is one of the first major universities to have such a focused offering (are there any others?).

To put things in perspective, here are some numbers from Wikipedia:

"In 2006 at the London Stock Exchange over 40% of all orders were entered by algorithmic traders, with 60% predicted for 2007. American markets and equity markets generally have a higher proportion of algorithmic trades than other markets, and estimates for 2008 range as high as an 80% proportion in some markets. Foreign exchange markets also have active algorithmic trading (about 25% of orders in 2006). Futures and options markets are considered to be fairly easily integrated into algorithmic trading, with about 20% of options volume expected to be computer generated by 2010. Bond markets are moving toward more access to algorithmic traders."

Keep in mind that any insight into what goes on inside Wall Street is probably just the tip of the iceberg. For example, pairs trading, a technique covered in the course, is an old one:

"The pairs trade or pair trading, also known as market neutral, was developed in the late 1980s by quantitative analysts and pioneered by Gerald Bamberger while at Morgan Stanley. With the help of others at Morgan Stanley at the time, including Nunzio Tartaglia, Bamberger found that certain securities, often competitors in the same sector, were correlated in their day-to-day price movements. When the correlation broke down, i.e. one stock traded up while the other traded down, they would sell the outperforming stock and buy the underperforming one, betting that the "spread" between the two would eventually converge. A notable pair trader was hedge fund Long-Term Capital Management."

So what is this a sign of? As the Wikipedia article suggests, algorithmic trading is becoming more and more common place, and so employers are expecting future employees to know about it straight out of university, just as they need to know other essentials such as finance or economics.

I wonder if this trend will contribute to making "glitches" like the Dow's recent drop and subsequent rise of nearly 1,000 points in the space of an hour more common or less common.

Java vs C++

2010-05-04T16:19:00.000-07:00

Here is a completely unscientific yet entertaining piece of anecdotal evidence regarding java vs c++ performance:

So, I'm thinking about ways to scale recommender systems, and as a baseline, I'm using Timely Development's (TD) implementation of Simon Funk's well described approach to the Netflix Challenge.

One of the things I'm playing with is Mahout, which aims to be a scalable machine learning library, and being built on top of Hadoop, it is written in java. I rewrote TD's implementation in java, with the only change worth mentioning being that I changed the way the data is structured from something like this:


struct Data
{
int         CustId;
short       MovieId;
BYTE        Rating;
float       Cache;
};
Data            m_aRatings[MAX_RATINGS];

To something like this:


int[] m_aRatings_CustId;
short[] m_aRatings_MovieId;
byte[] m_aRatings_Rating;
float[] m_aRatings_Cache;

m_aRatings_CustId = new int[Constants.MAX_RATINGS];
m_aRatings_MovieId = new short[Constants.MAX_RATINGS];
m_aRatings_Rating = new byte[Constants.MAX_RATINGS];
m_aRatings_Cache = new float[Constants.MAX_RATINGS];

The reason for doing so being that java does not have structs and using a class instead of a struct would use more memory than I have in my laptop. Even with this change, the java implementation still uses more memory than the c++ implementation, but as long as I could run it on my laptop, I didn't care. Oh yeah, I also changed some Windows specific things in TD's code so I could run it on Linux, but that has a negligible effect in what the overall code is doing.

Running TD's code using 10 features and reading only 99 of the 17,770 movie files finishes fairly quickly:


real    0m49.097s

I was now curios to see how much slower the java code would run when given the same parameters:


real    0m26.520s

That brought a smile to my face :)

Handling uncertainty and complexity in AI

2010-04-09T16:14:00.000-07:00

We've been seeing several approaches to handling uncertainty (using probabilistic models) and complexity (using structured/relational representations) in tandem in artificial intelligence systems, including Church, the topic of an MIT news article released a couple of weeks ago. Other systems include:

This paper, from Michael Jordan, provides another approach towards an expressive probabilistic representation, one that draws from Bayesian nonparametrics.

Enhancing Video Resolution With Stills

2010-03-19T16:11:00.000-07:00

Enhancing videos or images is nothing new to Hollywood:

But some of the stuff that can be done back here in the real world can be quite impressive. I just finished an image processing class and got a taste of some of what's possible.

In one of the projects in the class, my friend and I tried to reproduce some of the results from a recent paper, Enhancing and Experiencing Spacetime Resolution with Videos and Stills, and though we did not implement all the components described in the paper (such as the spacetime fusion in the diagram below), we were still able to obtain some decent results.

Here's a video, the well-known foreman sequence, downloaded from here:

We processed the original sequence to simulate the output of a hybrid camera (as described in the original paper), which gave us the video below:

The flickering observed is caused by the intermediate high-resolution frames, the entire sequence being something like image below (red represents high-res frames, while gray represents low-res frames). In the hybrid foreman sequence, the gray frames are simply the original frames down-sampled by a factor of 4.

Finally, we processed the simulated hybrid input with our system, and we got this (if you watch closely, you'll observe several leftover artifacts, areas where the video appears to warp):

Even with the leftover artifacts, I'd say we were to propagate a decent amount of info from the high-res frames to the low-res frames, wouldn't you?

Data Modeling vs Algorithmic Modeling

2010-03-16T12:24:00.000-07:00

An interesting paper (plus comments) contrasting the 2 differing approaches used in statistics and machine learning.

Least-Squares For Image Restoration

2010-02-20T09:54:00.000-08:00

Here's an interesting use of least-squares: image restoration. You can find a good mathematical description of the matlab code in the paper titled Image Restoration by an Iterative Damped Least-Squares Method with Two Types of Adaptation ("damped" meaning regularized and "adaptive" is implemented in the code by using edge detection for spatially dependent regularization: no regularization around edges while regularizing on smooth sections of the image).


clc;clear all;close all;

filename='some_image.png';

%read image from disk
f=imread(filename);

%convert image to grayscale
f=rgb2gray(f);

%convert grayscale values to double
f=im2double(f);

%display original image
figure; imagesc(f); axis image; colormap(gray);

%compute the motion blur filter
h=fspecial('motion',72,315);

%blur the image
g=imfilter(f,h,'conv','circular');

%display the blurred image
figure; imagesc(g); axis image; colormap(gray);

%add gaussian noise to the image
mean=0;
sd=10/255;
noise=mean+sd*randn(size(f,1),size(f,2));
n=g+noise;

%display the noise
figure; imagesc(noise); axis image; colormap(gray);

%display the blurred and noisy image
figure; imagesc(n); axis image; colormap(gray);

%perform gradient descent to recover f from n
lambda=.01;
step=1;
error_norm=1;
laplacian=[0  1  0;
           1 -4  1;
           0  1  0];

%the noisy image is our first estimate
f_hat_0=n;
relative_error=intmax;
iter=0;
restf=figure;
wf=figure;

%repeat until little improvement is made
while relative_error>.1
    
    %compute left side of gradient equation
    f_inter_left=imfilter(f_hat_0,h,'conv','circular')-n;
    f_inter_left=imfilter(f_inter_left,h','conv','circular');
    
    %perform edge detection
    [w,t]=edge(f_hat_0,'canny',[.2 .3],2);
    %rescale edge intensities
    w=1./(1+.01.*w);
    
    %display the edges image
    figure(wf);
    imagesc(w); axis image; colormap(gray);
    drawnow;
    
    %compute right side of gradient equation
    f_inter_right=imfilter(f_hat_0,laplacian,'conv','symmetric');
    f_inter_right=f_inter_right.*w;
    f_inter_right=imfilter(f_inter_right,laplacian','conv','symmetric');
    f_inter_right=lambda*f_inter_right;
       
    %compute the gradient
    gradient=f_inter_left+f_inter_right;
    
    %update current image estimate
    f_hat_1=f_hat_0-step*(gradient);
    
    %compute relative error
    relative_error=(norm(f_hat_1-f_hat_0,error_norm)/...
        norm(f_hat_0,error_norm))*100;
    
    %reset our next estimate
    f_hat_0=f_hat_1;
    iter=iter+1;
    disp(['iteration: ',num2str(iter),...
        ', error: ',num2str(relative_error),...
        ', gradient: ',num2str(norm(gradient,2))]);
    
    %display our current image estimate
    figure(restf);
    imagesc(f_hat_0); axis image; colormap(gray);
    drawnow;
end

We start with the original image (found here):

And add some noise to it (in addition to blurring it), ending up with:

After one iteration of steepest descent, here is what edge detection looks like on the image:

And the image has not changed much:

After a couple hundred iterations of steepest descent, here's what the edges look like:

And the image is much clearer:

Markov Logic

2010-02-07T12:26:00.000-08:00

This is probably one of the most exciting things I've come across in terms re-unifying the several clusters AI has been divided into, and that aims to build a general-purpose, consistent, and realistic system: Markov Logic.

You can get a high level overview from this video, and find additional details here.

Boosting Bayesian Networks

2010-01-30T10:37:00.000-08:00

I strongly favor the idea of using a joint distribution as a knowledge base (memory?) for reasoning agents, but I don't necessarily think Bayesian networks in their current incarnation are the "end all, be all" representation. Specifically, I think the acyclicity requirement of the DAG may be too constraining (lately I've been particularly intrigued by dependency networks and their properties).

Nevertheless, the idea of boosting Bayesian networks for better density estimation described in this paper, Boosting Density Estimation, is very interesting.

Now, having an ensemble of Bayesian networks might be harder to interpret from a data mining perspective, but from an artificial intelligence perspective, I imagine a "language center" that knows how to query that knowledge base and transform the complex joint distribution representation into some sort of natural language form (and vice-versa).

Probabilistic Reasoning

2010-01-23T09:24:00.000-08:00

Assuming you're familiar with Bayes' rule (if not, check this out), let's review a classical problem, the Monty Hall problem:

"Suppose you're on a game show, and you're given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what's behind the doors, opens another door, say No. 2, which has a goat. He then says to you, "Do you want to pick door No. 3?" Is it to your advantage to switch your choice?"

The first time I saw that problem was in Text Mining Application Programming, and the solution went something like this:

Initially, the probability of the car being behind any of the doors is equal:

p(p_1)=p(p_2)=p(p_3)=1/3
Once the contestant picks door number 1 and the host opens door number 2, the contestant should use Bayes' rule to update his/her beliefs about where the prize is. The probability of the host opening door number 2 is 0.5, since he can only choose one of the two doors the contestant has not selected:

p(o_2)=1/2

Furthermore, we have the following conditionals:
1. The probability of the host opening door 2, given that the prize is behind door 1, is 0.5, since it makes no difference if the host chooses door 2 or door 3:
  
  p(o_2|p_1)=1/2
2. The probability of the host picking door 2, given that the prize is behind door 2, is 0:
  
  p(o_2|p_2)=0
3. The probability of the host picking door 2, given that the prize is behind door 3, is 1:
  
  p(o_2|p_3)=1
So, applying Bayes' rule with the numbers we have, we get the following:

p(p_1|o_2)=p(o_2|p_1)*p(p_1)/p(o_2)=1/3
p(p_2|o_2)=p(o_2|p_2)*p(p_2)/p(o_2)=0
p(p_3|o_2)=p(o_2|p_3)*p(p_3)/p(o_2)=2/3

According to this, the contestant can double his/her chances of getting the car by switching from door 1 to door 3. Not an intuitive result, as far as I'm concerned, but nevertheless, the numbers seem to make sense.

Perhaps you're not convinced you should switch. Perhaps you think the problem is not well defined enough and leaves much to interpretation. In that case, here is a more detailed enunciation of the problem:

"Suppose you’re on a game show and you’re given the choice of three doors. Behind one door is a car; behind the others, goats. The car and the goats were placed randomly behind the doors before the show. The rules of the game show are as follows: After you have chosen a door, the door remains closed for the time being. The game show host, Monty Hall, who knows what is behind the doors, now has to open one of the two remaining doors, and the door he opens must have a goat behind it. If both remaining doors have goats behind them, he chooses one randomly. After Monty Hall opens a door with a goat, he will ask you to decide whether you want to stay with your first choice or to switch to the last remaining door. Imagine that you chose Door 1 and the host opens Door 2, which has a goat. He then asks you “Do you want to switch to Door Number 3?” Is it to your advantage to change your choice?"

Surprisingly, the analysis above remains and the advantage in switching also remains :)

Confused, still? Well, let's briefly look at another classical problem, the Three Prisoners problem:

"Three prisoners, 1, 2, and 3, have been tried for murder. Exactly one will be hanged tomorrow morning, but only the guard knows who. 1 asks the guard to give a letter to another prisoner — one who will be released. Later 1 asks the guard to whom he gave the letter. The guard answers '2'. 1 thinks, '2 will be released. Only 3 and I remain. My chances of dying have risen from 1/3 to 1/2.' "

Is prisoner 1 correct in his thinking?

I first came across this problem in Probabilistic Reasoning in Intelligent Systems, where the author explains that no, prisoner 1 is not correct in his thinking: his probability of dying remains 1/3.

Let's apply Bayes' rule to this problem:

p(g_1|i_2)=p(i_2|g_1)*p(g_1)/p(i_2)=?

The initial probability that 1 is guilty, p(g_1), is 1/3. The probability that 2 is innocent, p(i_2), is 2/3. Finally, the probability that 2 is innocent given that 1 is guilty, p(i_2|g_1), is 1. So, p(g_1|i_2)=1/2.

Wait... huh? I thought prisoner 1's probability of dying was to remain 1/3... what happened? Also, the Prisoners problem sounds so similar to the Monty Hall problem... assuming p(g_1|i_2)=1/3 is the correct answer, why is it that our beliefs should be updated given the new information in the Monty Hall problem, but not in the Prisoners problem?

It is at this point that I am reminded of the mathematical proof that 1=2:

Let a=b

Then a^2=ab

a^2+a^2=a^2+ab

2a^2=a^2+ab

2a^2-2ab=a^2+ab-2ab

2a^2-2ab=a^2-ab

2(a^2-ab)=1(a^2-ab)

Canceling the (a^2-ab) from both sides gives 1=2

The fallacy is in the last step, where we cancel (a^2-ab) from both sides. We cannot do that because (a^2-ab)=0 and we cannot divide by zero (see here). The moral of the story is that we can't blindly apply a method without thinking about what we're doing, a point that is sometimes obvious (as in 1=2) and sometimes much more subtle (as in the Prisoners problem).

The fallacy in our analysis of the Prisoners problem consists in not conditioning on the right piece of information. The probability that 2 is innocent, 2/3, is not what we should be conditioning on since that is not exactly the new information that we acquired. The information we acquired was that the guard said that 2 would be declared innocent. The guard could only provide one of two answers to prisoner 1:

2 will be released, or
3 will be released

So, the probability that the guard said 2 will be released, p(r_2), is 1/2, and that is what we should be conditioning on:

p(g_1|r_2)=p(r_2|g_1)*p(g_1)/p(r_2)=?

In this case, the probability that the guard would have said 2 was gonna be released given that 1 was guilty, p(r_2|g_1), is also 1/2 since he could have as well chosen 3 to be whom he'd give the letter to, so p(g_1|r_2)=1/3.

As we can see, our beliefs are updated by the Bayesian mechanism in both problems, but whereas in the Monty Hall problem the new information is relevant to our beliefs, in the Prisoners problem it isn't.

If you're still not convinced about the Monty Hall problem, check the page Wikipedia has for it, a surprisingly long page. There you will find a suggestion to consider a variation of the problem where rather than 3, there are 1,000,000 doors instead:

"In this case there are 999,999 doors with goats behind them and one door with a prize. The player picks a door. The game host then opens 999,998 of the other doors revealing 999,998 goats—imagine the host starting with the first door and going down a line of 1,000,000 doors, opening each one, skipping over only the player's door and one other door. The host then offers the player the chance to switch to the only other unopened door. On average, in 999,999 out of 1,000,000 times the other door will contain the prize, as 999,999 out of 1,000,000 times the player first picked a door with a goat. A rational player should switch. Intuitively speaking, the player should ask how likely is it, that given a million doors, he or she managed to pick the right one."

Boosting

2009-12-11T13:20:00.000-08:00

A short paper I wrote for a class, exploring boosting for classification, part of a bigger project I'm working on.

Also, here is some code for boosting MATLAB's classregtree:


function [ensemble,iter,w]=SAMME_trees(X,y,w,T,ensemble,t)
%m is the number of training instances.
m=size(X,1);
%number of classes.
k=10;

for iter=t:T-1 
    %train classifier on weighted dataset and compute its error.
    %sample according to w.
    randnum=rand(1,m);
    cW=cumsum(w);
    indices=zeros(1,m);
    for i=1:m,
        %find which bin the random number falls into
        idx=find(randnum(i)>cW, 1,'last')+1;
        if isempty(idx)
            indices(i)=1;
        else
            indices(i)=idx;
        end
    end
    
    [tree,error,diff_vec]=train(X(indices,:),y(indices,1),X,y);
    if error<=0 || error>=(1-(1/k))
        iter=iter-1;
        return
    end
    
    %calculate the weight of the classifier based on its error.
    alpha=log((1-error)/error)+log(k-1);
    
    %add the classifier and its weight to our ensemble.
    ensemble{iter,1}=alpha;
    ensemble{iter,2}=tree;
    
    %update the datapoint weights.
    w=w.*exp(alpha.*(diff_vec.*2-1));
    w=w./sum(w);
    
    disp(['boosting iteration: ',num2str(iter)]);     
end
%nothing here!
end

function [tree,error,diff_vec]=train(X,y,X_to_pred,y_to_pred)
%train classifier.
tree=classregtree(X,y,'method','classification','prune','off');
%predict everything.
y_hat_cell=tree(X_to_pred);
%convert everything to int vector.
y_hat=int32(str2num(char(y_hat_cell)));
%see what was missed.
diff_vec=ones(size(y));
diff_vec(y_to_pred==y_hat)=0;
error=sum(diff_vec)/size(diff_vec,1);

disp(['missclassification ratio: ',num2str(error)]);
end

Driven with:


clc;clear;
stream0=RandStream('mt19937ar','Seed',0);
RandStream.setDefaultStream(stream0);

out_iters=1;
in_iters=10;
accuracies=zeros(out_iters,in_iters);
accuracies_tolerance_one=zeros(out_iters,in_iters);
    
%read data.
source_data='../data/winequality-white.csv';
covariate_num=12;
data=textread(source_data,'','delimiter',';','headerlines',1);
header=textread(source_data,'%s',covariate_num,'delimiter',';');
k=10;
%m is the full number of data instances (rows).
m=size(data,1);
%n is the number of attributes (take one out, cause of target).
n=size(data,2)-1;
%separate instances into attributes and targets.
X=data(:,1:end-1);
y=int32(data(:,end));

for out_iter=1:out_iters
        
    %split data into train and test sets.
    probs=rand(1,m);
    prob_thresh=2/3;
    X_train=X(probs<=prob_thresh,:);
    y_train=y(probs<=prob_thresh,:);
    X_test=X(probs>prob_thresh,:);
    y_test=y(probs>prob_thresh,:);
    
    %w is the weight assigned to each data point.
    m_train=size(X_train,1);
    w=ones(m_train,1)./m_train;
    %will be storing rows of [alpha,tree].
    boost_iters=2;
    ensemble=cell(in_iters*boost_iters,2);
    t=0;
    
    for in_iter=1:in_iters
        %fit model.
        t=t+1;
        [ensemble,t,w]=SAMME_trees(X_train,y_train,w,boost_iters+t,...
            ensemble,t);
        
        %predict targets on test set.
        m_y_test=size(y_test,1);
        candidates=zeros(m_y_test,k);
        for i=1:t
            alpha=ensemble{i,1};
            tree=ensemble{i,2};
            y_pred_i=int8(str2num(char(tree(X_test))));
            for j=1:m_y_test
                candidates(j,y_pred_i(j,1))=...
                    candidates(j,y_pred_i(j,1))+alpha;
            end
        end
        [val,y_pred]=max(candidates,[],2);
        y_pred=int32(y_pred);
        
        %calculate confusion matrix and accuracies.
        [confus,legend]=confusionmat(y_test,y_pred)
        total=sum(sum(confus))
        accuracy=sum(diag(confus))*100.0/total
        accuracy_tolerance_one=sum(sum(spdiags(confus,[-1 0 1])))*100.0/total;
        
        accuracies(out_iter,in_iter)=accuracy;
        accuracies_tolerance_one(out_iter,in_iter)=accuracy_tolerance_one;
    end
end
accuracies=accuracies';
accuracies_tolerance_one=accuracies_tolerance_one';

Best Data Mining Algorithm

2009-12-04T07:54:00.000-08:00

Which data mining algorithm is the best one obviously depends on whatever you may take 'best' to mean, and picking a single one may be a futile exercise in the current state of affairs, but this paper (published in 2007) does a good job discussing the top 10 data mining algorithms identified by the IEEE International Conference on Data Mining (ICDM) in December 2006:

C4.5
k-Means
SVM
Apriori
EM
PageRank
AdaBoost
kNN
Naive Bayes
CART

Easy Learning of Bayesian Networks

2009-11-25T08:49:00.000-08:00

The other day I came across a very interesting idea, described in this paper, where the authors discuss learning a Bayesian network by first learning a dependency network.

Why would you want to do that? Well, a dependency network can be thought of as a Bayesian network where cycles are allowed in the graph, and so if you lift the acyclicity restriction, it turns out to be really easy (almost trivial) to learn a probabilistic graphical model: all you have to do is learn a conditional probability distribution (CPD) for each variable given the other variables. This is a supervised learning problem and we know many ways of solving it (any probabilistic discriminative method, particularly if coupled with a nice feature selection method, will do just fine). Furthermore, you can learn each node's CPD in parallel.

Then, to obtain a Bayesian network from the dependency network, all you have to do is remove the least significant edges participating in cycles in the dependency network. The paper describes a very easy way of doing so. However, I'm starting to wonder if I'd even wanna do that. Why not just stick with the dependency network? If causality matters to you, you may argue that causality is not cyclical, or at best simply hard to represent cyclically, but I think that ceases to be an issue when you consider the time dimension (more about this on a future post).

Causality aside, as you may have expected, a dependency network is not a panacea. A joint distribution doesn't factorize as nicely over a dependency network as it does over a Bayesian network. For example, if we had 2 variables, X and Y, with a mutual cyclic dependency (i.e. a bidirectional edge between the 2 nodes), we'd hope to be able to write something like: P(X,Y)=P(X|Y)P(Y|X), but if you simplify that out using the chain rule P(X,Y)=P(X|Y)P(Y), you'd end up with P(X,Y)=P(X)P(Y), which says the 2 are independent. Also, it is possible, and it does happen often (particularly if learning a single tree for each CPD, something that could be alleviated with surrogate splits), that you may end up with X being a significant variable in Y's conditional distribution, but Y not being significant in X's conditional distribution. No big deal, though, you can still recover the joint distribution and answer any probabilistic query using Gibbs Sampling, and that might just be good enough.

For illustration purposes, I used R's rpart package to learn a dependency network structure for this dataset (using only 10 of the 14 attributes), and this is what it came up with (rendered with Graphviz):

Reactive vs Deliberative Agents

2009-10-03T19:01:00.000-07:00

So, what exactly are agents? What does it mean for an agent to be reactive or deliberative and what are [dis]advantages of one over the other?

From the Intelligent Agent Book, an agent is anything that can be viewed as perceiving its environment through sensors and acting upon that environment through actuators. I would add (as does the author) that a most important component of the agent is its performance measure, the details of which are not under control of the agent, but instead are provided to it by its designer/creator (not unlike the reward function in reinforcement learning). It is the purpose of the agent, then, to act in its environment so as to maximize its performance measure.

It might be interesting to contrast this "designer/creator" approach of reinforcement learning to an approach based purely on evolutionary techniques such as genetic algorithms. Both approaches are in fact characterized by a performance function which must be maximized but a differentiating characteristic is that reinforcement learning agents learn from [trial-and-error] interaction with the environment whereas evolutionary techniques traditionally don't. Are these approaches mutually exclusive or are they instead complementary? Could reproduction be the link which establishes the connection between reinforcement learning and evolutionary techniques?

There are other properties of agents which may help in organizing them into some sort of taxonomy. One such property is the degree of autonomy a given agent possesses. Autonomy can be a measure of how much an agent relies of the prior knowledge of its designer rather than on its own precepts to act in its environment. The more an agent relies on its own precepts, the more autonomous it is. Therefore, simplistically, a reinforcement learning agent will tend to be more autonomous than a candidate solution (an agent, if you will) in a genetic algorithm.

Some examples? OK, well, the most obvious example of an agent is a human being, which has ears, eyes, etc, for sensors, and hands, legs, etc, for actuators. Another example would be a robot, which could have cameras and radars for sensors and wheels and robotic arms for actuators:

A more abstract example of an agent would be a stock trading software program, designed to "sense" a stream of stock price quotes and provided with the ability to buy or sell those stocks (its "actuators"), a simple instance of which can be found here.

All of these agents exhibit a certain degree of autonomy (the stock trading agent being the least autonomous; go look at the source code), and can be thought of as entities which sense and act in their environments, and as they do so, they build a model of what happens when they take a certain action while the environment is in a given state. This model, or abstract representation, of their environment could possibly help the agent "reason" about what action to take next so as to maximize its performance measure not only today, but throughout the lifetime of the agent. The degree to which a model evolves with the agent's experience is an indicator of the agent's level of autonomy (more plasticity indicates more autonomy).

A related and very interesting perspective is that described in The Hedonistic Neuron: A Theory of Memory, Learning, and Intelligence, by A. Harry Klopf. That name may sound familiar to you, especially if you've read Sutton and Barto's introduction to reinforcement learning. Harry Klopf is whom the book was dedicated to, and in the last section of chapter I, History of Reinforcement Learning, there is a mention of how influential Harry was on shaping the authors' ideas on reinforcement learning. To give you a sense of what "The Hedonistic Neuron" is about here is how the preface begins:

"At the level of whole organisms, it is commonly observed that animals pursue pleasure and avoid pain. At the level of the single neuron, it is well known that the individual brain cell receives two classes of synaptic inputs: excitation and inhibition. What is the relationship between these two descriptions, one mentalistic and at an organismic level, the other physicalistic and at a neuronal level? Is there no simple relationship, as most neuroscientists would be likely to assume today? Have we described phenomena at two such vastly different levels and in such vastly different terms that the complexity and subtlety of their relationship will defy elucidation for some time to come? Or, to go to the other extreme, might their relationship be the simplest one conceivable? Namely, might excitation and inhibition represent, in elementary physical terms, one and the same thing that pleasure and pain represent in complex mental terms? But, at this point, can we even say precisely what we would mean by such questions? And even if we can, are not such questions simplistic to the point of being unproductive? This book seeks to demonstrate that we can say what we mean and the we will benefit from asking such questions."

Now, the book was published in 1982, and much has changed in neuroscience since then, nevertheless, much still remains unclear about what it is that neurons do in cognitive terms (even if we do know a little more about how it is that they do it in chemical terms). So, don't expect to read about the latest and greatest explanations regarding the physiological and cognitive properties of the brain if you decide to read Harry's book, but do expect to come out with some very thought-provoking ideas. In fact, evidence of the book's transcendence is that today we have a hedonistic synapse hypothesis: "the assumptions are that operant conditioning in living organisms has an underlying neuro-chemical mechanism and that this mechanism could be implemented by means of each synapse optimizing its behavior to harvest the most of globally broadcasted reward" (check out this paper from Sebastian Seung). So, it may be synapses and not neurons that exhibit a goal-seeking behavior, but we may be splitting hairs here. The idea is that we could have a goal-seeking system composed of goal-seeking elements:

"The theory to be developed explains goal-seeking brain function in terms of goal-seeking neurons. It has been generally (and implicitly) assumed in the past that advanced (intelligent) goal-seeking brain function emerges from the interaction of non-goal-seeking neurons. Assuming such a passive, non-goal-seeking role for the neuron may not be valid. The single neuron is a remarkably complex and sophisticated cell and it may well play a more active role. Perhaps an analogy will help make the point clear. Consider that goal-seeking social systems (such as the United States) would probably remain mysterious if we assumed that the people making up the social system were non-goal-seeking in nature. There would probably be no way of explaining complex goal-seeking social system behavior (such as putting a person on the moon) in terms of the interactions of non-goal-seeking people. Are we any more likely to be successful in understanding goal-seeking nervous systems by assuming non-goal-seeking neurons?"

Perhaps embedded in the mechanism by which a neuron makes a binary decision (to fire or not to fire) is the essence of free-will? Could a neuron be described as an agent that senses a stream of inhibitory and excitatory synaptic inputs (possibly encoding the temporal and spatial patterns in the stream into some sort of model) and then deciding whether or not to fire (whether or not to act)? Could clusters of neurons conspire so as to better satisfy their individual hedonistic needs? What role would knowledge play in a world in which these ideas weren't so far-fetched?

An answer to that last question could be that knowledge is that which enables an agent to predict the evolution of an given process through time, in scenarios where the agent is simply observing the process or where the agent is attempting to influence the process; think Markov chain vs Markov decision process, respectively (although the Markov assumption, that knowledge of the present renders the past and the future independent, though useful, may be overly simplistic and I do not imply here that it is crucial for more general modeling purposes):

Markov Models		Do we have control over the state transitions?
Markov Models		No	Yes
Are the states completely observable?	Yes	Markov Chain	Markov Decision Process
Are the states completely observable?	No	Hidden Markov Model	Partially Observable Markov Decision Process

So, as an agent perceives and acts in its environment, it builds a model of how it sees the environment evolving, and uses that model to make decisions about what actions to take next so as to maximize its performance measure. The model, then, represents the agent's knowledge, and it is how the agent uses its model (how the agent "reasons") that differentiates reactive from deliberative agents.

The simplest kind of a reactive agent is one which possesses a direct mapping between states and actions, so that at each point in time, no matter what state the environment is in, there would be a function which given the current state, would output the best possible action to be taken by the agent, presumably one that would maximize the agent's expected reward from the environment. For a concrete example of this kind of agent, take a look at lecture 20 of CS229, where policy search algorithms are discussed (the kind of reinforcement learning algorithms that have been successfully applied to fly autonomous helicopters). Briefly, a policy is a function that maps states to actions, and the algorithm searches for a policy that maximizes the agent's expected return (very much like how gradient descent algorithms can be used to search for a set of parameters which maximize a function's fit to certain data).

So, an agent whose controller/brain was trained using a policy search algorithm would know exactly what to do on environment states it was specifically trained/simulated on, but it would be severely limited on anything beyond that. For complex tasks, training a controller sure beats hand-coding it, but for even more complex tasks, policy search may require extensive, possibly impractical, training. This training usually happens in a simulator, and the line between reactive and deliberative agents can be blurry (especially considering the generalization power of learning continuous functions), but for illustration purposes, imagine a simple reactive agent that is faced with an environment state it has never been trained on and therefore has no state-action mapping for this particular state: the agent would be at a loss of what to do. Granted, we could have some sort of default action for unfamiliar states, but I hope you get a feel for how restrictive this approach can be, or conversely, how much work you have to do in a simulator to train a good controller of this kind. On the other hand, given a state, policy search algorithms require minimal computation in order to select actions, so they can have very fast "reaction times."

Now what if we equip our agent with some more tools so it isn't as stuck if it finds itself in an unfamiliar situation? We could design our agent so that it has a model of the environment that allows it to "predict" how the environment will transition from state to state given a certain action by the agent. In addition, this agent could learn a state-value function based on rewards it receives from the environment when in a given state (that way it knows what states it prefers and hopefully will seek to be in those states). Voila! Now the agent can ponder/deliberate on the results of its actions and presumably select a plan (a sequence of actions), not just a single action. Algorithms that do this fall under the umbrella of value-function approximation (see lecture 17 of CS229), and evidently, agents driven by this kind of algorithms appear to be much more flexible than agents based on policy search. This deliberation, however, requires computational resources, in other words, it takes time to simulate several levels of state transitions resulting from action sequences. We have several options in terms of how to model the state-transition-function and value-function, giving rise to a broad spectrum of agents with probabilistic planning capabilities (see Sebastian Thrun's Probabilisitc Robotics for more details).

Can we do more to enhance our agents' deliberative capabilities? Actually, there are some very neat ideas in Causality that suggest we can. The reading can be somewhat dense, but the punchline is that the kinds of models and algorithms therein presented may allow an agent to infer the consequences of actions the agent might not have known a priori that it had the capability to execute. So, rather than asking itself "how would things turn out if I took this action?", an agent could ask itself "how would things turn out if I could take this action?", and if the outcome was favorable, the agent could work towards being able to take the hypothesized action.

The kind of knowledge that would allow an agent to ask and answer that hypothetical question is much more complex than the knowledge with which we have traditionally endowed classical reinforcement learning agents (state transition probabilities and possibly a richly structured graphical representation of a given time slice). It is the kind of knowledge that would require causal modeling and causal reasoning, the kind of knowledge which would allow an agent to answer interventional queries from non-experimental data and to answer conterfactual queries; basically, the kind of knowledge that would allow an agent to run complex thought experiments. That's some pretty powerful stuff!

Singularity University

2009-09-23T17:10:00.000-07:00

Cool video from the Singularity University:

Bayesian Network Approximate Inference

2009-09-12T13:18:00.000-07:00

Sampling methods for bayesian network inference, also known as particle-based methods (each sample being a particle), are just like any other Monte Carlo method: the more samples, the better the approximation (approximate inference, as exact inference, is NP-hard, that is, if you want arbitrarily accurate answers). A sample consists of a particular instantiation of the bayesian network (an assignment of values to each of the nodes in the network) and the approximation sought is that which answers a probabilistic query, which in its most basic form would look something like p(x|e)=?, where 'x' are unobserved nodes, 'e' are evidence/observed nodes, and '?' would be the posterior distribution of 'x' given 'e'.

Though there are many other kinds of algorithms for performing inference in probabilistic graphical models (PGM), I like sampling methods because they are very intuitive, they highlight the generative nature of PGM's, and finally, because they are more generally applicable than the other methods (even if often a bit challenging to get to work properly).

For illustration purposes, we will use Likelihood Weighting (LW), which is a special case of a general sampling framework called Importance Sampling, of which you can learn more from this video lecture.

In LW, to generate a particle, we sample the nodes in the network in an orderly fashion (ensuring parent nodes are sampled before their children), using the nodes' conditional probability distributions (CPD's), and setting evidence nodes to their observed values. Finally, we assign to the resulting particle a weight corresponding to the likelihood of seeing that particular particle. Having a bunch of weighted particles, we can then generate a conditional distribution for the query variables. More concretely (using BNT and the same network as last time):


clear;
clc;

%number of nodes
N=4; 
%graph adjacency matrix
dag=zeros(N,N);
%nodes order topologically
C=1; 
R=2;
S=3;
W=4;
%set grph edges
dag(C,[R S])=1;
dag(R,W)=1;
dag(S,W)=1;

false=1;
true=2;

%node sizes vector (all binary)
ns=2*ones(1,N); 
%crate the bayesian network
bnet = mk_bnet(dag,ns);
%create the nodes' CPD's
bnet.CPD{C}=tabular_CPD(bnet,C,[0.5 0.5]);
bnet.CPD{R}=tabular_CPD(bnet,R,[0.8 0.2 0.2 0.8]);
bnet.CPD{S}=tabular_CPD(bnet,S,[0.5 0.9 0.5 0.1]);
bnet.CPD{W}=tabular_CPD(bnet,W,[1 0.1 0.1 0.01 0 0.9 0.9 0.99]);

%create the LW inference engine
engine=likelihood_weighting_inf_engine(bnet);

%assume no observed nodes
evidence=cell(1,N);
samples=100;
[engine,loglik]=enter_evidence(engine,evidence,samples);
%what's the prior probability for 'wet grass'?
distr=marginal_nodes(engine,W);
distr.T

%now assume sprinkler is off, what's the
%conditional probability for 'wet grass'?
evidence{S}=false;
[engine,loglik]=enter_evidence(engine,evidence,samples);
distr=marginal_nodes(engine,W);
distr.T

%now assume it is raining, what's the
%conditional probability for 'wet grass' now?
evidence{R}=true;
[engine,loglik]=enter_evidence(engine,evidence,samples);
distr=marginal_nodes(engine,W);
distr.T

Which yields:


%prior for W implies grass is slightly more 
%likely to be wet than not
ans =

    0.3700
    0.6300

%once we learn the sprinkler is off,
%the grass is almost as likely to be wet
%as it is to be dry
ans =

    0.4778
    0.5222

%once we learn it is raining though, the
%grass becomes much more likely be wet
ans =

    0.1433
    0.8567

The magic happens when we make the call to 'enter_evidence', which looks like this:


function [engine, ll] = enter_evidence(engine, evidence, nsamples)
% ENTER_EVIDENCE Add the specified evidence to the network (likelihood_weighting)
% [engine, ll] = enter_evidence(engine, evidence, nsamples)
% evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value 
%(scalar or column vector)
%
% If nsamples is not specified, the value specified when the engine was created will be used.
% ll (log-likelihood) is set to [].

ll = [];
if nargin < 3, nsamples = engine.nsamples; end

bnet = bnet_from_engine(engine);
N = length(bnet.dag);
samples = cell(nsamples, N);
weights = zeros(1, nsamples);

ns = bnet.node_sizes;
original_evidence = evidence;
observed = ~isemptycell(original_evidence);
for s=1:nsamples
  evidence = original_evidence(:); % must be a column vector
  w = 1;
  for i=1:N
    ps = parents(bnet.dag, i);
    e = bnet.equiv_class(i);
    if observed(i)
      p = exp(log_prob_node(bnet.CPD{e}, evidence(i), evidence(ps)));
      w = w * p;
    else
      x = sample_node(bnet.CPD{e}, evidence(ps));
      evidence{i} = x;
    end
  end
  samples(s,:) = evidence;
  weights(s) = w;
end                 

engine.samples = samples;
engine.weights = weights;

log_prob_node uses a node's CPD to compute the [log] probability of the observed value for a given observed node conditioned on the values of its parents and this value is then used to calculate the weight for a given particle. sample_node samples values from unobserved nodes given the values of their parents. The particles are then aggregated and normalized in the call to marginal_nodes. Below is a fraction of the particles and their weights generated in this example:


%sprinkler is off and it is raining
%nodes are in the specified topological order
[1]    [2]    [1]    [2]
[2]    [2]    [1]    [2]
[2]    [2]    [1]    [2]

%weights associated with each of those particles
0.1000    
0.7200    
0.7200

Probabilistic Graphical Models

2009-09-09T10:18:00.001-07:00

So, I recently got a copy of Daphne Koller's Probabilistic Graphical Models (PGM), and it is massive: 1200 pages of pretty much everything you've ever wanted to know about bayesian [and markov] networks.

To be honest, I would have liked to have seen some stuff on modeling Markov Decision Processes with Dynamic Bayesian Networks, but the book already covers so much material that I'm willing to overlook that detail (she does mention some good pointers on DBN's for MDP's at the end of chapter 23).

There may be gentler introductions into the field out there (Learning Bayesian Networks, for example), but if you want to go beyond the introduction and build a good, solid foundation, PGM is the way to go.

The content is structured within 4 main sections:

Representation
Inference
Learning
Actions and Decisions

And you can take a look at the rest in the table of contents, but those 4 sections are a very descriptive way of summarizing what PGM's are all about. Briefly, we assume that a knowledge base can be represented as a joint probability distribution. This is by no means a far-fetched assumption: even if you've never taken a probability class, you use the concepts in your daily discourse. Consider the following examples, from Pearl's Probabilistic Reasoning in Intelligent Systems (PRIS):

Likelihood: "Tim is more likely to fly than to walk"
Conditioning: "If Tim is sick, he can't fly"
Relevance: "Whether Tim flies depends on whether he is sick"
Causation: "Being sick caused Tim's inability to fly"

Believe it or not, when you use language like that, you are in fact making probabilistic assertions :)

Pearl wrote PRIS at a time when most people were trying to build knowledge systems based on logic, but that approach turned out to be riddled with problems, a big one being the logic systems' inability to handle uncertainty. So, building on the work of those who came before him, he successfully made the case for and launched us into this new era of structured probabilistic reasoning systems.

The reason for adding the qualifier "structured" and not simply "probabilistic reasoning systems," is that naively manipulating a joint probability distribution can be a very computationally expensive proposition. So, we exploit the structure inherent in the problem domain. More specifically, we exploit conditional independences in our joint probability distribution to drastically reduce the number of computations we need to perform (basically, as you head out the door for a short evening walk, you don't need to know what the traffic is like in Japan in order to decide whether or not to take an umbrella with you in New York). It turns out graphs (networks) are a very natural data structure with which to exploit these conditional independences and write computer algorithms for.

PRIS is a very interesting read, even though it may be slightly dated by now. Highlights include several discussions on the merits of probabilistic systems over logic systems, for example, a discussion on Polya, who "argued that the process of discovery, even in as formal a field as mathematics, is guided by nondeductive inference mechanisms, entailing a lot of guesswork" (as opposed to a system that would have all the basic facts and would reach all other conclusions through deductive inference). The discussion is about Polya's so-called "patterns of plausible inference," which was his term for the patterns he identified as representing these nondeductive inference rules that could potentially serve as reasoning mechanisms in some reasoning system. Interestingly, after he had identified a bunch of these patterns and attempted to formalize them into a coherent set of reasoning rules, he [wisely] basically said "screw it," and settled upon the calculus of probabilities as a meta-pattern from which all his patterns could be derived.

So, we have PGM's as an efficient way to represent joint probability distributions, and for a mental image, you can just think of a graph whose nodes are strategically linked to each other so as to maximally exploit conditional independences; the nodes representing the random variables of interest in our domain. Each node has attached to it a conditional probability distribution (CPD), as is the case in a bayesian network, or more generally, a factor (as in the case of a markov network), and we say that the joint probability distribution factors into the CPD's and over the graph. Pictorially, this:


joint(:,:,1,1) =

    0.2000    0.0200
    0.0900    0.0010

joint(:,:,2,1) =

    0.0050    0.0005
    0.0360    0.0004

joint(:,:,1,2) =

         0    0.1800
         0    0.0090

joint(:,:,2,2) =

    0.0450    0.0495
    0.3240    0.0396

Splits into each of the CPD's (depicted below as tables next to each node), and each CPD goes to its node as specified by the structure of the graph:

The joint distribution can be easily generated with this script, using BNT:


N = 4; 
dag = zeros(N,N);
C = 1; S = 2; R = 3; W = 4;
dag(C,[R S]) = 1;
dag(R,W) = 1;
dag(S,W)=1;

false = 1; true = 2;
ns = 2*ones(1,N); % binary nodes

%bnet = mk_bnet(dag, ns);
bnet = mk_bnet(dag, ns, 'names', {'cloudy','S','R','W'}, 'discrete', 1:4);
names = bnet.names;
%C = names{'cloudy'};
bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);


CPD{C} = reshape([0.5 0.5], 2, 1);
CPD{R} = reshape([0.8 0.2 0.2 0.8], 2, 2);
CPD{S} = reshape([0.5 0.9 0.5 0.1], 2, 2);
CPD{W} = reshape([1 0.1 0.1 0.01 0 0.9 0.9 0.99], 2, 2, 2);
joint = zeros(2,2,2,2);
for c=1:2
  for r=1:2
    for s=1:2
      for w=1:2
 joint(c,s,r,w) = CPD{C}(c) * CPD{S}(c,s) * CPD{R}(c,r) * ...
     CPD{W}(s,r,w);
      end
    end
  end
end

joint

The 'representation' section of the book goes into all the gory details, not just on bayesian networks, but markov networks as well, and deeper into how to represent the CPD's, template-based representations, etc, etc.

OK, so we got this great way to represent knowledge, now how do we use it? Well, we use it to answer probabilistic queries, and the algorithms for doing so are the topic of the 'inference' section of the book, where we infer/reason with our knowledge base.

Now, it can be a pain to create this models by hand, which would entail:

Create the graph structure being careful to maximize our ability to exploit conditional independences
Parameterize the graph with CPD's

Fortunately, there exist algorithms to ease that burden, and those are discussed in the 'learning' section, where the goal is to automatically learn/infer the structure and parameterization of the network from raw data.

Finally, the last section, 'actions and decisions', brings it all together and provides insights into issues one might have to deal with in order to build an intelligent agent based on these principles. Topics discussed include causality, utilities and decisions, and structured decision problems.

Machine Learning and Computational Finance

2009-08-24T15:04:00.000-07:00

Videos for the recent Advances in Machine Learning for Computational Finance conference are posted on the videolectures.net site:

Advances in Machine Learning for Computational Finance

Good stuff!

Trading Strategy Implementation

2009-08-20T18:12:00.000-07:00

Disclaimer: run this strategy at your own peril; you will most likely lose all your money. This is purely for illustrative purposes.

OK, now that I've got that disclaimer out of the way, we can continue. Seriously though, it's so easy to get connected and start trading that I am sure many people have gotten caught up in the excitement and whatever little knowledge they may have acquired from some infomercial-like seminar or book, that they went in head-first and ended up with a story to tell but no profits to show. Make sure you understand what you're getting into very well if you decide to experiment with trading (let alone automated trading), although you are probably better off just reading Malkiel's Random Walk Guide To Investing and following the 10 rules proposed there. Most of that is nicely summarized in Dilbert's investment advice:

Make a will.
Pay off your credit cards.
Get term life insurance if you have a family to support.
Fund your 401(k) to the maximum.
Fund your IRA to the maximum.
Buy a house if you want to live in a house and you can afford it.
Put six months’ expenses in a money market fund.
Take whatever money is left over and invest 70% in a stock index fund and 30% in a bond fund through any discount broker and never touch it until retirement.
If any of this confuses you, or you have something special going on (retirement, college planning, tax issues) hire a fee-based financial planner, not one who charges a percentage of your portfolio.

That's it!

Now, if really you wanna try your hand at automated trading systems, whatever your motivation may be, use a demo or paper trading account as much as possible and educate yourself on the markets you will be trading (an MBA should be a good baseline). If you're a student and have a professor interested, you could use IB's student trading lab. The strategy below can actually be run against Interactive Brokers' own demo software (which does not require opening an account, just download the software), and to do that, you can follow Max's and Domenic's tutorials.

The strategy below is a mean-reverting strategy where we imagine the prices to "bounce" between upper and lower price bands, based on one from Ernest Chan, and which connects to IB's Trader Workstation application through their ActiveX API. Max's tutorial should guide you through all that. Now, unlike Ernest's code, I don't use the matlab2IB library, so as long as you have Matlab on a Windows box you should be good to go. It is a technical analysis (TA) trading strategy which uses Bollinger bands, and although TA sounds a little too much like astrology for my taste, the mechanics of an automated trading system are well illustrated.

The code is contained in 2 files, 'trader.m' and 'driver.m', and 'driver.m' is the one you execute (sorry about the crappy syntax highlighting... Matlab is not a choice in syntaxhighlighter):


%driver.m

%clear workspace variables and command window
clear;
clc;

%state holds info that will be persisted througout the execution,
%and that needs to be accessed from elsewhere (trader)
global state; 
state=containers.Map();

%exSeq denotes more or less the execution sequence
exSeq={'preInit','init','accountUpdates','accountUpdatesEnd',...
    'marketData','histData','updateHistData','placeOrders','exit'};
currStep='preInit';

%leave ip blank for localhost
%socket port as defined in tws app
%ib api webinar says up to 8 concurrent clients are supported,
%which can be achieved by increasing clientId
conInfo=struct('ipAddress','','port',7496,'clientId',0);

%pacific standard time
timeZone='PST';

%create instance of tws activex control and pass event handler 'trader'
tws=actxcontrol('TWS.TwsCtrl.1',[0 0 0 0],figure('Visible','off'),...
    'trader');

%all output should include as prefix this anonynous function
dispPrefix=@() [datestr(now,13),' driver: '];

%connect to tws app
disp([dispPrefix(),'connecting...']);
tws.connect(conInfo.ipAddress,conInfo.port,conInfo.clientId);

while true
    switch currStep
        
        case 'preInit'
            %here we just get confirmation of the connection and 
            %some useful data such as the nextValidId
            disp([dispPrefix(),'in preInit...']);
            
            if isKey(state,'nextValidId')
                %this is the condition that tells us we're good
                %to move to the next step (see nextValidId event 
                %handled by trader)
                currStep='init';
            end
            
        case 'init'         
            %here we add any kind of initialization or configuration
            %parameters for this strategy
            disp([dispPrefix(),'in init...']);   
            
            params=struct();
            
            %symbol for the front contract
            params.localSymbol='ESU9';
            state('localSymbol')=params.localSymbol;
            
            %IB symbol for this future
            params.symbol='ES';
            state('symbol')=params.symbol;
            
            params.securityType='FUT';
            params.exchange='GLOBEX';
            params.currency='USD';
            
            %desired order size
            params.numberOfContracts=1;
            
            %number of ticks used to compute the middle band, 
            %a N-period simple moving average (MA)
            params.N=20;
            
            %the upper band at KA times an N-period standard deviation 
            %above the middle band (MA+KA*sigma)
            params.KA=.5;
            
            %the lower band at KB times an N-period standard deviation 
            %below the middle band (MA-KB*sigma)
            params.KB=.5;    
            
            %KA and KB are generally 2 according to wkipedia, but that all
            %depends on the strategy... this is one value you can play with
            %to see more action happen...
            
            currStep='accountUpdates';
            
        case 'accountUpdates'
            %request account updates
            disp([dispPrefix(),'in accountUpdates...']); 
            
            %the particular event that this will kick off is the 
            %'updatePortfolioEx' event, which will tell us if we
            %currently have any positions in the security we'll be trading
            tws.reqAccountUpdates(1,'');  
            
            currStep='accountUpdatesEnd';
            
        case 'accountUpdatesEnd'
            %wait until we get confirmation that those were all
            %available account updates
            disp([dispPrefix(),'in accountUpdatesEnd...']);  
            
            if isKey(state,'doneAccountUpdates')
                %this batch of updates has completed, so let's ask to
                %stop receiving account position updates;
                %even though tws' activex api doc says to use '2' to stop,
                %'0' is what seems to work (1=subscribe, 0=unsub?)
                tws.reqAccountUpdates(0,'');
                
                %set next step to be executed
                currStep='marketData';
                
                %if there was nothing for tws to nitify us of, 
                %then we need to initialize the value in state here
                if ~isKey(state,'pos')
                    state('pos')=0;
                end
                
                disp([dispPrefix(),'starting position in security: ',...
                    num2str(state('pos'))]);
            end
            
        case 'marketData'
            %here we are gonna ask tws to start sending us market data
            disp([dispPrefix(),'in marketData...']);   
            
            %tws has these factory (create*) methods for all the
            %COM objects we need to pass as arguments
            contract=tws.createContract();
            
            %string: this is the symbol of the underlying asset.
            contract.symbol=params.symbol;
            
            %string: this is the security type. Valid values are:
            %STK, OPT, FUT, IND, FOP, CASH, BAG
            contract.secType=params.securityType;  
            
            %string: the expiration date. Use the format YYYYMM if
            %applicable
            contract.expiry='';  
            
            %double: the strike price
            contract.strike=0; 
            
            %string: specifies a Put or Call. 
            %Valid values are: P, PUT, C, CALL.
            contract.right='';   
            
            %string: allows you to specify a future or option contract 
            %multiplier. This is only necessary when multiple 
            %possibilities exist.
            contract.multiplier='';    
            
            %string: the order destination, such as Smart.
            contract.exchange=params.exchange; 
            
            %string: specifies the currency. Ambiguities may require that
            %this field be specified, for example, when SMART is the
            %exchange and IBM is being requested (IBM can trade in
            %GBP or USD). Given the existence of this kind of
            %ambiguity, it is a good idea to always specify the currency.
            contract.currency=params.currency;       
            
            %string: this is the local exchange symbol of the 
            %underlying asset.
            contract.localSymbol=params.localSymbol;
            
            %string: identifies the listing exchange for the 
            %contract (do not list SMART).
            contract.primaryExch=params.exchange;  
            
            %integer: if set to true, contract details requests and 
            %historical data queries can be performed pertaining to 
            %expired contracts.
            contract.includeExpired=0;  
            
            %object: dynamic memory structure used to store the leg 
            %definitions for this contract.
            contract.comboLegs='';  
            
            %integer: the unique contract identifier.
            contract.conId=0;

            %request market data ('0' indicates we do not want a snapshot,
            %but wanna continue to receive a stream of data)
            tws.reqMktDataEx(contract.conId,contract,'',0);
            
            currStep='histData';
            
        case 'histData'
            %request historical data
             disp([dispPrefix(),'in histData...']);
             
             if ~isKey(state,'primedHistData')
                 %we only need to request historical data once, after that,
                 %we'll just keep updating it with market data
                 
                 %request historical data
                 tws.reqHistoricalDataEx(contract.conId,contract,...
                     [datestr(now,'yyyymmdd HH:MM:SS'),' ',timeZone],...
                     '1 D','1 min','TRADES',1,2);
                 
                 %start tic from the time we request all initial histData
                 tic;
                 
                 state('primedHistData')=true;
             else
                 if isKey(state,'doneHistData')
                     %we've already requested historical data, 
                     %so no need to do it again, we're now ready to 
                     %process this data
                     currStep='updateHistData';
                 else
                     %do nothing, still waiting on all the data,
                     %which is getting handled by the 'historicalData'
                     %event
                 end
             end
             
        case 'updateHistData'
            %update our historical data based on the market 
            %data that we subscribed to earlier
            disp([dispPrefix(),'in updateHistData...']);
            
            elapsed=toc;
            %every X elapsed seconds, get last tick data and update
            %our price history vector, keeping the vector's
            %length constant, once there is an update,
            %then we move forward to placeOrders
            
            %the elapsed time thershold is also another value you can
            %play with to see more action happen...
            if (elapsed >= 5)
                tic;
                histDataVec=state('histDataVec');
                histDataVec(1:end-1)=histDataVec(2:end);
                
                %if by chance we haven't received any market data yet,
                %this call below should fail
                histDataVec(end)=state('last');
                
                state('histDataVec')=histDataVec;
                
                %calculate moving average
                ma=mean(histDataVec(end-params.N+1:end));
                
                %calculate moving standard deviation
                mstd=std(histDataVec(end-params.N+1:end));
                
                %calculate deviation of ask and bid prices from the
                %moving average middle band.  All we're doing here is
                %solving for K on the 'ask' and 'bid' prices:
                %ask=(MA+K*sigma)
                askK=(state('ask')-ma)/mstd;
                %bid=(MA+K*sigma)
                bidK=(state('bid')-ma)/mstd;
                
                %move onto next step
                currStep='placeOrders';
            end
            
        case 'placeOrders'
            %place orders if the signals we're watching tells us to do so
            disp([dispPrefix(),'in placeOrders...']);
            
            %here we decide whether to exit or go back to updateHistData,
            %which could be based on time (trading during trading hours,
            %or a predefined number of trades); we illustrate with a 
            %predefined number of trades
            if ~isKey(state,'placeOrdersCount')
                state('placeOrdersCount')=1;
                currStep='updateHistData';
            else
                if state('placeOrdersCount') > 50
                    currStep='exit';
                else
                    state('placeOrdersCount')=state('placeOrdersCount')+1;
                    currStep='updateHistData';
                end
            end
            
            pos=state('pos');
            orderId=state('nextValidId');
            order=tws.createOrder();
            order.orderId=orderId;
            order.clientId=conInfo.clientId;
            order.action=''; %need below!
            order.totalQuantity=0; %need below!
            order.orderType='MKT';
            order.transmit=1;
            order.ocaGroup='';
            order.lmtPrice=0;
            order.auxPrice=0;
            order.timeInForce='DAY';
            
            disp([dispPrefix(),...
                'pos: ',num2str(pos),...
                ', askK: ',num2str(askK),...
                ', bidK: ',num2str(bidK),...
                ', ma: ',num2str(ma),...
                ', bid: ',num2str(state('bid')),...
                ', ask: ',num2str(state('ask'))]);
            
            if (pos == 0 && askK < -params.KA) 
                
                %if we are below the lower band and have nothing, go long
                
                order.action='BUY';
                order.totalQuantity=params.numberOfContracts;
                tws.placeOrderEx(orderId,contract,order);
                disp([dispPrefix(),'placed entry only...']);
                pos=params.numberOfContracts;
                
            elseif (pos < 0 && askK < -params.KA) 
                
                %if we are below the lower band and are short, buy 2,
                %one for establishing the long position and 
                %another one to cover the short
                
                order.action='BUY';
                order.totalQuantity=2*params.numberOfContracts;
                tws.placeOrderEx(orderId,contract,order);
                disp([dispPrefix(),'placed exit and entry...']);
                pos=params.numberOfContracts;
                
            elseif (pos == 0 && bidK > params.KB) 
                
                %if we are above the upper band and have nothing, go short
                
                order.action='SELL';
                order.totalQuantity=params.numberOfContracts;
                tws.placeOrderEx(orderId,contract,order);
                disp([dispPrefix(),'placed entry only...']);
                pos=-params.numberOfContracts;
                
            elseif (pos > 0 && bidK > params.KB) 
                
                %if we're above the upper band and are currently long,
                %sell 2, one for establishing the short position and 
                %anothe rone to take profits on the previous long
                
                order.action='SELL';
                order.totalQuantity=2*params.numberOfContracts;
                tws.placeOrderEx(orderId,contract,order);
                disp([dispPrefix(),'placed exit and entry...']);
                pos=-params.numberOfContracts;
                
            end
            
            state('pos')=pos;
            state('nextValidId')=orderId+1;
            
        case 'exit'
            %get ready to leave
            disp([dispPrefix(),'in exit...']);
            
            %maybe should cancel marketdata request, but we'll let
            %disconnecting take care of that for now
            pos=state('pos');
            orderId=state('nextValidId');
            order=tws.createOrder();
            order.orderId=orderId;
            order.clientId=conInfo.clientId;
            order.action=''; %need below!
            order.totalQuantity=0; %need below!
            order.orderType='MKT';
            order.transmit=1;
            order.ocaGroup='';
            order.lmtPrice=0;
            order.auxPrice=0;
            order.timeInForce='DAY';
            
            %should really be looking and waiting for order confirmations,
            %but the pauses will do the job for now...
            if (pos > 0)
                %cover the longs
                order.action='SELL';
                order.totalQuantity=pos;
                tws.placeOrderEx(orderId,contract,order);
                disp([dispPrefix(),'placed final order...']);
                pause(5);
            elseif (pos < 0)
                %cover the shorts
                order.action='BUY';
                order.totalQuantity=abs(pos);
                tws.placeOrderEx(orderId,contract,order);
                disp([dispPrefix(),'placed final order...']);
                pause(5);
            end
            
            pos=0;
            state('pos')=pos;
            state('nextValidId')=orderId+1;
            
            break;
            
        otherwise
            disp([dispPrefix(),'invalid currStep: ',state('currStep')]);
            
    end
    
    %the pause gives the activex component events a chance to be handled
    pause(2);
end

%disconnect from the tws app 
tws.disconnect();
disp([dispPrefix(),'disconnecting...']);
pause(2);

This is 'trader.m':


%trader.m

function trader(varargin)

%as per matlab's documentation on arguments passed to event handlers,
%'end-1' is the position the event payload will be at
arg=varargin{end-1};

%very little code should go in the switch statement,
%all logic should be dispatched to subfunctions on this file
switch arg.Type     
    case 'errMsg'
        %EVERY message should be well understood
        custom_disp(['in errMsg: ',arg.errorMsg]);   
    case 'updateAccountValue'
       %no need to do anything here
    case 'updateAccountTime'
        %no need to do anything here
    case 'updatePortfolioEx'
        updatePortfolioEx(arg);
    case 'accountDownloadEnd'
        accountDownloadEnd();
    case 'nextValidId'
        nextValidId(arg);
    case 'openOrderEnd'
        custom_disp('in openOrderEnd');
    case 'tickPrice'
        tickPrice(arg);
    case 'tickSize'
        %no need to do anything here
    case 'tickString'
        %no need to do anything here
    case 'tickGeneric'
        %no need to do anything here
    case 'historicalData'
        historicalData(arg);
    case 'updatePortfolio'
        %this guy has been replaced by updatePortfolioEx, so
        %no need to do anything here
    otherwise
        custom_disp(['unknown event: ',arg.Type]);    
end
end

function custom_disp(textMessage)
disp([datestr(now,13),' trader: ',textMessage]);
end

function nextValidId(arg)
%we get this message soon after having established a connection
custom_disp('in nextValidId');
global state;
state('nextValidId')=arg.id;
end

function updatePortfolioEx(arg)
custom_disp('in updatePortfolioEx');
global state;
if strcmp(arg.contract.symbol,state('symbol')) &&...
        strcmp(arg.contract.localSymbol,state('localSymbol'))
    state('pos')=arg.position;
end
end

function accountDownloadEnd()
custom_disp('in accountDownloadEnd');
global state;
state('doneAccountUpdates')=true;
end

function tickPrice(arg)
%as we get updates for this guy, shove the info into our state
custom_disp('in tickPrice');
global state;
switch arg.tickType
    case 1 %bid
        state('bid')=arg.price;
    case 2 %ask
        state('ask')=arg.price;
    case 4 %last
        state('last')=arg.price;
end
end

function historicalData(arg)
%all we're doing here is gathering the historical data messages
%into a vector
%custom_disp('in historicalData');
global state;

if ~isKey(state,'histDataVec')
    %create the vector and position indicator 
    
    %we shouldn't get more than a day's worth of data
    state('histDataVec')=zeros(24*60,1,'double');
    state('histDataVecN')=0;
end

%-1 is shown on the last message in the histData request
if (arg.close == -1)
    %resize the data vector to the number of messages we actualy received
    state('doneHistData')=true;
    histDataVec=state('histDataVec');
    N=state('histDataVecN');
    state('histDataVec')=histDataVec(1:N,1);
else
    %Matlab arrays have pass-by-value semantics, so this thing
    %that we're doing here is very inneficient (the array is getting
    %cloned everytime).  Just a heads up...
    
    %keep track of the pieces we're interested in
    histDataVec=state('histDataVec');
    N=state('histDataVecN')+1;
    histDataVec(N,1)=arg.close;
    state('histDataVec')=histDataVec;
    state('histDataVecN')=N;
end

end

And the output window with all the tracing would look something like this:


14:14:14 driver: connecting...
14:14:15 driver: in preInit...
14:14:15 trader: in nextValidId
14:14:15 trader: in openOrderEnd
14:14:15 trader: in errMsg: Market data farm connection is OK:ibdemo
14:14:15 trader: in errMsg: HMDS data farm connection is OK:demohmds
14:14:17 driver: in preInit...
14:14:19 driver: in init...
14:14:21 driver: in accountUpdates...
14:14:22 trader: in accountDownloadEnd
14:14:23 driver: in accountUpdatesEnd...
14:14:23 driver: starting position in security: 0
14:14:25 driver: in marketData...
14:14:25 trader: in tickPrice
14:14:25 trader: in tickPrice
14:14:25 trader: in tickPrice
14:14:25 trader: in tickPrice
14:14:25 trader: in tickPrice
14:14:25 trader: in tickPrice
14:14:25 trader: in tickPrice
14:14:26 trader: in tickPrice
14:14:26 trader: in tickPrice
14:14:27 driver: in histData...
14:14:29 trader: in tickPrice
14:14:29 driver: in histData...
14:14:31 driver: in updateHistData...
14:14:32 trader: in tickPrice
14:14:33 driver: in updateHistData...
14:14:35 driver: in placeOrders...
14:14:35 driver: pos: 0, askK: -0.62025, bidK: -1.5063, ma: 966.1, bid: 965.25, ask: 965.75
14:14:35 driver: placed entry only...
14:14:36 trader: in tickPrice
14:14:37 trader: unknown event: openOrder1
14:14:37 trader: unknown event: openOrder2
14:14:37 trader: unknown event: openOrder3
14:14:37 trader: unknown event: openOrder4
14:14:37 trader: unknown event: openOrderEx
14:14:37 trader: unknown event: orderStatus
14:14:37 trader: in tickPrice
14:14:37 trader: in tickPrice
14:14:37 trader: unknown event: execDetails
14:14:37 trader: unknown event: execDetailsEx
14:14:37 trader: unknown event: openOrder1
14:14:37 trader: unknown event: openOrder2
14:14:37 trader: unknown event: openOrder3
14:14:37 trader: unknown event: openOrder4
14:14:37 trader: unknown event: openOrderEx
14:14:37 trader: unknown event: orderStatus
14:14:37 trader: unknown event: openOrder1
14:14:37 trader: unknown event: openOrder2
14:14:37 trader: unknown event: openOrder3
14:14:37 trader: unknown event: openOrder4
14:14:37 trader: unknown event: openOrderEx
14:14:37 trader: unknown event: orderStatus
14:14:37 driver: in updateHistData...
...
14:19:45 trader: in tickPrice
14:19:45 trader: in tickPrice
14:19:45 trader: unknown event: execDetails
14:19:45 trader: unknown event: execDetailsEx
14:19:45 trader: unknown event: openOrder1
14:19:45 trader: unknown event: openOrder2
14:19:45 trader: unknown event: openOrder3
14:19:45 trader: unknown event: openOrder4
14:19:45 trader: unknown event: openOrderEx
14:19:45 trader: unknown event: orderStatus
14:19:45 trader: unknown event: openOrder1
14:19:45 trader: unknown event: openOrder2
14:19:45 trader: unknown event: openOrder3
14:19:45 trader: unknown event: openOrder4
14:19:45 trader: unknown event: openOrderEx
14:19:45 trader: unknown event: orderStatus
14:19:45 trader: in tickPrice
14:19:46 trader: in tickPrice
14:19:47 trader: in tickPrice
14:19:48 trader: in tickPrice
14:19:49 trader: in tickPrice
14:19:49 driver: disconnecting...

The code is pretty self-explanatory, and is a nice baseline to get started on some simple strategies. Anything more complicated than that and I would recommend some sort of state pattern implementation.

I ran it against the demo account and despite the fact that my parameters were optimized to minimize my waiting time (very low params.KA and params.KB and an arbitrary number of cycles before I liquidated any current positions, not to mention the fact that only market orders were used), it only lost about $100 (you start with about $500,000 in the demo account and each transaction cost about $1,000). This shouldn't by any means sound encouraging, since those prices aren't even real, they're simulated.

Anyways, this should give you an idea of the whole thing. You may wanna look into what futures and e-minis are all about or change the strategy to trade some stock that you may be familiar with.

Another good set of resources are IB's webinars on their APIs and the sample applications they provide that hook up to their APIs.

Enjoy!

Interactive Brokers, ActiveX & Matlab

2009-08-17T11:53:00.000-07:00

So, last time I mentioned a few useful sites with info about algorithmic trading in general and about using matlab as a trading platform. Ernest Chan, in particular, offers a wealth of information on his blog and on his website, which includes a premium content section that you can access if you have a copy of his quantitative trading book.

There's some good stuff in that premium content section, including an article mentioned on a post of his that discusses a matlab implementation of a technical analysis trading strategy (using Bollinger bands). The strategy is simple and illustrative enough, but the problem is that if you want to execute it and play around with it, you need the matlab2ib library to hook up to interactive brokers. Fortunately, the creators of the matlab2ib library (which appears to be a thin shell around IB's own java API) offer a free 30-day trial version of their product; unfortunately however, this limited-time trial version is also limited in functionality and happens to lack certain functions used in the Bollinger band example.

What I wanted was a quick way to get a feel for the workflow in the Trader Workstation application and its API so that I could then nicely structure a lightweight and reusable template to try out a few of my own strategies. As I mentioned in my previous post, I have an interest in becoming more than simply familiar with matlab, and so my platform was decided. On the other hand, R is particularly attractive for this since it offers for free much of the functionality that you would otherwise have to pay for in matlab. For example, matlab offers the following relevant toolboxes:

Which is all really impressive stuff (I mean, they've got everything!), but if you're operating on a budget and have the know-how, a little elbow grease goes a long way in R:

Furthermore, there is even a package for connecting R to IB, IBrokers, which looks very much alive, although last time I looked at it there was no order handling functionality yet.

Another possible venue could have been python, which also has IB connectivity through ibpy, a library interesting on its own merits due to the approach used by its creator, Troy Melhase, who wrote a java2python utility for converting IB's java api to python source using ANTLR.

Personally, though, as much as I'm tempted otherwise, at this point I wanna mess with infrastructure stuff as little as possible and get to the meat of the system, the trading strategies, as quickly as possible. The approach that seemed to satisfy that requirement best was to use IB's ActiveX API from matlab, so I resurrected my dusty old windows box and wrote this little post in the time that it took for it to boot up :)

Next time we'll actually look at some source and explore Trader Workstation's workflow.

Automated Trading Systems

2009-08-09T14:53:00.001-07:00

An automated trading system seems like a good project in which to bring together all the ideas we've been discussing in this blog so far, perhaps with an aim to participate in the trading competition Interactive Brokers (IB) has been holding and will hopefully hold again this year.

Matlab appears to be a widely popular choice for these kinds of things, and plus it won't hurt to get some experience on that platform since it is also pretty popular in the engineering classes at UC Santa Cruz, where I'll be transferring later this year. Classes I'm particularly looking forward to this fall include:

And I'm sure they'll also be the subject of a few blog posts. Back on the topic of automated trading systems, here are a few links that can help jumpstart the project (video tutorials, code, etc) if you're also interested in building such a system:

Paper trading is a nice, safe way to get started, but one of the really interesting things has to be being able to move the market based on one's trading activity (the difference between observing a markov chain and influencing it in a markov decision problem). IB offers a paper trading account to their customers.

The following quote from Assistant U.S. Attorney Joseph Facciponti acting on behalf of Goldman Sachs in the recent snafu where an ex-Goldman employee, Sergey Aleynikov, allegedly uploaded proprietary trading software from Goldman to a server in Germany prior to moving to a rival company, sheds some light on the otherwise extremely secretive industry:

"The bank has raised the possibility that there is a danger that somebody who knew how to use this program could use it to manipulate markets in unfair ways..."

My irony meter exploded when I read that...

Anyways, there are a lot of other venues in which to exert control over markets including market making (see here) and how the Fed attempts to influence the greater economy by manipulating very short term interest rates (here's a video resource on the Fed's open market operations). I wonder the extent to which such control mechanisms are automated.

Cognitive Reflection Test

2009-08-07T11:11:00.001-07:00

Speaking of control and optimal strategy, more often than not, the correct [or optimal] answer isn't the most obvious one. Nathan Rothschild gave us an example of that, and I am also reminded of a study that makes an interesting analysis based on an individual's response to the following questions:

A bat and a ball cost $1.10 in total. The bat costs a dollar more than the ball. How much does the ball cost?
If it takes 5 machines 5 minutes to make 5 widgets, how long would it take 100 machines to make 100 widgets?
In a lake, there is a patch of lily pads. Every day, the patch doubles in size. If it takes 48 days for the patch to cover the entire lake, how long would it take for the patch to cover half of the lake?

Try to answer them before reading further...

In the study, Shane Frederick calls the 3-question test the "Cognitive Reflection Test" (CRT), and claims that it measures the "ability or disposition to resist reporting the response that first comes to mind." You will have to read the paper for the correct answers and for more details, but briefly, the nature of the questions is such that if you resist the temptation to blurt out the obvious and ponder the question for a little bit, you will almost inevitably arrive at the correct answer. To a lot of people, however, investing the time in thinking about the problem can be a painful experience, and the alternative, a quick and apparently good enough answer, can provide a relief to this pain.

The study goes further and shows that performance on the CRT is correlated with an individual's ability or tendency to forgo immediate reward in place of a higher future reward, a quality that can be essential in optimal decision making.

Of course it helps if you have seem similar problems before and are able to reason by analogy; however, there is the implication that you have [willingly] endured through those similar problems before and learned the pattern. Clearly, there is a powerful feedback mechanism at work here, were your good performance increases the likelihood you will practice your skill, which will in turn improve your performance, and so on. Some believe that this mechanism (and a few other important details such as a drive, favorable environment, good socio-economic circumstances, etc), if exploited well, is one of the main factors in great performance (genius-like work) vs average performance (see Geoff Colvin's Talent is Overrated), and I tend to agree with that (most people have the same foundational ability but there is a time element, similar to how continuous compounding works, that needs to work with this feedback mechanism).

Anyways, the study also mentions an interesting question:

What would you prefer, a sure $500 or a 15% chance of $1,000,000?

There was also correlation of people's performance on the CRT with their answer to this question, specifically, more people who scored well on the CRT chose the 15% gamble, whereas most people who did not score so well chose on the CRT chose the $500.

Decision theory tells us that the correct way to answer this question is to calculate the expected value of each choice:

$500 * 1.0 = $500
$1,000,000 * 0.15 = $150,000

And the thing is that after playing this particular game enough times, having always chosen the second option as a policy will have yielded a higher profit. What if you only play the game once? Well, then there can be some controversy, and that's why some people see this particular question (one-shot game) as a case where personal preference is what matters (what's better: blue or green?) and there being no necessarily wrong answers; however, if one looks at this situation as an indicative of an individual's general tendency towards choice under uncertainty throughout the course of their life, even if they never face this exact problem again, having this general policy and decision-making strategy will statistically yield higher profits/rewards.

This kind of thinking is one of the strengths of reinforcement learning methods and bayesian influence networks, which when automated (given a proper model), can work much faster and objectively than any human today possibly could. Such simulation systems may also be more likely to identify unintended consequences.

Now, what happens when the decisions we make influence other people's decisions, which in turn may influence our decisions? How do we find optimal strategies then? Well, game theory has a lot to say about that and it may very well be the subject of a future post.

Time Series Modeling, Forecasting & Control II

2009-08-04T16:27:00.001-07:00

So, continuing with our discussion on time series analysis, I would like to focus on what I consider to be the end-goal of such enterprise: control. Modeling and forecasting (I think 'estimating' would be a better, more encompassing term than 'forecasting') are just things we need to do in order to be able to control.

By control, I mean influencing a system (entirely or partially) so as to maximize some notion of reward received from that system. There is nothing really new here, these ideas have been the core of a several related fields:

reinforcement learning: "reinforcement learning is a sub-area of machine learning concerned with how an agent ought to take actions in an environment so as to maximize some notion of long-term reward."
optimal control: "optimal control deals with the problem of finding a control law for a given system such that a certain optimality criterion is achieved."
operations research: "concerned with determining the maxima (of profit, assembly line performance, crop yield, bandwidth, etc) or minima (of loss, risk, etc.) of some objective function."
artificial intelligence: "the study and design of systems that perceive their environment and take actions which maximize their chances of success."

Now, given all these fields and the countless minds which have been dedicated to solving what is essentially the same problem, have we found a solution yet? Or have we at least found some sort of general framework or unifying theme that may help guide us in our search? Well, unfortunately (or fortunately?), a perfect solution still evades us, but I think we have a pretty good general framework for developing good enough (and sometimes optimal) solutions: bayesian networks.

The term 'bayesian networks' may not be the best (perhaps 'probabilistic graphical models' is better, more general), but then again, what's in a name? In any event, I think bayesian networks are the best general framework we have for dealing with these problems (if you know of a better one, do let me know), and this post will illustrate that a bit. This is not a tutorial or anything like that, rather my goal is to pique your curiosity enough so that you may decide to look further into bayesian networks.

Let's start with some motivating problems:

There are eight types of blood (A, B, AB and O, each coming in + or -), which may be stored until it is six weeks old, at which point it has to be discarded. Some blood types may be substituted for others (as shown in the diagram). For example, O- blood may be used for anyone, while B+ blood may be used only with patients with B+ blood or AB+ blood.
Each week, new donations are made, and new demands for blood arise. The problem is to determine which blood type to use for each demand type in a given week. Do you use the O- blood now, or hold it for later when it might be used for a patient with A- blood, which can only use A- blood or O- blood?
A pig breeder is growing pigs for a period of four months and subsequently selling them. During this period the pig may or may not develop a certain disease. If the pig has the disease at the time it must be sold, the pig must be sold for slaughtering, and its expected market price is then 300 DKK (Danish kroner). If it is disease free, its expected market price as a breeding animal is 1000 DKK.

Once a month, a veterinary doctor sees the pig and makes a test for presence of the disease. If the pig is ill, the test will indicate this with probability 0.80, and if the pig is healthy, the test will indicate this with probability 0.90. At each monthly visit, the doctor may or may not treat the pig for the disease by injecting a certain drug. The cost of an injection is 100 DKK.

A pig has the disease in the first month with probability 0.10. A healthy pig develops the disease in the subsequent month with probability 0.20 without injection, whereas a healthy and treated pig develops the disease with probability 0.10, so the injection has some preventive effect. An untreated pig that is unhealthy will remain so in the subsequent month with probability 0.90, whereas the similar probability is 0.50 for an unhealthy pig that is treated. Thus spontaneous cure is possible, but treatment is beneficial on average. The problem here is to come up with an optimal course of action course of action, a strategy, which maximizes the expected utility for the pig breeder.

How would you solve these 2 problems? Is there a method that can find a satisfactory answer that can be applied to both problems? I think you know the answer already :)

The blood management problem is from Warren Powell's book "Approximate dynamic programming: solving the curses of dimensionality" and there is some material online about it, including a matlab implementation of an approximate dynamic programming solution. The second problem is from a paper by Steffen Lauritzen, and there is an an implementation of the algorithm proposed by Steffen in the matlab toolbox for bayesian networks developed by Kevin Murphy, which includes as example a solution to the pig problem.

Granted, methods such as those under the umbrella of reinforcement learning and approximate dynamic programming are extremely powerful, maybe as powerful as bayesian networks, however there appears to be a general consensus in the research community that regards having a more explicit model of the environment as a generally superior approach. To fuel the controversy (or confusion, really, since we really just don't quite know which is in fact better: model-free or model-based), there is a place for models in certain reinforcement learning architectures; in fact, there is a varying degree of model complexity in the wide spectrum of reinforcement learning techniques out there.

Personally, and this is purely qualitative speculation on my part, I think that we could get more mileage out of using bayesian networks as the framework onto which we integrate reinforcement learning ideas rather than the other way around (Lauritzen's Single Policy Updating can be thought of as a sort of Policy Iteration used for solving Markov Decision Problems). Why do I think so? Well, for one, I think just about every machine learning technique we know of can be framed in a bayes net perspective, not just those that require sequential decision making and trading-off rewards in the present vs rewards in the future. Furthermore, there is a variety of reasoning patterns that are naturally facilitated by the bayes net formulation:

diagnostic reasoning: reasoning from symptoms to causes
predictive reasoning: reasoning from new information about causes to new beliefs about effects
intercausal reasoning: reasoning about mutual causes of a common effect

Yes, this 'reasoning patterns' taxonomy can be a little sketchy if looked at too closely, specially if the reader is in any way 'allergic' to notions of causality such as those discussed by Pearl, but when we get down to it, they are simply computations we do on a joint probability distribution; and tell me, do we have a more efficient way to represent and compute on high-dimensional joint probability distributions than bayesian networks?

So, what's the best way to jump into bayesian nets? Well, assuming you've gone through AIMA, I would recommend starting with Neapolitan's Learning Bayesian Networks. As you go through it, download and play with the bayes net toolbox for matlab, a preview of which is included below (solution for the pig breeder problem):


% pigs model from Lauritzen and Nilsson, 2001
% from BNT/examples/limids

seed = 0;
rand('state', seed);
randn('state', seed);

% we number nodes down and to the right
h = [1 5 9 13];
t = [2 6 10];
d = [3 7 11];
u = [4 8 12 14];

N = 14;
dag = zeros(N);

% causal arcs
for i=1:3
 dag(h(i), [t(i) h(i+1)]) = 1;
 dag(d(i), [u(i) h(i+1)]) = 1;
end
dag(h(4), u(4)) = 1;

% information arcs
fig = 2;
switch fig
case 0,
 % no info arcs
case 1,
  % no-forgetting policy (figure 1)
  for i=1:3
    dag(t(i), d(i:3)) = 1;
  end
case 2,
 % reactive policy (figure 2)
 for i=1:3
   dag(t(i), d(i)) = 1;
 end
case 7,
 % omniscient policy (figure 7: di has access to hidden state h(i-1))
 dag(t(1), d(1)) = 1;
 for i=2:3
   %dag([h(i-1) t(i-1) d(i-1)], d(i)) = 1;
   dag([h(i-1) d(i-1)], d(i)) = 1; % t(i-1) is redundant given h(i-1)
 end
end


ns = 2*ones(1,N);
ns(u) = 1;

% parameter tying
params = ones(1,N);
uparam = 1;
final_uparam = 2;
tparam = 3;
h1_param = 4;
hparam = 5;
dparams = 6:8;

params(u(1:3)) = uparam;
params(u(4)) = final_uparam;
params(t) = tparam;
params(h(1)) = h1_param;
params(h(2:end)) = hparam;
params(d) = dparams;

limid = mk_limid(dag, ns, 'chance', [h t], 'decision', d, 'utility', u, ...
               'equiv_class', params);

% h = 1 means healthy, h = 2 means diseased
% d = 1 means don't treat, d = 2 means treat
% t = 1 means test shows healthy, t = 2 means test shows diseased

if 0
 % use random params
 limid.CPD{final_uparam} = tabular_utility_node(limid, u(4));
 limid.CPD{uparam} = tabular_utility_node(limid, u(1));
 limid.CPD{tparam} = tabular_CPD(limid, t(1));
 limid.CPD{h1_param} = tabular_CPD(limid, h(1));
 limid.CPD{hparam} = tabular_CPD(limid, h(2));
else
 limid.CPD{final_uparam} = tabular_utility_node(limid, u(4), [1000 300]);
 % costs have negative utility!
 limid.CPD{uparam} = tabular_utility_node(limid, u(1), [0 -100]);

 % h  P(t=1) P(t=2)
 % 1  0.9   0.1
 % 2  0.2   0.8
 limid.CPD{tparam} = tabular_CPD(limid, t(1), [0.9 0.2 0.1 0.8]);

 % P(h1)
 limid.CPD{h1_param} = tabular_CPD(limid, h(1), [0.9 0.1]);

 % hi di P(hj=1) P(hj=2),  j = i+1, i=1:3
 % 1  1  0.8     0.2
 % 2  1  0.1     0.9
 % 1  2  0.9     0.1
 % 2  2  0.5     0.5
 limid.CPD{hparam} = tabular_CPD(limid, h(2), [0.8 0.1 0.9 0.5 0.2 0.9 0.1 0.5]);
end

% Decision nodes get assigned uniform policies by default
for i=1:3
 limid.CPD{dparams(i)} = tabular_decision_node(limid, d(i));
end

engines = {};
engines{end+1} = global_joint_inf_engine(limid);
engines{end+1} = jtree_limid_inf_engine(limid);
%engines{end+1} = belprop_inf_engine(limid, 'max_iter', 1*N, 'filename', ...
                                   %fname, 'tol', 1e-3);

exact = [1 2];
%approx = 3;
approx = [];

max_iter = 1;
order = d(end:-1:1);
%order = d(1:end);

NE = length(engines);
MEU = zeros(1, NE);
niter = zeros(1, NE);
strategy = cell(1, NE);
for e=1:NE
 [strategy{e}, MEU(e), niter(e)] = solve_limid(engines{e}, 'max_iter', ...
     max_iter, 'order',  order);
end
MEU

% check results match those in the paper (p. 22)
direct_policy = eye(2); % treat iff test is positive
never_policy = [1 0; 1 0]; % never treat
tol = 1e-0; % results in paper are reported to 0dp
for e=exact(:)'
 switch fig
  case 2, % reactive policy
   assert(approxeq(MEU(e), 727, tol));
   assert(approxeq(strategy{e}{d(1)}(:), never_policy(:)))
   assert(approxeq(strategy{e}{d(2)}(:), direct_policy(:)))
   assert(approxeq(strategy{e}{d(3)}(:), direct_policy(:)))
  case 1, assert(approxeq(MEU(e), 729, tol));
  case 7, assert(approxeq(MEU(e), 732, tol));
 end
end


for e=approx(:)'
 for i=1:3
   approxeq(strategy{exact(1)}{d(i)}, strategy{e}{d(i)})
   dispcpt(strategy{e}{d(i)})
 end
end

Which yields a strategy for a maximum expected utility of 727 DKK, a pretty good approximation to the optimal maximum expected utility of 729 DKK. You should really go read Lauritzen's paper, but just in case you're curious, here is what an optimal strategy would look like:

Never treat in the first month.
Treat in the second month if and only if test results are positive in the first month and in the second month.
Treat in the last month if and only if the test results are positive in the third month or in both of the first and second months.

An interesting exercise would be to now solve the blood management problem under this framework. We shall leave that for a future post.

Time Series Modeling, Forecasting & Control

2009-07-29T09:43:00.000-07:00

Most of the algorithms we've discussed so far are for data which is considered to be independently and identically distributed, or IID (we did look briefly at Markov models here and here); in other words, each data point is sampled from the same invariant distribution and where a given sample is independent from any samples that came before it. This is easily visualized by thinking of a data generating machine/process whose core engine is a Gaussian distribution. Once we turn the machine on and it begins to generate data, we may see something like this:


> rnorm(10,mean=10,sd=1)
[1] 10.685401  9.262564  9.965892  9.637854 11.069165  8.789038 10.151829
[8] 10.851841  9.615401 11.004101

What if we have a different kind of generating machine? Say we have a state machine, where each data point generated depends stochastically on the current state of the machine:


seed<-0
for (val in rnorm(10,mean=10,sd=1)) 
{
 seed<-seed+val
 cat(seed,"\n")
}
8.715963 
17.36178 
25.38322 
35.13913 
45.95719 
56.685 
66.37728 
75.36139 
85.55207 
95.67641

Or how about a machine that is evolving through time according to some deterministic (or nondeterministic) function, which really is just another kind of state machine where the state is a time counter:


samples<-rnorm(10,mean=10,sd=1)
for (i in 1:10) 
{
 sample<-samples[i]*i
 cat(sample,"\n")
}
11.94815 
21.08698 
24.32918 
40.33526 
48.37734 
73.4921 
62.65447 
70.49196 
84.39426 
83.01182

In all these cases, we would not be well served by the IID assumption when trying to infer the generating process from data. What we do instead is to incorporate the sequential data dependency into our models. There are several standard ways of doing so, the most common, in order of increasing generality, being:

White Noise Processes
y_t = e_t; where e_t has a certain invariant distribution such as a Normal distribution with mean 'mu' and standard deviation 'sd'.

Moving Average (MA) Processes
y_t = e_t + w_1 * e_t-1; which is simply a white noise processes plus feedback from the noise at time t-1. This could be extended to also contain feedback from even more white noise signals from the past (t-2, t-3, etc).

Auto Regressive (AR) Processes
y_t = e_t + w_1 * y_t-1; where now we feed back not the previous noise, but the previous value of the series. This can be seen as a special case of an MA process where we add up all previous white noise terms. This model can also be generalized by feeding back more past values of the time series (t-2, t-3, etc).

Auto Regressive Moving Average (ARMA) Processes
y_t = e_t + w_1 * y_t-1 + w_2 * e_t-1; where we combine AR and MA processes, in order to seek better fit with fewer parameters than having a pure AR or MA model. We could also extend this model by adding further feedback components from time series values and/or noise.

Auto Regressive Integrated Moving Average (ARIMA) Processes
The previous models were for stationary processes, whereas ARIMA is a model for nonstationary processes. Take stock prices, for example, which follow a nonstationary process. If we subtract price at t-1 from price at t, we end up with another time series which may be stationary. If so, we could then apply some sort of ARMA model, say ARMA(p,q), which would be equivalent to saying we used an ARIMA(p,1,q) model (the 1 is for first difference). There are other variations, including SARIMA (seasonal) and VARIMA (vector), and they're all nicely structured under the Box-Jenkins framework.

All these models, once parameterized at their meta level (choosing 'p' and/or 'q'), can then be fit using maximum likelihood methods. Furthermore, rather than having a simple point-estimate prediction of the expectation, we can get a probabilistic distribution for the predicted value, in other words E(y|x) vs P(y|x). Speaking of which, I wonder what would be the result of introducing probabilities earlier in children's math curriculum, say as early as when they learn their multiplication tables. Then, they could be taught not that 3*5=15, but that E(3*5)=15. Perhaps I'm exaggerating a little, but just because we don't see it, it does not mean it is not there. These things pop up all the time in math (when dealing with 2-dimensional problems, the 3rd dimension is always lurking; or when we convert to the domain of complex numbers to solve a differential equation and then drop the imaginary part at the end); it is always helpful to learn to frame problems under new frameworks or imagine what's "behind" the paper, so to speak (projecting data points to a higher-dimensional space in order to linearly separate nonlinear data is just another example).

Anyways, all of the above can be done in R using the arima command and a variety of ancillary commands (acf, etc), but what I would like to explore further is a different kind of model for sequential data: linear dynamical systems, specifically, Kalman filters.

One obvious way to seek improvement in the temporal models described above is to include additional (exogenous) variables we think may have an effect on our variable of interest. For example, if we are trying to model the evolution of a stock price through time, we may add to the model, in addition to feedback lags from the stock price, feedback lags from macroeconomic variables such as inflation or gross national product, or maybe even feedback lags from the market price, resulting in a spatio-temporal model (an example of this is the so called ARMAX model).

We can think of the methods we've looked at so far as modeling time series having a certain (possibly zero) momentum component and a random component. Another way to visualize a time series is to think of it as a sequence of observed measurement from a process that transitions from state to state at each time interval, including the possibility that the state may not be completely observable. A Markov chain is an example of a dynamic, stationary process that evolves this way. A hidden Markov model is one in which the state is not observable (all we see are observable measurements) and is considered to be discrete-valued. The analogous to a hidden Markov model but with continuous latent variables is a linear dynamical system, which can be tracked/monitored/filtered and predicted with a Kalman filter, also known as a linear quadratic estimator (LQE).

That last name, LQE, is interesting because we can pair it with a linear quadratic regulator (LQR) and end up with an optimal solution to the linear quadratic gaussian control problem (LQG), which is what lectures 18 and 19 are mostly about. However, the really interesting part is that as useful and widely applicable as the Kalman filter is, it can be thought of as a special case of a dynamic bayesian network and so any improvements one may want to make to a Kalman filter can be accomplished by looking at it from the dynamic bayesian network perspective. We will look at these things further in a soon to follow part 2 of this post.

Capital Asset Pricing Model

2009-07-26T13:51:00.000-07:00

In our last post we looked at how Markowitz suggested we could optimize portfolios to maximize return and minimize risk. The optimal strategy turned out to be holding some combination of the tangency portfolio and the risk free asset (depending on one's specific risk profile). Naturally, we overlooked many details; however, the interested reader can refer to the following books for more information:

Rupert's book includes implementations in SAS and MATLAB, and is generally more accessible than Lai and Xing's Statistical Models and Methods for Financial Markets. In addition, here is a page that includes R code for most of the examples and illustrations in Rupert's.

Under certain idealized but plausible assumptions, elucidated as premises for the Capital Asset Pricing Model (CAPM), the market portfolio is the tangency portfolio. The S&P 500 being a commonly used proxy for the market portfolio implies that a reasonable strategy under this framework is to invest into some combination of some market index such as an ETF for the S&P 500 (SPY) and the risk free asset (Treasury bonds). Other indices include much more than just 500 companies, such as the Wilshire 5000.

There are 2 books from Burton Malkiel that further discuss these passive investment strategies and many other interesting bits of information, including:

An illustration of the difficulty of timing the market, first described by Charles Ellis (author of Winning the Loser's Game), where an investor who stayed invested in the S&P 500 index from 1982-2000 made an 18% return, whereas if that investor had missed the best 30 days in the market, the investor would have realized an 11% return.
Discussions on several market booms and busts, which should be enlightening for those who believe this current economic recession/depression is the first or the last one we will ever see. This should also be a good starting point for a reader interested in researching how previous crises were handle, internationally and in the US, particularly before the Federal Reserve System existed.

The two books are:

The second one is a more recent, summarized version of the first one, but the first one includes many interesting discussions on technical analysis, fundamental analysis, modern portfolio theory, behavioral finance, efficient market theory, etc.

Anyway, what I wanted to look at today was CAPM, and the perspective we will take is that of someone who is familiar with linear regression and wonders what insights can be gleaned from runing regression analysis to try and explain the return of a particular risky asset as a function of the market return (with some index as proxy for the market portfolio).

The code below shows that we will be using SPDR's for both the risky assets and the market index. Our risk free asset will be the 3-month treasury bill rate:


library(quantmod)
symbols<-c("XLY","XLP","XLE","XLF","XLV","XLI","XLB","XLK","XLU","SPY")
fromDate<-'1999-01-01'

#download all data available for the given symbols from yahoo finance
setDefaults(getSymbols.yahoo,from=fromDate) 
getSymbols(symbols,return.class="timeSeries")

#get 3-month t-bill rate
getSymbols("DTB3",return.class="timeSeries",src="FRED")

#gather all data into one multivariate time series, using only the
#close price, which is adjusted for splits and dividends
allAdj<-get(symbols[1])[,6]
for(symbol in symbols[-1]) allAdj<-cbind(allAdj,get(symbol)[,6])
allAdj<-cbind(allAdj,DTB3)

#adjust the column names
for(i in seq(along=symbols)) colnames(allAdj)[i]<-strsplit(colnames(allAdj)[i],".Adjusted")

#clean and save to disk
allAdj<-na.omit(allAdj)
write.csv(allAdj,file="data_spy.csv")

Once we have our data, we can run regression of a given risky asset's excess return against the market's excess return (excess from the risk free asset):


library(fPortfolio)
snpData<-readSeries(file="data_spy.csv",header=TRUE,sep=",")

rfree<-as.numeric(snpData$DTB3)
market<-as.numeric(snpData$SPY)
risky<-as.numeric(snpData$XLF)

market<-diff(log(market))-(rfree/(100*252))[2:length(rfree)]
risky<-diff(log(risky))-(rfree/(100*252))[2:length(rfree)]

fit<-lm(risky~market)
summary(fit)

Which yields:


Call:
lm(formula = risky ~ market)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.2019700 -0.0051045 -0.0002983  0.0046492  0.2107481 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.269e-05  2.774e-04  -0.298    0.766    
market       1.354e+00  1.972e-02  68.664   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.01424 on 2634 degrees of freedom
Multiple R-squared: 0.6416, Adjusted R-squared: 0.6414 
F-statistic:  4715 on 1 and 2634 DF,  p-value: < 2.2e-16

We end up with an intercept of -8.269e-05 and a slope of 1.354. Another metric of interest is R squared, for which we have a value of 0.6416. R squared can be thought of as a measure of how well the model explains the predicted data, and in this case, 1 - 0.64 = 0.36 means that 36% of the behavior of the excess returns of the risky asset are not explained by the market. So 36% is a measure of the amount of nonsystematic risk (sigma) in this particular risky asset (this can and should be diversified away in a good portfolio, which is another use of CAPM since it allows for a computationally efficient way to estimate the values needed in Markowitz Portfolio Optimization theory). The intercept and slope are also commonly known as alpha and beta under the CAPM framework, and together with the sigma, they are considered significant characteristics of a particular asset.

Under CAPM, a non-zero alpha would indicate the degree of mispricing in a given asset (if greater than zero underpriced, otherwise overpriced), which according to the efficient market hypothesis (an assumption of CAPM) should not happen (something that was confirmed in our little exercise). Finally, beta measures the asset's sensitivity to market movements (if the return on the market goes up by 1%, then the return on the asset would go up by beta%), so in CAPM, it is beta and not the variance in the returns that is the significant measure of [systematic] risk in the asset.

Although useful, CAPM has received a lot of heat, especially since it is so easy set up and test (you can read Malkiel's or Ruppert's book for details). In addition, there are extensions to the model, where, following the regression mindset, we simply add more regressors, giving rise to the multifactor pricing model under Arbitrage Pricing Theory, an example of which is the Fama-French three factor model.

For now, we will stick with the plain-vanilla CAPM and look at what happens when we recognize that beta is a time-varying quantity and how we may be able to track it. Until then, here is an interesting interview with Willian Sharpe, one of the creators of CAPM.

Markowitz Portfolio Optimization

2009-07-15T13:06:00.000-07:00

Having looked at a method for modeling stock returns, we can now look at how those ideas form the basis for a very significant contribution from Harry Markowitz to modern portfolio theory. In essence, investors want to:

Maximize returns, specifically, the mean of the return
Minimize risk, which we define to be the standard deviation of the return

And although these two goals are generally at odds with each other, Markowitz came up with a framework for structuring the problem so as to facilitate making the trade-off decision fairly objectively. Specifically, the framework gives us a way to mathematically handle something that we are all intuitively familiar with: the need for diversification.

We illustrate the ideas with an example, where an investor with a time horizon of 1 month is trying to find a portfolio composed of some combination of the 9 select sector SPDR ETF's that the S&P 500 has been divided into.

The data we need for our example can be found online for free (Yahoo! Finance), and downloading it directly from within R is a breeze using the quantmod library:


library(quantmod)
symbols<-c("XLY","XLP","XLE","XLF","XLV","XLI","XLB","XLK","XLU")

#download all data available for the given symbols from yahoo finance
setDefaults(getSymbols.yahoo,from='1900-01-01') 
getSymbols(symbols,return.class="timeSeries")

#gather all data into one multivariate time series, using only the
#close price which is adjusted for splits and dividends
allAdj<-get(symbols[1])[,6]
for(symbol in symbols[-1]) allAdj<-cbind(allAdj,get(symbol)[,6])

#adjust the column names
for(i in 1:9) colnames(allAdj)[i]<-strsplit(colnames(allAdj)[i],".Adjusted")

#save to disk
write.csv(allAdj,file="spdrAdj.csv")

Now, they key property that we want to exploit is that by mixing assets whose returns are not correlated, we can reduce the standard deviation of the return of the mix itself. In addition, we can do so without necessarily sacrificing our demands of a high expected return. The expectation (mean) is a linear operator and so expectations of linear combinations of expectations are easily computed. The variance, however, is a little different.

Having downloaded the data into our R session, we need to estimate the expectations and variance/covariances of the returns of the ETFs in our portfolio. Having that information, we then need to decided what percentage of our funds we will dedicate to each candidate ETF. This vector of weights assigned to each ETF is what defines the portfolio. In one scenario, the one we will consider in the code example, in addition to the sum of all weights needing to add up to one, the weights must also be positive, a so called long-only portfolio (an alternative scenario would consider portfolios where shorting is permitted).

The calculations for estimating the returns' expectations and their variances/covariances are easy. The hard part is making reliable assumptions about how stationary the time-series processes are, how much data to use, how robust our estimates are, etc. In this example, we make the very basic assumption that the processes are stationary, but in future posts we will consider more realistic (dynamic) models for drift and volatility. Optimizing the portfolio mix under a long-only constraint can be done numerically with quadratic programming, which R has a nice library for. When we do not have the long-only constraint, the optimization can be done analytically with the method of Lagrange multipliers. More details can be found on chapter 3 of Statistical Models and Methods for Financial Markets. A very useful R library that provides all this functionality out-of-the-box is RMetrics, which is what we use below:


library(fPortfolio)
library(fEcofin)
library(corpcor)

#the two functions simply fix a small problem in how RMetrics
#was plotting the tangency line in 'tangencyLines' (neither 
#the intercept nor the slope considered the risk free rate 
#if one was given)

myTangencyLines<-function(object, 
 return = c("mean", "mu"), 
 risk = c("Cov", "Sigma", "CVaR", "VaR"), 
 auto = TRUE, ...) 
{
    return = match.arg(return)
    risk = match.arg(risk)
    data = getSeries(object)
    spec = getSpec(object)
    constraints = getConstraints(object)
    tgPortfolio = tangencyPortfolio(data, spec, constraints)
    assets = frontierPoints(tgPortfolio, return = return, risk = risk, 
        auto = auto)
    slope = (assets[2]-getRiskFreeRate(spec))/assets[1]
    abline(getRiskFreeRate(spec), slope, ...)
    invisible(list(slope = slope, assets = assets))
}

myTailoredFrontierPlot<-function(object, 
 risk = c("Cov", "Sigma", "CVaR", "VaR"), 
 mText = NULL, 
 col = NULL, xlim = NULL, ylim = NULL, 
 twoAssets = FALSE) 
{
    offset = 0.1
    risk <- match.arg(risk)
    if (is.null(xlim)) {
        if (risk == "Cov") {
            xmax = max(sqrt(diag(getCov(object))))
        }
        if (risk == "Sigma") {
            xmax = max(sqrt(diag(getSigma(object))))
        }
        if (risk == "CVaR") {
            alpha = getAlpha(object)
            quantiles = colQuantiles(getSeries(object), prob = alpha)
            n.max = which.max(-quantiles)
            r = getSeries(object)[, n.max]
            r = r[r < quantiles[n.max]]
            xmax = -mean(r)
        }
        if (risk == "VaR") {
            xmax = max(-colQuantiles(getSeries(object), prob = alpha))
        }
        xlim = c(0, xmax)
        Xlim = c(xlim[1] - diff(xlim) * offset, xlim[2] + diff(xlim) * 
            offset)
    }
    if (is.null(ylim)) {
        ylim = range(getMean(object))
        Ylim = c(ylim[1] - diff(ylim) * offset, ylim[2] + diff(ylim) * 
            offset)
    }
    frontierPlot(object, pch = 19, risk = risk, xlim = Xlim, 
        ylim = Ylim)
    if (is.null(mText)) 
        mText = getTitle(object)
    mtext(mText, side = 3, line = 0.5, font = 2)
    grid()
    abline(h = 0, col = "grey")
    abline(v = 0, col = "grey")
    data = getData(object)
    spec = getSpec(object)
    constraints = getConstraints(object)
    mvPortfolio = minvariancePortfolio(data, spec, constraints)
    minvariancePoints(object, risk = risk, auto = FALSE, pch = 19, 
        col = "red")
    tangencyPoints(object, risk = risk, pch = 19, col = "blue")
    myTangencyLines(object, risk = risk, col = "blue")
    xy = equalWeightsPoints(object, risk = risk, pch = 15, col = "grey")
    text(xy[, 1] + diff(xlim)/20, xy[, 2] + diff(ylim)/20, "EWP", 
        font = 2, cex = 0.7)
    if (is.null(col)) 
        col = rainbow(6)
    xy = singleAssetPoints(object, risk = risk, cex = 1.5, col = col, 
        lwd = 2)
    text(xy[, 1] + diff(xlim)/20, xy[, 2] + diff(ylim)/20, rownames(xy), 
        font = 2, cex = 0.7)
    if (twoAssets) {
        twoAssetsLines(object, risk = risk, lty = 3, col = "grey")
    }
    invisible(object)
}

#this is where the 'action' begins

djiData<-readSeries(file="spdrAdj.csv",header=TRUE,sep=",")

by<-unique(timeLastNdayInMonth(time(djiData), 5))
djiData<-aggregate(djiData,by,mean)

djiData.ret<-100*returns(djiData,method="discrete")

bmp("returns.bmp")
plot(djiData.ret)
dev.off()

bmp("prices.bmp")
plot(djiData)
dev.off()

bmp("corrMatrix.bmp")
assetsCorImagePlot(djiData.ret)
dev.off()

globminSpec<-portfolioSpec()
globminPortfolio<-minvariancePortfolio(
 data=djiData.ret,
 spec=globminSpec,
 constraints="LongOnly")
print(globminPortfolio)

tgSpec<-portfolioSpec()
setRiskFreeRate(tgSpec)<-0.15
tgPortfolio<-tangencyPortfolio(
 data=djiData.ret,
 spec=tgSpec,
 constraints="LongOnly")
print(tgPortfolio)

pfront<-portfolioFrontier(
 data=djiData.ret,
 spec=tgSpec,
 constraints="LongOnly")

bmp("frontier.bmp")
myTailoredFrontierPlot(pfront)
dev.off()

The output includes the minimum variance portfolio:


Title:
 MV Minimum Variance Portfolio 
 Estimator:         covEstimator 
 Solver:            solveRquadprog 
 Optimize:          minRisk 
 Constraints:       LongOnly 

Portfolio Weights:
   XLY    XLP    XLE    XLF    XLV    XLI    XLB    XLK    XLU 
0.0000 0.6629 0.0683 0.0000 0.2367 0.0000 0.0000 0.0000 0.0321 

Covariance Risk Budgets:
   XLY    XLP    XLE    XLF    XLV    XLI    XLB    XLK    XLU 
0.0000 0.6629 0.0683 0.0000 0.2367 0.0000 0.0000 0.0000 0.0321 

Target Return and Risks:
  mean     mu    Cov  Sigma   CVaR    VaR 
0.1472 0.1472 3.0633 3.0633 7.9035 5.2951

The tangency portfolio:


Title:
 MV Tangency Portfolio 
 Estimator:         covEstimator 
 Solver:            solveRquadprog 
 Optimize:          minRisk 
 Constraints:       LongOnly 

Portfolio Weights:
   XLY    XLP    XLE    XLF    XLV    XLI    XLB    XLK    XLU 
0.0000 0.0000 0.9999 0.0000 0.0001 0.0000 0.0000 0.0000 0.0000 

Covariance Risk Budgets:
XLY XLP XLE XLF XLV XLI XLB XLK XLU 
  0   0   1   0   0   0   0   0   0 

Target Return and Risks:
   mean      mu     Cov   Sigma    CVaR     VaR 
 0.7914  0.7914  5.3748  5.3748 11.4222  7.4833

And some nice plots; monthly returns:

Monthly price series:

Correlation matrix of ETF monthly returns:

Efficient frontier with minimum variance and tangency portfolios:

The slope of the line that connects the risk-free asset (return of 0.15% on today's 4-week Treasury Bill) on the vertical axis is called the Sharpe ratio, and can be thought of as a reward-to-risk ratio (we will revisit this number in future posts).

A remarkable result of all this is that an optimal portfolio (in the sense that the Sharpe ratio is maximized) can be achieved by some mix of the risk-free asset with the tangency portfolio of risky assets.

Stock Returns

2009-07-06T08:10:00.000-07:00

A good (concise) reference for statistical models and methods for financial markets is this book, which goes by just that name. One of the authors teaches at Stanford, and the other one at Columbia University, which incidentally is where the winner from the last Interactive Brokers Collegiate Trading Olympiad is from. As I've mentioned before, I find the financial markets a particularly appealing environment for experimenting with autonomous agents, and so for the next few posts we will be reviewing how some of the machine learning methods described in CS229 can be used to analyze financial data.

The first example illustrates the use of linear regression to model the daily log returns of the stock of a given company as a function of the log returns of several other companies. The data and code can be downloaded from the book's website.

Why returns? Well, returns (given by the change in price of assets) are what dictate the profit (or loss) resulting from investing in stocks (neglecting transaction fees, taxes, etc), and so it makes sense to study their behavior. The big question is whether or not stock returns are predictable. The short answer is that the statistical evidence suggests they are not. The long answer is, of course, a little more complicated than that and is the subject of much debate and research (for a good, non-technical overview check Malkiel's "A Random Walk Down Wall Street."

Why log returns? Well, log returns (nearly equivalent to returns on a small scale) exhibit the additive property, which allows us to add single-period log returns to get the total multi-period log return (as opposed to multiplying the single-period returns). This additive property allows us to make use of the central limit theorem and say that the log returns (single and multi-period) are independent and identically, normally, distributed. Assuming that this holds, then the price movement could readily be described as a geometric random walk, which is the discrete-time equivalent of the geometric Brownian motion, a widely used gaussian process model of stock prices in option pricing.

At this point I would like to remind the reader of a quote attributed to George Box: "All models are wrong, but some are useful." Essentially, always keep in mind that models are (or try to be) useful representations of an observed process, they are not the process, and so there is pretty much always room for improvement. For example, later on we will look at time series models of stock returns that do not make the i.i.d. assumption.

To illustrate these ideas, we can simulate geometric random walks in R. We imagine the process to be described by a drift term and a volatility term. Drift accounts for the upward historical trend observed in stock prices (the only term needed for modeling returns in risk-free assets such as savings accounts or treasury bonds), where as volatility accounts for how much the price fluctuates form the path described by just the drift. Thus, a generative process would go something like this:

Start at some price Pinitial
Set Pnow=Pinitial
Draw a log return 'r' from a normal distribution with drift 'mu' and volatility 'sd'
Set Pnow = Pnow * exp(r)
Repeat from step 3 until done

In the example below, we examine an asset with annual mean (drift) 10% and sd (volatility) 20%, not unreasonable for large stocks in the U.S. for the last 80 years or so:


#grw_impl.r

#for repeatability
set.seed(123)

#simulation parameters
mu<-0.1
sd<-0.2
Pinitial<-100
N<-80
t<-1:N

#create placeholders for time series
price_series1<-1:N
price_series2<-1:N

#sample vector from normal distribution
sample_returns1<-rnorm(t,mu,sd)
sample_returns2<-rnorm(t,mu,sd)

#seed both series with same initial price
Pnow1<-Pinitial
Pnow2<-Pinitial

#run both series forward
for(i in 1:N)
{
  price_series1[i]<-Pnow1*exp(sample_returns1[i])
  Pnow1<-price_series1[i]

  price_series2[i]<-Pnow2*exp(sample_returns2[i])
  Pnow2<-price_series2[i]
}

#plot
bmp("price_realization_samples.bmp")

ymin<-min(Pinitial,price_series1,price_series2)
ymax<-max(Pinitial,price_series1,price_series2)
plot(c(0,t),c(Pinitial,price_series1),type='l',xlab="year",ylab="price",xlim=c(-10,N),ylim=c(ymin,ymax))
lines(c(0,t),c(Pinitial,price_series2),lty=2)
legend(-9.5,ymax,c("price_series1","price_series2"),lty=c(1,2),merge=TRUE,bg='gray90')

dev.off()

#pretend this is unknown time series and calculate drift and volatility
returns1<-rep(NA,N)
returns2<-rep(NA,N)

for(i in 2:N)
{
  returns1[i]<-log(1+(price_series1[i]-price_series1[i-1])/price_series1[i-1])
  returns2[i]<-log(1+(price_series2[i]-price_series2[i-1])/price_series2[i-1])
}

cat("\nsingle period drift 1:\t",mean(returns1,na.rm=TRUE),"\n")
cat("\nsingle period drift 2:\t",mean(returns2,na.rm=TRUE),"\n")

cat("\nsingle period volatility 1:\t",sd(returns1,na.rm=TRUE),"\n")
cat("\nsingle period volatility 2:\t",sd(returns2,na.rm=TRUE),"\n")

Driven with:


source("grw_impl.r")

Yielding:


single period drift 1:  0.1053144 

single period drift 2:  0.09006281 

single period volatility 1:  0.1854532 

single period volatility 2:  0.1906028

Notice that despite the 2 very different realization paths in the figure, the 2 series came from the same generating process and both display similar drift and volatility.

Now, could we ascertain from these 2 sample paths that the log returns are normally distributed? A qualitative way to determine this is to plot a histogram of the log returns and see what it looks like:


bmp("returnHist.bmp")
par(mfrow=c(2,1)) 
hist(returns1)  
hist(returns2)
dev.off()

Another qualitative way would be a normal probability plot:


par(mfrow=c(2,1)) 
qqnorm(returns1);qqline(returns1)
qqnorm(returns2);qqline(returns2)
dev.off()

A quantitative way of determining normality is to look at the sample skewness and kurtosis to see if those values near those of a normal distribution (zero for both):


library(moments)
cat("\nsample skewness 1:\t",skewness(returns1,na.rm=TRUE),"\n")
cat("\nsample skewness 2:\t",skewness(returns2,na.rm=TRUE),"\n")
cat("\nsample kurtosis 1:\t",kurtosis(returns1,na.rm=TRUE)-3,"\n")
cat("\nsample kurtosis 2:\t",kurtosis(returns2,na.rm=TRUE)-3,"\n")

From wikipedia: "The "minus 3" at the end of this formula is often explained as a correction to make the kurtosis of the normal distribution equal to zero." The excess skewness and kurtosis (in either direction) gives us a measure of how much the sample data deviates from a normal distribution:


sample skewness 1:  0.03947952 

sample skewness 2:  0.2006779 

sample kurtosis 1:  -0.1816579 

sample kurtosis 2:  -0.383789

Of course everything fits almost perfectly since we carefully generated the data (even if just little of it); real stock returns will most likely not look as pretty, but fairly close.

So, now that we've talked a little about log returns, we come back to the original example, which attempted to model the log returns of a specific company, MSFT, as a linear function of the log returns of several other similar companies (AAPL, DELL, YHOO, etc). The code is self explanatory, and the exercise is divided in 3 parts: regression for full model, variable selection, and regression diagnostics.

Regression for full model

R's 'lm' method is used to fit a linear model using all candidate companies in the data frame. Notable functions include 'plot', which, with the data frame as argument, plots the entire scatter plot matrix - also 'cor' - which outputs a matrix of correlations that is handy to have while looking at the scatter plot matrix.

Variable Selection

Up until now, we have been looking simply at the MSE for comparing regression methods; however statisticians have developed a whole slew of metrics that can be used to determine the quality of a model and guide the feature selection process. Calling 'summary' on the object returned by 'lm' provides the following:


Call:
lm(formula = msft ~ aapl + adbe + adp + amd + dell + gtw + hp + 
    ibm + orcl + sunw + yhoo)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0380497 -0.0028902 -0.0002195  0.0027663  0.0286453 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.782e-05  1.634e-04   0.170   0.8649    
aapl         4.180e-02  1.624e-02   2.573   0.0102 *  
adbe         7.821e-02  1.389e-02   5.630 2.22e-08 ***
adp          4.891e-02  2.266e-02   2.158   0.0311 *  
amd          1.920e-02  1.069e-02   1.796   0.0727 .  
dell         1.771e-01  2.127e-02   8.325  < 2e-16 ***
gtw          3.959e-02  9.919e-03   3.991 6.97e-05 ***
hp           2.703e-02  1.772e-02   1.526   0.1273    
ibm          2.302e-01  2.753e-02   8.360  < 2e-16 ***
orcl         1.264e-01  1.541e-02   8.204 5.73e-16 ***
sunw        -2.057e-02  1.213e-02  -1.696   0.0902 .  
yhoo         2.583e-02  1.238e-02   2.085   0.0372 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.005753 on 1243 degrees of freedom
Multiple R-squared: 0.5726, Adjusted R-squared: 0.5688 
F-statistic: 151.4 on 11 and 1243 DF,  p-value: < 2.2e-16

'Std. Error', a measure of precision, is the standard deviation of the estimate. 't-value' is the t-statistic for testing that the specified coefficient is zero; it is obtained by dividing the estimate by its standard error and it's used in the calculation of the p-value. 'Pr(>|t|)' is the p-value. If the calculated p-value is below the threshold chosen for statistical significance (usually the 0.10, the 0.05, or 0.01 level), then the null hypothesis (that the coefficient is zero) is rejected in favor of the alternative hypothesis. R gives you a nifty little reminder by adding a graphical representation of the p-value for each estimate: the more stars, the better. One could use this information for tuning the model, for example, by removing statistically insignificant regressors (i.e. HP, AMD, SUNW) in a backward elimination procedure (alternatively, R includes the 'stepAIC' function, which performs stepwise model selection by Akaike Information Criterion).

Regression Diagnostics

The techniques mentioned above to test for normality (normal probability plot, skewness and kurtosis, etc) are also useful here, since under the probabilistic view of least squares regression, we assume the errors (residuals) to be normally distributed.

Supervised Learning Method Comparison

2009-07-01T17:34:00.000-07:00

"Which supervised learning method works best for what?" - an empirical survey by Rich Caruana (Video & Slides).

Abstract:
"Decision trees are intelligible, but do they perform well enough that you should use them? Have SVMs replaced neural nets, or are neural nets still best for regression, and SVMs best for classification? Boosting maximizes margins similar to SVMs, but can boosting compete with SVMs? And if it does compete, is it better to boost weak models, as theory might suggest, or to boost stronger models? Bagging is simpler than boosting -- how well does bagging stack up against boosting? Breiman said Random Forests are better than bagging and as good as boosting. Was he right? And what about old friends like logistic regression, KNN, and naive bayes? Should they be relegated to the history books, or do they still fill important niches?
In this talk we compare the performance of ten supervised learning methods on nine criteria: Accuracy, F-score, Lift, Precision/Recall Break-Even Point, Area under the ROC, Average Precision, Squared Error, Cross-Entropy, and Probability Calibration. The results show that no one learning method does it all, but some methods can be "repaired" so that they do very well across all performance metrics. In particular, we show how to obtain the best probabilities from max margin methods such as SVMs and boosting via Platt's Method and isotonic regression. We then describe a new ensemble method that combines select models from these ten learning methods to yield much better performance. Although these ensembles perform extremely well, they are too complex for many applications. We'll describe what we're doing to try to fix that. Finally, if time permits, we'll discuss how the nine performance metrics relate to each other, and which of them you probably should (or shouldn't) use."

Strategy Optimization

2009-06-29T04:14:00.000-07:00

A particularly thought-provoking chapter in AIMA is chapter 17, "Making Complex Decisions," which elegantly introduces many powerful concepts in a way that highlights their interrelations. The list below is an excerpt from the TOC:

Sequential Decision Problems
Partially observable MDPs
Decision-Theoretic Agents
Decisions with Multiple Agents: Game Theory
Mechanism Design

Of particular interest to me is how these topics relate to strategy discovery and optimization.

A strategy is a plan of action designed to achieve a particular goal, and strategy discovery can be thought of as a subset of strategy optimization since in theory, we could start from a trivial, arbitrary, possibly bad strategy and somehow optimize it until we have "discovered" a satisfactory strategy. Perhaps "strategy design" would more generally describe what I'm talking about, but I prefer to term "optimization" for all the sentiments it evokes.

We have already seen how we can use reinforcement learning methods to come up with satisfactory strategies for controllers in settings such as the inverted pendulum problem. In that context, we called the strategy a policy and defined the optimal policy to be a mapping from states to actions that would maximize the controller's expected return from the environment. What other kinds of scenarios may be these concepts and methods be useful for? Consider the following 2 examples:

Legend has it that Nathan Rothschild, a London Banker in the 19th century, acted cleverly on early news of the outcome of the battle of Waterloo and made a fortune. Consols would rise as a result of Napoleon's defeat, but rather than buying Consols anticipating the increase, Nathan sold his existing holdings. Other investors, unaware of the outcome in Waterloo, got wind of this and assumed Nathan was acting on information that the battle was lost. Nathan was known at the time to have many connections and was expected to have access to early information, so these other investors also sold their holdings, driving the prices further down. At the proverbial last minute, Nathan bought at bargain prices. When news of the battle came out, the prices soared, and Nathan realized a nice profit.
A more modern story is that of how Porsche cornered the market in VW stock and made a killing in 2008. You can read about it here, here, and here.

These and many other accounts, in finance and in other fields, are examples of successful strategies that while on one hand are amazing because of their ingenuity, on the other hand appear evident in hindsight, suggesting that recognizable causal mechanisms were exploited rather than being successful because of good fortune. Furthermore, after reading several such accounts, there appear to be common traces among such successful strategies (even across disparate fields) which suggests that there may exist a template strategy (or a reusable process for arriving at such strategy) that may be customized to specific scenarios, leading to the question: "could such strategy optimization processes be automated?"

The strategy optimization process can indeed be automated, and a framework for doing so draws from fields such as decision theory, game theory, system dynamics, and control theory, among others, making heavy use of modeling and simulation, not unlike our inverted pendulum example. Take the London Banker, for instance, who if ambivalent about the outcomes resulting from initially selling or buying, could run simulations of each strategy to determine the most profitable one. Of course, the success of this exercise would primarily depend on the quality of the models for each independent system component and their interrelations, some of which may include a model for information flow, a model for Nathan and traders other than Nathan, and a model of the "physics" of the environment (i.e. what actions may or may not be taken by market participants at any given time).

And therein lies the rub: it is very hard to model complex systems, specially when it involves modeling human behavior. Nevertheless, there exist methods to ease the burden and that is something I expect to dive into deeper in future posts.

Lec17 - Fitted Value Iteration

2009-05-31T10:50:00.000-07:00

Video Lecture - 17

Below is a python implementation for controlling an inverted pendulum in a cart (described here and animated here) using the so called bang-bang control regime (alternatively, we could have used some measure of proportional control), in which the state space is discretized. The task is further defined in problem set #4 (whose solution includes a matlab implementation) and uses system dynamics defined here:


#polecart_disc.py
from pylab import *
from numpy import *

#Displays the pole-cart graph.
class show_cart_params: pass
show_cart_args=show_cart_params()
show_cart_args.first=True
show_cart_args.trial=0
def show_cart(state_vector):
    figure(1)
    length=3
    x=state_vector[0,0];x_dot=state_vector[1,0]
    theta=state_vector[2,0];theta_dot=state_vector[3,0]
    x1=x
    y1=0
    x2=x+length*sin(theta)
    y2=length*cos(theta)
    cart_linex1,cart_liney1=[x-.4,x+.4],[-.25,-.25]
    cart_linex2,cart_liney2=[x+.4,x+.4],[-.25,0]
    cart_linex3,cart_liney3=[x+.4,x-.4],[0,0]
    cart_linex4,cart_liney4=[x-.4,x-.4],[0,-.25]
    base_linex,base_liney=[x,x],[-.5,-.25]
    axis([-3,3,-0.5,3.5])
    global show_cart_args
    
    if show_cart_args.first:
        show_cart_args.first=False 
        show_cart_args.pend_line,=plot([x1,x2],[y1,y2],c='black')
        show_cart_args.cart_line1,=plot(cart_linex1,cart_liney1,c='cyan')
        show_cart_args.cart_line2,=plot(cart_linex2,cart_liney2,c='cyan')
        show_cart_args.cart_line3,=plot(cart_linex3,cart_liney3,c='cyan')
        show_cart_args.cart_line4,=plot(cart_linex4,cart_liney4,c='cyan')
        show_cart_args.base_line,=plot(base_linex,base_liney,c='cyan')

    show_cart_args.pend_line.set_data([x1,x2],[y1,y2])
    show_cart_args.cart_line1.set_data(cart_linex1,cart_liney1)
    show_cart_args.cart_line2.set_data(cart_linex2,cart_liney2)
    show_cart_args.cart_line3.set_data(cart_linex3,cart_liney3)
    show_cart_args.cart_line4.set_data(cart_linex4,cart_liney4)
    show_cart_args.base_line.set_data(base_linex,base_liney)
    if show_cart_args.trial!=num_failures:
        title('trial #%d'%num_failures)
        show_cart_args.trial=num_failures
    draw()

#Returns a discretized value (an index) for a continuous
#state vector: x is divided into 3 "boxes", x_dot into 3,
#theta into 6, and theta_dot into 3.
def get_state(state_vector):
    one_degree=0.0174532
    six_degrees=0.1047192
    twelve_degrees=0.2094384
    fifty_degrees=0.87266
    total_states=163
    state=0
    x=state_vector[0,0];x_dot=state_vector[1,0]
    theta=state_vector[2,0];theta_dot=state_vector[3,0]

    if -2.4>x or x>2.4 or -twelve_degrees>theta or theta>twelve_degrees:
        #went past boundaries
        state=total_states-1
    else:
        if -1.5>x:
            state=0
        elif x<1.5:
            state=1
        else:
            state=2

        if -.5>x_dot:
            pass
        elif .5>x_dot:
            state+=3
        else:
            state+=6

        if -six_degrees>theta:
            pass
        elif -one_degree>theta:
            state+=9
        elif 0>theta:
            state+=18
        elif one_degree>theta:
            state+=27
        elif six_degrees>theta:
            state+=36
        else:
            state+=45

        if -fifty_degrees>theta_dot:
            pass
        elif fifty_degrees>theta_dot:
            state+=54
        else:
            state+=108

    return state

#Simulates pole-cart dynamics.
def cart_pole(action,state_vector):
    x=state_vector[0,0];x_dot=state_vector[1,0]
    theta=state_vector[2,0];theta_dot=state_vector[3,0]
    
    gravity=9.8
    masscart=1.0
    masspole=.3
    total_mass=masspole+masscart
    length=.7
    polemass_length=masspole*length
    force_mag=10.0
    tau=.02
    fourthirds=1.3333333333333

    action_flip_prob=0.0
    force_noise_factor=0.0
    no_control_prob=0.0

    if action_flip_prob>random.random():
        action=1-action

    if action>0:
        force=force_mag
    else:
        force=-force_mag

    force=force* \
           (1-force_noise_factor+random.random()*2*force_noise_factor)

    if no_control_prob>random.random():
        force=0

    costheta=cos(theta)
    sintheta=sin(theta)

    temp=(force+polemass_length*theta_dot*theta_dot*sintheta)/total_mass
    thetaacc=(gravity*sintheta-costheta*temp)/(length*(fourthirds \
            -masspole*costheta*costheta/total_mass))
    xacc=temp-polemass_length*thetaacc*costheta/total_mass

    new_x=x+tau*x_dot
    new_x_dot=x_dot+tau*xacc
    new_theta=theta+tau*theta_dot
    new_theta_dot=theta_dot+tau*thetaacc

    return array([new_x,new_x_dot,new_theta,new_theta_dot]).reshape(4,1)

ion()
display_started=False
num_states=163
gamma=.995
tolerance=.01
no_learning_threshold=20
time=0
time_steps_to_failure={}
num_failures=0
time_at_start_of_current_trial=0
transition_counts=zeros((2,num_states,num_states))
transition_probs=ones((2,num_states,num_states))/num_states
reward_counts=zeros((num_states,2))
reward=zeros((num_states,1))
value=random.ranf((num_states,1))*0.1
consecutive_no_learning_trials=0
x=0.0;x_dot=0.0;theta=0.0;theta_dot=0.0
state_vector=array([x,x_dot,theta,theta_dot]).reshape(4,1)
state_index=get_state(state_vector)

#The main loop in this animation, will continue to run
#until the number of consecutive trials in which value
#iteration converged within one iteration is whatever
#no_learning_threshold is defined to be.
while no_learning_threshold>consecutive_no_learning_trials:
    
    #Calculate the expected value of taking each of the
    #2 actions we have at our disposal. If one action is
    #not clearly better than the other, choose randomly.
    score1=dot(transition_probs[0,state_index,:],value)
    score2=dot(transition_probs[1,state_index,:],value)
    if score1>score2:
        action=0
    elif score2>score1:
        action=1
    else:
        if random.random()>.5:
            action=1
        else:
            action=0

    #Take the action and see what state we end up in next.
    state_vector=cart_pole(action,state_vector)
    time+=1
    new_state_index=get_state(state_vector)

    if display_started:
        show_cart(state_vector)

    #This here defines the reward function for this task: -1
    #if the pole falls or the cart has gone beyond the specified
    #boundary, and 0 otherwise.  The last state is considered
    #to be the one that satisfies these conditions.
    if new_state_index==num_states-1:
        R=-1.0
    else:
        R=0.0

    #Update transition counts and reward counts.  We keep track of
    #2 things in reward_counts: the number of times we received a
    #reward in that state and the total reward we received in that state
    transition_counts[action,state_index,new_state_index]= \
        transition_counts[action,state_index,new_state_index]+1.0
    reward_counts[new_state_index,1]=reward_counts[new_state_index,1]+1.0
    reward_counts[new_state_index,0]=reward_counts[new_state_index,0]+R

    if new_state_index==num_states-1:

        #If we ended up in the terminal state, this trial is over,
        #so let's now update our statistics.
        #First, turn the transition counts into probabilities,
        for a in range(2):
            for s in range(num_states):
                den=sum(transition_counts[a,s,:])
                if den>0:
                    transition_probs[a,s,:]=transition_counts[a,s,:]/den

        #then calculate the expected reward in each state using
        #the total reward for that state and how many times
        #we received a reward in that state.
        for s in range(num_states):
            if reward_counts[s,1]>0:
                reward[s,0]=reward_counts[s,0]/reward_counts[s,1]

        iterations=0
        new_value=zeros((num_states,1))
        while True:
            #Now update value estimates by performing
            #value iteration until convergence, which is defined
            #as our value estimates not changing below the
            #specified tolerance.
            iterations=iterations+1
            for s in range(num_states):
                value1=dot(transition_probs[0,s,:],value)
                value2=dot(transition_probs[1,s,:],value)
                new_value[s,0]=max(value1,value2)
            new_value=reward+gamma*new_value
            diff=max(abs(value-new_value))
            value=new_value+0#copy
            if tolerance>diff:
                break

        #if we converged in 1 iteration, keep track of that because
        #once this happens no_learning_threshold times in a row,
        #our experiment is over.
        if iterations==1:
            consecutive_no_learning_trials+= 1
            #watch the last trial
            if consecutive_no_learning_trials>no_learning_threshold-2:
                display_started=True
        else:
            consecutive_no_learning_trials=0
            display_started=False

        #num_failures can be thought of as the trial number
        num_failures+=1
        #update how many time steps we lasted in this trial
        time_steps_to_failure[num_failures]= \
            time-time_at_start_of_current_trial
        time_at_start_of_current_trial=time
        #reset our state vector to some random location
        x=random.uniform(-1.1,1.1);x_dot=0.0;theta=0.0;theta_dot=0.0
        state_vector=array([x,x_dot,theta,theta_dot]).reshape(4,1)
        state_index=get_state(state_vector)
    else:
        state_index=new_state_index

#Plot learning curve:
figure(2)
plot(time_steps_to_failure.values(),'k')
title('Time steps per trial')
ioff()
show()

The output includes an animated plot of the cart on its last trial and a plot showing how many time-steps the system was balanced within constraints on each trial:

Rather than discretizing the state space, we could directly estimate the value as a function of a continuous feature vector which is derived from the "perceived" continuous state vector. This approach, fitted value iteration, is described here, and is just one of the many methods applied to the inverted pendulum problem:


#polecart_cont.py
import sys
from pylab import *
from numpy import *

#To add more artificial features just append them to feats.
def get_feat_vec(x,x_dot,theta,theta_dot):
    feats=[1.0,
           x,
           x_dot,
           theta,
           theta_dot,
           x_dot**2,
           theta_dot**2,
           theta*x]
    return array(feats).reshape(1,len(feats))

#Displays the cart pole graph.
class show_cart_params: pass
show_cart_args=show_cart_params()
show_cart_args.first=True
show_cart_args.trial=0
def show_cart(x,x_dot,theta,theta_dot):
    figure(1)
    length=3
    x1=x
    y1=0
    x2=x+length*sin(theta)
    y2=length*cos(theta)
    cart_linex1,cart_liney1=[x-.4,x+.4],[-.25,-.25]
    cart_linex2,cart_liney2=[x+.4,x+.4],[-.25,0]
    cart_linex3,cart_liney3=[x+.4,x-.4],[0,0]
    cart_linex4,cart_liney4=[x-.4,x-.4],[0,-.25]
    base_linex,base_liney=[x,x],[-.5,-.25]
    axis([-3,3,-0.5,3.5])
    global show_cart_args
    
    if show_cart_args.first:
        show_cart_args.first=False 
        show_cart_args.pend_line,=plot([x1,x2],[y1,y2],c='black')
        show_cart_args.cart_line1,=plot(cart_linex1,cart_liney1,c='cyan')
        show_cart_args.cart_line2,=plot(cart_linex2,cart_liney2,c='cyan')
        show_cart_args.cart_line3,=plot(cart_linex3,cart_liney3,c='cyan')
        show_cart_args.cart_line4,=plot(cart_linex4,cart_liney4,c='cyan')
        show_cart_args.base_line,=plot(base_linex,base_liney,c='cyan')

    show_cart_args.pend_line.set_data([x1,x2],[y1,y2])
    show_cart_args.cart_line1.set_data(cart_linex1,cart_liney1)
    show_cart_args.cart_line2.set_data(cart_linex2,cart_liney2)
    show_cart_args.cart_line3.set_data(cart_linex3,cart_liney3)
    show_cart_args.cart_line4.set_data(cart_linex4,cart_liney4)
    show_cart_args.base_line.set_data(base_linex,base_liney)
    if show_cart_args.trial!=num_failures:
        title('trial #%d'%num_failures)
        show_cart_args.trial=num_failures
    draw()

#Returns whether or not the system is out of bounds.
def is_terminal(x,theta):
    return -2.4>x or x>2.4 or -twelve_degrees>theta or theta>twelve_degrees

#Simulates cart pole dynamics.
def cart_pole(action,x,x_dot,theta,theta_dot):   
    gravity=9.8
    masscart=1.0
    masspole=.3
    total_mass=masspole+masscart
    length=.7
    polemass_length=masspole*length
    force_mag=10.0
    tau=.02
    fourthirds=1.3333333333333

    action_flip_prob=0.0
    force_noise_factor=0.0
    no_control_prob=0.0

    if action_flip_prob>random.random():
        action=1-action

    if action>0:
        force=force_mag
    else:
        force=-force_mag

    force=force* \
           (1-force_noise_factor+random.random()*2*force_noise_factor)

    if no_control_prob>random.random():
        force=0

    costheta=cos(theta)
    sintheta=sin(theta)

    temp=(force+polemass_length*theta_dot*theta_dot*sintheta)/total_mass
    thetaacc=(gravity*sintheta-costheta*temp)/(length*(fourthirds \
            -masspole*costheta*costheta/total_mass))
    xacc=temp-polemass_length*thetaacc*costheta/total_mass

    new_x=x+tau*x_dot
    new_x_dot=x_dot+tau*xacc
    new_theta=theta+tau*theta_dot
    new_theta_dot=theta_dot+tau*thetaacc

    return get_feat_vec(new_x,new_x_dot,new_theta,new_theta_dot)

gamma=0.995
feat_vec=get_feat_vec(0,0,0,0)
n=shape(feat_vec)[1]
twelve_degrees=0.2094384
m=100000
w_vec=zeros((n,1))
feat_mat=ones((m,n))
diff=sys.maxint
new_diff=diff-1
epochs=0

#generate a bunch of state samples to fit the value function
x=random.uniform(-1.1,1.1);x_dot=0.0;theta=0.0;theta_dot=0.0
for i in range(m):
    if random.random()>.5:
        action=1
    else:
        action=0    
        
    feat_vec=cart_pole(action,x,x_dot,theta,theta_dot)
    feat_mat[i,:]=feat_vec

    if is_terminal(feat_vec[0,1],feat_vec[0,3]):
        x=random.uniform(-1.1,1.1);x_dot=0.0
        theta=random.uniform(-twelve_degrees,twelve_degrees);theta_dot=0.0
    else:
        x=feat_vec[0,1];x_dot=feat_vec[0,2]
        theta=feat_vec[0,3];theta_dot=feat_vec[0,4]

#Iteratively score each of the generated states and use the
#resulting feature,target tuple to estimate the parameters
#of a linear function.
while 10>epochs and diff>new_diff:
    diff=new_diff
    epochs+=1

    y_vec=zeros((m,1))   
    for i in range(m):
        feat_vec=feat_mat[i,:].reshape(1,n)
        x=feat_vec[0,1];x_dot=feat_vec[0,2]
        theta=feat_vec[0,3];theta_dot=feat_vec[0,4]

        if is_terminal(x,theta):
            R=-1.0
        else:
            R=0.0         
        
        feat_vec0=cart_pole(0,x,x_dot,theta,theta_dot)
        score1=R+gamma*dot(feat_vec0,w_vec)
        
        feat_vec1=cart_pole(1,x,x_dot,theta,theta_dot)
        score2=R+gamma*dot(feat_vec1,w_vec)
        
        if score1>score2:
            action=0
            score=score1
        elif score2>score1:
            action=1
            score=score2
        else:
            if random.random()>.5:
                action=1
                score=score2
            else:
                action=0
                score=score1

        y_vec[i,0]=score

    result=linalg.lstsq(feat_mat,y_vec)
    w_new=result[0]
    new_diff=result[1]
    w_vec=w_new+0#copy


ion()
time=0
time_steps_to_failure={}
num_failures=0
time_at_start_of_current_trial=0
x=random.uniform(-1.1,1.1);x_dot=0.0;theta=0.0;theta_dot=0.0
display_started=False
max_num_failures_to_run=250

#Run the learned controller a few times to plot and compare
#against the discretized version.
while max_num_failures_to_run>num_failures:
    
    feat_vec0=cart_pole(0,x,x_dot,theta,theta_dot)
    score1=dot(feat_vec0,w_vec)
    
    feat_vec1=cart_pole(1,x,x_dot,theta,theta_dot)
    score2=dot(feat_vec1,w_vec)
    
    if score1>score2:
        action=0
    elif score2>score1:
        action=1
    else:
        if random.random()>.5:
            action=1
        else:
            action=0

    feat_vec=cart_pole(action,x,x_dot,theta,theta_dot)
    x=feat_vec[0,1];x_dot=feat_vec[0,2]
    theta=feat_vec[0,3];theta_dot=feat_vec[0,4]
    time+=1

    if display_started:
        show_cart(x,x_dot,theta,theta_dot)

    if is_terminal(x,theta):
        num_failures+=1
        if num_failures>max_num_failures_to_run-2:
            display_started=True

        time_steps_to_failure[num_failures]= \
            time-time_at_start_of_current_trial
        time_at_start_of_current_trial=time

        x=random.uniform(-1.1,1.1);x_dot=0.0;theta=0.0;theta_dot=0.0

#Plot learning curve:
figure(2)
plot(time_steps_to_failure.values(),'k')
title('Time steps per trial')
ioff()
show()

Notice that this controller lasts much longer (although obviously there are nonlinearities not captured by our value function). It also performs the same on all trials, which has to do with how we are initializing the system (not completely random, but within an 'x' range that it can stay within, even though it allows the pole to fall below the specified constraint):

Pacman and Reinforcement Learning

2009-05-27T13:05:00.000-07:00

UC Berkeley's approach to teaching their undergraduate Artificial Intelligence class should be an example to all institutions out there offering a similar course. Not only have they not moved away from AIMA, they've designed the course with a strong focus on programming projects, giving students a chance to truly internalize the concepts learned. This year's class used the game Pacman as a running example to build intelligent agents driven by strategies of complexity ranging from search algorithms to probabilistic models, and it's so well done I would not be surprised if this turned out to be one of the most popular compsci undergraduate classes at Berkeley. Now, putting something together like this obviously requires a great deal of effort, but John DeNero, the course designer and instructor, actually encourages instructors at other universities to get in touch with him and use the materials for their own classes. So, the road is already paved; there should be no reason for other institutions not to go down it (University of Utah is already doing it).

Some of the code for the Pacman game is based on a implementation from LiveWires, and you can download their stuff for further documentation on the graphics, layout, and general spirit of the animation; however, this is not necessary for completing the programming projects.

Since we're talking about reinforcement learning, I thought it might nice to fill in the stubs for their RL project and experiment a little.

The first thing to do is to implement an agent that uses value iteration to compute the optimal policy on the grid world, and it does not take much to adapt the implementation in our last post to work in this scenario. Here is what the __init__ method in ValueIterationAgent would look like:


  def __init__(self, mdp, discount = 0.9, iterations = 100):
    """
      Your value iteration agent should take an mdp on
      construction, run the indicated number of iterations
      and then act according to the resulting policy.
    
      Some useful mdp methods you will use:
          mdp.getStates()
          mdp.getPossibleActions(state)
          mdp.getTransitionStatesAndProbs(state, action)
          mdp.getReward(state, action, nextState)
    """
    self.mdp = mdp
    self.discount = discount
    self.iterations = iterations
    self.values = util.Counter() # A Counter is a dict with default 0
     
    "*** YOUR CODE HERE ***"   
    #iteratively estimate opt value function
    values=dict([(s,0.0) for s in mdp.getStates()])
    for i in range(iterations):
      V=values.copy()
      for s in mdp.getStates():
        if not mdp.isTerminal(s):  
          now=mdp.getReward(s,None,None)
          actions=mdp.getPossibleActions(s)
          futures=[0.0]*len(actions)
          for a in actions:
            futures.append(sum([p*V[next] for (next,p) in
                                mdp.getTransitionStatesAndProbs(s,a)]))
          maxfuture=max(futures)
          values[s]=now+discount*maxfuture
    for k,v in values.items():
      self.values[k]=v
        
    #opt policy is greedy wrt opt value function
    pi={}
    for s in mdp.getStates():
      if not mdp.isTerminal(s):
        maxexpect=-sys.maxint
        for a in mdp.getPossibleActions(s):
          expect=sum([p*self.values[next] for (next,p) in
                      mdp.getTransitionStatesAndProbs(s,a)])
          if expect>maxexpect:
            maxexpect=expect
            pi[s]=a
      else:
        pi[s]=()
    self.pi=pi

For which the command below yields:


python gridworld.py -a value -i 100 -k 10

But the real fun begins when we move on to designing agents that learn from interaction with the environment rather than from a model that is just handed to them. Q-learning is the particular approach used for the rest of the exercise, and the 2 most salient methods in the QLearningAgent are update and getAction (which should really just call getPolicy, but no biggie):


  def getAction(self, s):
    """
      What action to take in the current state. With
      probability self.epsilon, we should take a random
      action and take the best policy action otherwise.
    
      After you choose an action make sure to
      inform your parent self.doAction(state,action) 
      This is done for you, just don't clobber it
       
      HINT: You might want to use util.flipCoin(prob)
      HINT: To pick randomly from a list, use random.choice(list)
    """  
    # Pick Action
    legalActions = self.getLegalActions(s)
    action = None
    "*** YOUR CODE HERE ***"
    if not util.flipCoin(self.epsilon):    
      maxactions=[]
      maxQ=max([self.getQValue(s,a) for a in self.getLegalActions(s)])
      for a in legalActions:
        if self.getQValue(s,a)>=maxQ:
          maxactions.append(a)
          maxQ=self.getQValue(s,a)
    else:
      maxactions=legalActions
        
    action=random.choice(maxactions)
    
    # Need to inform parent of action for Pacman (do not delete this line)
    self.doAction(s,action)    
    
    return action
  
  def update(self, s, a, nexts, r):
    """
      The parent class calls this to observe a 
      state = action => nextState and reward transition.
      You should do your Q-Value update here
      
      NOTE: You should never call this function,
      it will be called on your behalf
    """
    "*** YOUR CODE HERE ***"
    if len(self.getLegalActions(nexts))==0:
      self.Q[s,a]=r
    else:
      self.Q[s,a]=self.getQValue(s,a)+ \
                   self.alpha* \
                   (r+ \
                    self.gamma* \
                    max([self.getQValue(nexts,nexta) for nexta in self.getLegalActions(nexts)])- \
                    self.getQValue(s,a))

Watching this code run on the grid world is really amazing for 2 reasons, first because it works :) but more importantly because of how well it works given that there is no explicit model of the environment. This is a little bit like neural nets, whose "knowledge" is implicit in the weights assigned to connections between nodes, well similarly, a q-learning agent's "knowledge" is implicit in its Q function, and the dynamics of the environment (provided it is well explored) are reflected in the values assigned to each state-action pair. Sutton and Barto's book has a good discussion on this towards the end of part II and most of part III.

After the grid world, you get a chance to play with the crawler, which is an animation of a one-legged robot that drags itself around your screen, good to get an intuition on several of the algorithm parameters, especially epsilon, which helps balance the agent's dedication to exploration vs exploitation:

It's easy to start thinking about all sorts of things as you control the knobs and watch this little thing crawling on your screen: can we design an ultimate reward signal (or even a few) such as survival or continuation of the species? The robot could then be looking for food or a mate, and if it did mate, could its offspring "perform" better by moving itself faster, perhaps combining the learned strategies of its ancestors (think genetic algorithms)? Or the knobs that we control (epsilon, learning rate, discount), do they have equivalents in the brain, some neurotransmitter or other chemicals perhaps? How could a model of the environment speed learning up? What kind of model? I'm sure most people playing around with this game would ask themselves these questions (and many more) and that is one of the objectives of a class like this: to stimulate the intellect to all that is possible (and as time goes by, seemingly more and more within grasp) and give students the tools to answer those questions; an objective that is more than fulfilled in this case. Plus, we aren't even done yet, there is still the matter of equipping Pacman with q-learning abilities.

Once you get to the ApproximateQAgent, the most important methods would be getQValue and update:


  def getQValue(self, s, a):
    """
      Should return Q(state,action) = w * featureVector
      where * is the dotProduct operator
    """
    "*** YOUR CODE HERE ***"
    #return PacmanQAgent.getQValue(self, s, a)
    feats=self.featExtractor.getFeatures(s,a)
    return sum([feats[k]*self.W[k] for k,v in feats.items()])
    
  def update(self, s, a, nexts, r):
    """
       Should update your weights based on transition  
    """
    "*** YOUR CODE HERE ***"
    #PacmanQAgent.update(self, s, a, nexts, r)
    feats=self.featExtractor.getFeatures(s,a)
    if len(self.getLegalActions(nexts))==0:
      for k,v in feats.items():
        self.W[k]=r*feats[k]
    else:
      for k,v in feats.items():
        corr=r+self.gamma*self.getValue(nexts)- \
              self.getQValue(s,a)
        self.W[k]=self.W[k]+ \
                   self.alpha*corr*feats[k]

And that should be enough to run on some of the bigger grids with a decent chance of winning:


python pacman.py -p ApproximateQAgent -a extractor=SimpleExtractor -x 50 -n 60 -l mediumClassic 
Beginning 50 episodes of Training
Training Done (turning off epsilon and alpha)
---------------------------------------------
Pacman died! Score: -100
Pacman emerges victorious! Score: 1327
Pacman died! Score: -70
Pacman emerges victorious! Score: 1346
Pacman died! Score: -274
Pacman emerges victorious! Score: 591
Pacman died! Score: -180
Pacman emerges victorious! Score: 701
Pacman died! Score: 87
Pacman emerges victorious! Score: 1335
Average Score: 476.3
Scores:        -100, 1327, -70, 1346, -274, 591, -180, 701, 87, 1335
Win Rate:      5/10 (0.50)
Record:        Loss, Win, Loss, Win, Loss, Win, Loss, Win, Loss, Win

The last part of the project encourages you to build better feature extractors, which if done well should considerably improve Pacman's performance. This issue of generalizing in the state space is a recurring one not only in reinforcement learning, but throughout AI, and we've seen a few techniques as we learned about supervised and unsupervised learning. Next time we'll go back to cs229 and look at strategies for generalizing in continuous state-spaces.

Lec16 - Value and Policy Iteration

2009-05-25T08:04:00.000-07:00

Both value and policy iteration are algorithms for solving an MDP formulation of a reinforcement learning problem where the environment dynamics are completely known. Resources for this kind of scenario include:

chapter 17 of AIMA
the dynamic programming section of the RL book
the RL notes from cs229

Below we have some simple implementations of such algorithms running against a simple grid world described in AIMA:


#mdp_imp.py
import sys,operator,random

class Env:
    def __init__(self,states,rewards,terminals,dynamics,actions):
        self.states=states
        self.rewards=rewards
        self.terminals=terminals
        self.actions=actions
        self.dynamics=dynamics
        self.R=dict(zip(states,rewards))
        self.allowablestates=[states[i]
                              for i in range(len(states))
                              if rewards[i] is not None]

    def T(self,state,action):
        trans=[]
        for prob,move in self.dynamics:
            trans.append((prob,self.go(state,move(self,action))))
        return trans

    def go(self,state,action):
        next=tuple(map(operator.add,state,action))
        if next in self.allowablestates:
            return next
        else:
            return state

    def desired(self,action):
        return action

    def left(self,action):
        idx=self.actions.index(action)-1
        return self.actions[idx]

    def right(self,action):
        idx=(self.actions.index(action)+1)%len(self.actions)
        return self.actions[idx]
        
class Agent:
    def __init__(self,actions,gamma):
        self.actions=actions
        self.gamma=gamma

def display(agt,env,pi,desc):
    dirs={(1,0):'^',(0,1):'>',(-1,0):'v',(0,-1):'<',(0,0):'.'}
    print '\n'+desc
    for state in env.allowablestates:
        print 'state: %s, move: %s, action: %s'% \
              (state,dirs[pi[state]],pi[state])

def policyIter(agt,env):
    #create some random policy
    pi={}
    for s in env.terminals: pi[s]=(0,0)
    for s in set(env.allowablestates)-set(env.terminals):
        pi[s]=agt.actions[random.randint(0,len(agt.actions)-1)]

    #iteratively improve policy
    done=False
    while not done:
        done=True
        
        #first compute value function for current policy using modified
        #version of policy iteration (policy and iters are constant)
        V=dict([(s,0.0) for s in env.states])
        for s in env.terminals: V[s]=env.R[s]
        iters=20
        for i in range(iters):
            for s in set(env.allowablestates)-set(env.terminals):
                now=env.R[s]
                future=sum([p*V[next] for (p,next) in env.T(s,pi[s])])                
                V[s]=now+agt.gamma*future
        
        #then make policy greedy towards computed value function
        for s in set(env.allowablestates)-set(env.terminals):
            maxexpect=sum([p*V[next] for (p,next) in env.T(s,pi[s])])
            for a in agt.actions:
                expect=sum([p*V[next] for (p,next) in env.T(s,a)])
                if expect>maxexpect and a!=pi[s]:
                    done=False
                    maxexpect=expect
                    pi[s]=a
    return pi

def valueIter(agt,env):
    #initialize value function for each state
    V=dict([(s,0.0) for s in env.states])
    for s in env.terminals: V[s]=env.R[s]

    #iteratively estimate opt value function
    delta=sys.maxint
    while (delta>1e-5):
        delta=-sys.maxint
        for s in set(env.allowablestates)-set(env.terminals):
            now=env.R[s]
            maxfuture=max([sum([p*V[next] for (p,next) in env.T(s,a)])
                           for a in agt.actions])
            value=now+agt.gamma*maxfuture
            delta=max(delta,abs(value-V[s]))
            V[s]=value
            
    #opt policy is greedy wrt opt value function
    pi={}
    for s in env.terminals: pi[s]=(0,0)
    for s in set(env.allowablestates)-set(env.terminals):
        maxexpect=-sys.maxint
        for a in agt.actions:
            expect=sum([p*V[next] for (p,next) in env.T(s,a)])
            if expect>maxexpect:
                maxexpect=expect
                pi[s]=a
    return pi

Driven with:


#mdp_inter.py
from mdp_imp import *

#the states in the environment and
#their (x,y) coordinates
states=[(2,0),(2,1),(2,2),(2,3),
        (1,0),(1,1),(1,2),(1,3),
        (0,0),(0,1),(0,2),(0,3)]

#rewards obtained for simply entering
#each state. 'None' denotes an unreachable state
rewards=[-0.04,-0.04,-0.04,+1.00,
         -0.04, None,-0.04,-1.00,
         -0.04,-0.04,-0.04,-0.04]

#states considered terminal (absorbing)
terminals=[(2,3),(1,3)]

#given an action, the environment dynamics specify that
#with .8 probability it will have the intended effect (move
#the agent in the desired direction), with .1 probability
#the agent will go left, and with .1 probability the
#agent will go right. If the action takes the agent
#outside the boundaries of the environment or into a 'None'
#state, then the agent will remain where it is
dynamics=[(0.8,Env.desired),(0.1,Env.left),(0.1,Env.right)]

#the only actions the agent can take are to move
#north,east,south,west
#denoted below with unit vectors (y,x)
actions=[(1,0),(0,1),(-1,0),(0,-1)]

#gamma denotes the discount factor, indicative of an
#agent's preference of present rewards over future rewards
gamma=.9

agt=Agent(actions,gamma)
env=Env(states,rewards,terminals,dynamics,actions)
display(agt,env,valueIter(agt,env),'Value Iteration')
display(agt,env,policyIter(agt,env),'Policy Iteration')

With output:


Value Iteration
state: (2, 0), move: >, action: (0, 1)
state: (2, 1), move: >, action: (0, 1)
state: (2, 2), move: >, action: (0, 1)
state: (2, 3), move: ., action: (0, 0)
state: (1, 0), move: ^, action: (1, 0)
state: (1, 2), move: ^, action: (1, 0)
state: (1, 3), move: ., action: (0, 0)
state: (0, 0), move: ^, action: (1, 0)
state: (0, 1), move: >, action: (0, 1)
state: (0, 2), move: ^, action: (1, 0)
state: (0, 3), move: <, action: (0, -1)

Policy Iteration
state: (2, 0), move: >, action: (0, 1)
state: (2, 1), move: >, action: (0, 1)
state: (2, 2), move: >, action: (0, 1)
state: (2, 3), move: ., action: (0, 0)
state: (1, 0), move: ^, action: (1, 0)
state: (1, 2), move: ^, action: (1, 0)
state: (1, 3), move: ., action: (0, 0)
state: (0, 0), move: ^, action: (1, 0)
state: (0, 1), move: >, action: (0, 1)
state: (0, 2), move: ^, action: (1, 0)
state: (0, 3), move: <, action: (0, -1)

Note that both algorithms find the same optimal policy for this problem.

Lec16 - Reinforcement Learning

2009-05-21T17:50:00.000-07:00

Video Lecture - 16

Reinforcement learning (RL) is a particularly interesting field, and it starts from how it is defined: RL is concerned with how an agent ought to take actions in an environment so as to maximize some notion of long-term reward. Artificial intelligence (AI) itself is similarly defined: AI is concerned with the study and design of intelligent agents, where an intelligent agent is a system that perceives its environment and takes actions which maximize its chance of success. So, the problem of RL is as general as the problem of AI; however, because AI has always been a general umbrella for any approach that tried to synthetically reproduce any aspect of what one might consider the whole of human intelligence, it has diverged and fragmented in many more directions than has RL (one can easily get a feel for this by reading the wikipedia entries for RL and AI). Despite its also tumultuous past, the narrower focus of RL has allowed it to remain close to its core formulation of learning by trial and error through interaction with an environment that provides feedback, which has allowed for a fairly stable framework into which ideas from other AI threads could be incorporated (supervised and unsupervised learning systems can be thought of as components of a reinforcement learning system).

Some good resources (in addition to cs229) include AIMA, Sutton and Barto's book, and a couple others that discuss the same principles of RL but in other fields: approximate dynamic programming, neuro-dynamic programming, and dynamic programming and optimal control.

What are some of the things RL has been used for? Well, some pretty amazing things:

Autonomous helicopter flight

Legged robots

Inverted pendulum (a toy problem, but illustrative nevertheless)

Another interesting video, where AIBO learns using RL

Game playing

Automated trading systems

Models of biological neural systems

That last one is interesting because not only do our machine learning methods benefit from an understanding of how the brain works, but those very same methods may help guide our study of the brain. Here's another example of this from some of the work that Tom Mitchell is doing:

Next time, we'll look at some code and illustrate how some of these systems might be implemented.

Hidden Markov Model III

2009-05-18T10:53:00.000-07:00

Last time we used an existing HMM to showcase examples of the inference scenarios in which HMM's might be useful. Today we explore how we might infer the HMM itself from data.

The parameters to be estimated are the transition probabilities and the evidence probabilities, and the inputs we have to work with include the observation sequence and the number of latent states, similar to how the number of gaussians was an input to the mixture of gaussians algorithm. In fact, the standard algorithm for learning a hidden markov model is an instance of the Expectation-Maximization algorithm known as the Baum-Welch algorithm, which uses a slightly modified version of the Forward-Backward algorithm:


#hmm_imp.py
from numpy import *

############################################################
#hmmAlgorithms
############################################################

def fwd(tProb,eProb,prior,ev):
    nstates=shape(prior)[0]
    T=shape(ev)[1]-1
    fMatrix=zeros((nstates,T))

    #our initial forward message equals the prior
    #state distribution (time=0)
    fMatrix[:,0]=prior.reshape(nstates,)

    for t in range(1,T):
        fTemp=zeros((nstates,nstates))
        
        for i in range(nstates):
            for j in range(nstates):
                #for every state i at time t,
                #calculate the probability that each
                #state j at t-1 would have lead to i at t,
                #weighted by the belief that we were in
                #state j at t-1, and weighted by the
                #probability that state i would have
                #produced the evidence we observed at t
                fTemp[i,j]= \
                    tProb[j,i]*fMatrix[j,t-1]*eProb[i,ev[0,t]]

        #collapse fTemp from right to left, summing contributions
        #of each state j at time t-1 to a given state i at t
        fMatrix[:,t]=normEachRow(fTemp.sum(1))

    return fMatrix

#same as fwd, w/o norm
def likeli(tProb,eProb,prior,ev):
    nstates=shape(prior)[0]
    T=shape(ev)[1]-1
    fMatrix=zeros((nstates,T))

    #our initial forward message equals the prior
    #state distribution (time=0)
    fMatrix[:,0]=prior.reshape(nstates,)

    for t in range(1,T):
        fTemp=zeros((nstates,nstates))
        
        for i in range(nstates):
            for j in range(nstates):
                #for every state i at time t,
                #calculate the probability that each
                #state j at t-1 would have lead to i at t,
                #weighted by the belief that we were in
                #state j at t-1, and weighted by the
                #probability that state i would have
                #produced the evidence we observed at t
                fTemp[i,j]= \
                    tProb[j,i]*fMatrix[j,t-1]*eProb[i,ev[0,t]]

        #collapse fTemp from right to left, summing contributions
        #of each state j at time t-1 to a given state i at t
        fMatrix[:,t]=fTemp.sum(1)

    return fMatrix

#same as fwd, but backwards 
def bwd(tProb,eProb,prior,ev):
    nstates=shape(prior)[0]
    T=shape(ev)[1]-1

    #by initializing the entire matrix to ones,
    #the last backward message will be set to all ones,
    #which indicates our diffuse prior belief on our
    #final state.  If we knew what state we ended up on,
    #we could use that information here.
    bMatrix=ones((nstates,T))

    for t in range(T-1,0,-1):
        fTemp=zeros((nstates,nstates))

        #for every state i at time t-1,
        #calculate the probability that it would
        #have lead to state j at time t, weighted
        #by the probabilistic belief that we are
        #in state j at time t, also weighted by
        #the likelihood that j would emit the evidence
        #it emitted at time t
        for i in range(nstates):
            for j in range(nstates):
                fTemp[i,j]=tProb[i,j]*bMatrix[j,t]*eProb[j,ev[0,t]]

        #collapse fTemp from right to left, summing contributions
        #of each state j at time t to a given state i at t-1
        bMatrix[:,t-1]=normEachRow(fTemp.sum(1))

    #allow the backward message to be neutral
    #for the prior state distribution (time=0)
    bMatrix[:,0]=ones((1,nstates))
    return bMatrix

#uses only filtered estimates
def viterbi(tProb,eProb,prior,ev):
    nstates=shape(prior)[0]
    T=shape(ev)[1]-1
    fMatrix=zeros((nstates,T))
    pMatrix=zeros((nstates,T),int)

    #our initial forward message equals the prior
    #state distribution (time=0)
    fMatrix[:,0]=prior.reshape(nstates,)

    for t in range(1,T):
        fTemp=zeros((nstates,nstates))
        
        for i in range(nstates):
            for j in range(nstates):
                #for every state i at time t,
                #calculate the probability that each
                #state j at t-1 would have lead to i at t,
                #weighted by the belief that we were in
                #state j at t-1, and weighted by the
                #probability that state i would have
                #produced the evidence we observed at t
                fTemp[i,j]= \
                    tProb[j,i]*fMatrix[j,t-1]*eProb[i,ev[0,t]]

        #sweep fTemp from right to left, choosing the most likely
        #of all states j at time t-1 to reach state i at t
        fMatrix[:,t]=fTemp.max(1)
        pMatrix[:,t]=fTemp.argmax(1) 
        fMatrix[:,t]=normEachRow(fMatrix[:,t])

    return fMatrix ,pMatrix

def baumwelch(ev,nstates):
    counter=1
    T=shape(ev)[1]-1
    evclass=set(filter(lambda i: i is not None,ev.tolist()[0]))
    nclass=len(evclass)
    tProb=normEachRow(random.rand(nstates,nstates))
    eProb=normEachRow(random.rand(nstates,nclass))
    prior=normEachRow(ones((nstates,)))
    gamma=zeros((nstates,nstates,T))
    eps=zeros((nstates,T))
    
    tProb_old=tProb*0
    eProb_old=eProb*0  
    while (abs(tProb_old-tProb)).max()>1e-3 and \
          (abs(eProb_old-eProb)).max()>1e-3:      
        tProb_old=tProb+0
        eProb_old=eProb+0

        ###########
        #E-step:
        ###########
        
        f=fwd(tProb,eProb,prior,ev)
        b=bwd(tProb,eProb,prior,ev)

        #update gamma
        #for all time t,
        #compute the probability that we transition
        #from state i at time t to state j at time t+1,
        #weighted by the beliefs about where we are at
        #each time and the probability that the evidence
        #seen would have been emitted at time t+1, divided
        #by the smoothed likelihood of the observed sequence
        #observed up until time t would have been generated
        #by the hmm
        for t in range(T-1):
            for i in range(nstates):
                for j in range(nstates):
                    gamma[i,j,t]= \
                        f[i,t]*tProb[i,j]* \
                        b[j,t+1]*eProb[j,ev[0,t+1]]/ \
                        sum(f[i,t]*b[j,t])

        #update eps
        #which is our smoothed state distribution
        #for each time step
        eps=f*b
        eps=transpose(normEachRow(transpose(eps)))

        ###########
        #M-step:
        ###########
                    
        #updated tProb_new
        for i in range(nstates):
            for j in range(nstates):
                #num is the weighted number of times that
                #i has transitioned to j in our evidence sequence
                num=sum([gamma[i,j,t] for t in range(T)])

                #den is the weighted number of times that
                #i has transitioned to any other state in
                #each time t of our observation sequence
                den=0
                for k in range (nstates):
                    den+=sum([gamma[i,k,t] for t in range(T)])
                tProb[i,j]=num/den
                                   
        #update eProb_new
        #weighted probabilitic belief about each state
        #combined with the probability that each state
        #will emit the given evidence
        for j in range(nstates):
            for k in range(nclass):
                num=sum([eps[j,t]
                         for t in range(T)
                         if ev[0,t]==k])
                den=sum([eps[j,t]
                         for t in range(T)])
                eProb[j,k]=num/den
                    
        #update prior
        #sum each state here, which we will normalize later.
        #how about doing some viterbu here and then
        #looking at the most likely path and getting the
        #prior distirbution from there?
        for j in range(nstates):
            prior[j,]=sum([eps[j,t] for t in range(T)])

        tProb=normEachRow(tProb)
        eProb=normEachRow(eProb)
        prior=normEachRow(prior)
        
        print 'iter: %d'%counter
        counter+=1
        
    return prior,eProb,tProb

############################################################
#helperFunctions
############################################################

def normEachRow(pDist):
    dim=shape(pDist)
    rows=dim[0]
    if len(dim)>1 and rows>1:
        for row in range(rows):
            pDist[row,:]=pDist[row,:]/sum(pDist[row,:])
    else:
        pDist=pDist/sum(pDist)
    return pDist

def printMatrix(matrix,msg,start=0,dec=True):
    if not dec:
        format='%d\t'
    else:
        format='%1.3f\t'
    print '\n'+msg
    rows,cols=shape(matrix)
    if cols>10: cols=11
    for col in range(start,cols): print 't=%d\t'%col,
    print '\n',
    for row in range(rows):
        for col in range(start,cols):
            print format%matrix[row,col],
        print '\n',

def accu(a,b):
    assert shape(a)==shape(b)
    total=float(shape(a)[1])
    c=0
    for i in range(shape(a)[1]):
        if a[0,i]==b[0,i]:
            c=c+1
    return c/total

def load(filename):
    inf=open(filename,'r')
    inf.readline()
    data=[]
    for i in inf.readlines():
        x=i.split()
        y=x[0].split(",")
        data.append(y)
    return data

def backtrack(vMatrix,pMatrix):
    rows,cols=shape(vMatrix)
    path=zeros((1,cols),int)
    
    #the most likely state at the last time step
    #is the one with the highest probability in
    #the vMatrix at that time step:
    lastPathCol=shape(path)[1]-1
    lastvMatrixCol=shape(vMatrix)[1]-1
    lastpMatrixCol=shape(pMatrix)[1]-1
    path[0,lastPathCol]=vMatrix[:,lastvMatrixCol].argmax(0)

    #once we have the most likely last state's index, we 
    #use it to look up in the backpointer matrix to
    #trace our way back the most likely path:
    while lastPathCol>0:
        path[0,lastPathCol-1]= \
                pMatrix[path[0,lastPathCol],lastpMatrixCol]
        lastPathCol-=1
        lastpMatrixCol-=1
                                        
    return path

Driven with:


#hmm_inter.py
from hmm_imp import *

############################################################
#weatherModel
############################################################

#state map
stateMap={'sunny':0,'rainy':1,'foggy':2}
stateIndex={0:'sunny',1:'rainy',2:'foggy'}
nStates=len(stateMap)

#observation map
obsMap={'no':0,'yes':1}
obsIndex={0:'no',1:'yes'}

#prior probability on weather states
#P(sunny)=0.50
#P(rainy)=0.25
#P(foggy)=0.25
pi=array([0.50,
          0.25,
          0.25]).reshape(nStates,1)

#transition probabilities of weather states
#                    tomorrrow
#    today     sunny  rainy  foggy
#    sunny      0.8    0.05   0.15
#    rainy      0.2    0.60   0.20 
#    foggy      0.2    0.30   0.50
A=array([[0.8, 0.05, 0.15],
         [0.2, 0.60, 0.20],
         [0.2, 0.30, 0.50]])

#conditional probabilities of evidence given weather
#            P(umbrella=no|weather)  P(umbrella=yes|weather) 
#    sunny          0.9                     0.1 
#    rainy          0.2                     0.8 
#    foggy          0.7                     0.3 
B=array([[0.9, 0.1],
         [0.2, 0.8],
         [0.7, 0.3]])

############################################################
#session
############################################################

#data (evidence has the dummy 'None' entry for time=0)
evidence=array([None,0,0,0,1,0,0,1,1,0,1,0])
classPath=array([2,2,2,1,0,2,1,1,2,1,1])
classPath=classPath.reshape(1,len(classPath))
evidence=evidence.reshape(1,len(evidence))

#compute forward messages
fMatrix=fwd(A,B,pi,evidence)

#normalize probabilities
fMatrix=transpose(normEachRow(transpose(fMatrix)))

#show state distribution as a result of filtering
printMatrix(fMatrix,'fwd matrix')

#compute backward messages
bMatrix=bwd(A,B,pi,evidence)

#show backward messages
printMatrix(bMatrix,'bwd matrix')

#compute smoothed estimate (hindsight)
sMatrix=fMatrix*bMatrix
sMatrix=transpose(normEachRow(transpose(sMatrix)))

#show smooth matrix
printMatrix(sMatrix,'smooth matrix',start=1)

#by summing over the un-normalized state distribution
#resulting from filtering, we obtain the likelihood
#of the hmm having generated the given observation sequence
fMatrix=likeli(A,B,pi,evidence)
lkh=sum(fMatrix[:,shape(fMatrix)[1]-1])
print '\nLikelihood of obs. seq. given original hmm: %f '%lkh

def doViterbi(tProb,eProb,prior, \
              evidence_seq,state_seq,printPath=False):
    print '\nViterbi...'
    
    #compute viterbi matrix and backpointer matrix
    vMatrix,pMatrix=viterbi(tProb,eProb,prior,evidence_seq)

    #show viterbi matrix and backpointer matrix
    printMatrix(vMatrix,'viterbi matrix',start=1)
    printMatrix(pMatrix,'backpointer matrix',start=1,dec=False)

    #backtrack backpointer matrix to find most likely
    #sequence of states
    estPath=backtrack(vMatrix,pMatrix)
    assert shape(state_seq)==shape(estPath)

    if printPath:
        print '\nEstimated class path:\n',estPath
        print '\nActual class path:\n',state_seq
    
    print '\nViterbi decoder accuracy: %1.2f'% \
          accu(state_seq,estPath)

doViterbi(A,B,pi,evidence,classPath,printPath=True)

#load more data from disk
data=load('weather-test1-1000.txt')

#adding 'None' as per note at beginning of session
observations=[None]
states=[]
for c,o in data:
  observations.append(obsMap[o])
  states.append(stateMap[c])

observations=array(observations)
observations=observations.reshape(1,len(observations))
states=array(states)
states=states.reshape(1,len(states))

doViterbi(A,B,pi,observations,states)

print '\nBaum-Welch...'

random.seed(135)
prior,eProb,tProb=baumwelch(observations,3)

print '\nprior:\n',prior
print '\ntProb:\n',tProb
print '\neProb:\n',eProb

doViterbi(tProb,eProb,prior,observations,states)

With output:


fwd matrix
t=0 t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0.500 0.667 0.728 0.752 0.282 0.558 0.679 0.241 0.082 0.447 0.142 
0.250 0.074 0.042 0.034 0.424 0.120 0.055 0.466 0.722 0.198 0.600 
0.250 0.259 0.231 0.214 0.294 0.322 0.265 0.293 0.197 0.356 0.258 

bwd matrix
t=0 t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
1.000 0.448 0.346 0.180 0.404 0.272 0.101 0.152 0.317 0.145 1.000 
1.000 0.228 0.276 0.450 0.241 0.312 0.533 0.473 0.292 0.493 1.000 
1.000 0.324 0.378 0.370 0.355 0.416 0.365 0.376 0.391 0.361 1.000 

smooth matrix
t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0.747 0.718 0.589 0.356 0.470 0.352 0.100 0.083 0.223 0.142 
0.042 0.033 0.067 0.318 0.115 0.151 0.601 0.672 0.335 0.600 
0.210 0.249 0.344 0.326 0.415 0.496 0.300 0.245 0.441 0.258 

Likelihood of obs. seq. given original hmm: 0.000654 

Viterbi...

viterbi matrix
t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0.754 0.858 0.862 0.485 0.737 0.856 0.485 0.198 0.480 0.196 
0.063 0.017 0.012 0.242 0.061 0.019 0.242 0.594 0.240 0.589 
0.183 0.125 0.126 0.273 0.202 0.125 0.273 0.209 0.280 0.215 

backpointer matrix
t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0 0 0 0 0 0 0 0 0 0 
1 2 0 0 1 2 0 1 1 1 
2 0 0 0 2 0 0 2 1 2 

Estimated class path:
[[0 0 0 0 0 0 0 1 1 1 1]]

Actual class path:
[[2 2 2 1 0 2 1 1 2 1 1]]

Viterbi decoder accuracy: 0.36

Viterbi...

viterbi matrix
t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0.754 0.858 0.862 0.485 0.737 0.856 0.485 0.198 0.480 0.196 
0.063 0.017 0.012 0.242 0.061 0.019 0.242 0.594 0.240 0.589 
0.183 0.125 0.126 0.273 0.202 0.125 0.273 0.209 0.280 0.215 

backpointer matrix
t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0 0 0 0 0 0 0 0 0 0 
1 2 0 0 1 2 0 1 1 1 
2 0 0 0 2 0 0 2 1 2 

Viterbi decoder accuracy: 0.61

Baum-Welch...
iter: 1
iter: 2
[...]
iter: 26
iter: 27

prior:
[  7.33520128e-01   2.75132722e-05   2.66452359e-01]

tProb:
[[  6.88340698e-01   9.00867784e-06   3.11650293e-01]
 [  7.81393712e-01   1.87048374e-05   2.18587583e-01]
 [  8.65993461e-01   4.23514682e-05   1.33964187e-01]]

eProb:
[[ 0.68360127  0.31639873]
 [ 0.13036119  0.86963881]
 [ 0.66894045  0.33105955]]

Viterbi...

viterbi matrix
t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0.693 0.693 0.693 0.679 0.693 0.693 0.679 0.679 0.693 0.679 
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 
0.307 0.307 0.307 0.321 0.307 0.307 0.321 0.321 0.307 0.321 

backpointer matrix
t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 
0 0 0 0 0 0 0 0 0 0 
2 2 2 2 2 2 2 2 2 2 
0 0 0 0 0 0 0 0 0 0 

Viterbi decoder accuracy: 0.49

It helps to look at a good diagram of the trellis that results as we unroll the HMM through time; an example of such a diagram can be found in DHS, towards the end of chapter 3. Another good resource to get a good grasp of the HMM algorithms is an interactive spreadsheet that can be found here.

You can check out the estimated transition and emission probability matrices and try to make out how the entries correspond to our original state and evidence vocabulary. In practice, in systems that do things like recognize speech, a lot of human guidance is used when estimating hidden markov models. Running the viterbi algorithm on the observation sequence with the estimated model returns an accuracy of 0.49, which sounds better than it is, since one of the states (the sunny state) is much more prevalent in our dataset (a prior of 0.50), and so if you always guessed that state you would get a score close to the prior of that state given enough data.

Completely unsupervised learning is hard; even humans will come to varying answers in a manual clustering exercise. Reinforcement learning (what we'll be focusing on the next few posts) allows us to use at least some sort of feedback, although not as much as supervised learning, and so should yield better results than learning without any kind of supervision (provided, of course, that a proper reward system is engineered).

Hidden Markov Model II

2009-05-15T08:51:00.000-07:00

Today we look at some of the things we can use an HMM for:

filtering/monitoring
predicting
smoothing
computing likelihood of observation sequence
decoding observation sequence into most likely state sequence

Some good resources for HMM inference algorithms include:

this paper on hidden markov models from the handouts section in cs229.
chapter 15 from AIMA. This is the best introductory resource I know that discusses the inference scenarios listed above, both from a general perspective and from the perspective of specific implementations (hidden markov models, kalman filters, and dynamic bayesian networks).
Berkeley's undergraduate AI class website, which includes an HMM project whose domain model was derived from these notes. This is the domain model (a very simple weather model) that was used in our implementations below.

The Berkeley AI class (cs188) website also includes implementations of these inference algorithms, but in my implementations I sought a simpler presentation to facilitate following chaper 15 of AIMA (too many loops in the solutions posted for cs188 obscure what's going on). I did reuse their helper functions and data, which can be downloaded from the projects section of the website.


#hmm.py

######################################################################
# Weather Model
######################################################################

# state map
weatherStateMap   = {'sunny' : 0, 'rainy' : 1, 'foggy' : 2}
weatherStateIndex = {0 : 'sunny', 1 : 'rainy', 2 : 'foggy'}

# observation map
weatherObsMap   = {'no' : 0, 'yes' : 1}
weatherObsIndex = {0 : 'no', 1 : 'yes'}

# prior probability on weather states
# P(sunny) = 0.5  P(rainy) = 0.25  P(foggy) = 0.25
weatherProb = [0.5, 0.25, 0.25]

# transition probabilities of weather states
#                    tomorrrow
#    today     sunny  rainy  foggy
#    sunny      0.8    0.05   0.15
#    rainy      0.2    0.6    0.2 
#    foggy      0.2    0.3    0.5
weatherTProb = [ [0.8, 0.05, 0.15], [0.2, 0.6, 0.2], [0.2, 0.3, 0.5] ]

# conditional probabilities of evidence (observations) given weather
#                          sunny  rainy  foggy 
# P(umbrella=no|weather)    0.9    0.2    0.7
# P(umbrella=yes|weather)   0.1    0.8    0.3
weatherEProb = [ [0.9, 0.2, 0.7], [0.1, 0.8, 0.3] ]

######################################################################
# Helper Functions
######################################################################

# Using the prior probabilities and state map, return:
#     P(state)
def getStatePriorProb(prob, stateMap, state):
   return prob[stateMap[state]]

# Using the transition probabilities and state map, return:
#     P(next state | current state)
def getNextStateProb(tprob, stateMap, current, next):
   return tprob[stateMap[current]][stateMap[next]]

# Using the observation probabilities, state map, and observation map,
# return:
#     P(observation | state)
def getObservationProb(eprob, stateMap, obsMap, state, obs):
   return eprob[obsMap[obs]][stateMap[state]]

# Normalize a probability distribution
def normalize(pdist):
   s = sum(pdist)
   for i in range(0,len(pdist)):
      pdist[i] = pdist[i] / s
   return pdist

def loadData(filename):
   input = open(filename, 'r')
   input.readline()
   data = []
   for i in input.readlines():
      x = i.split()
      y = x[0].split(",")
      data.append(y)
   return data

def accuracy(a,b):
   total = float(max(len(a),len(b)))
   c = 0
   for i in range(min(len(a),len(b))):
      if a[i] == b[i]:
         c = c + 1          
   return c/total

######################################################################
# HMM inference algorithms 
######################################################################

# Filtering/monitoring.
def filterfwd(stateMap,stateIndex,obsMap,obsIndex,pdist,tprob,eprob,obs):
   pdist_new=[0.0]*len(stateMap)
   for i in range(len(pdist_new)):

      #predict state at 't' based on state belief at 't-1'
      #p(x_t|e_1:t-1) = sum(p(x_t|x_t-1) * p(x_t-1|e_1:t-1)) 
      pdist_new[i]=sum([tprob[j][i]*pdist[j]
                        for j in range(len(stateMap))])
      
      #update prediction with current observation/evidence
      #p(x_t|e_1:t) ~ p(e_t|x_t) * p(x_t|e_1:t-1) 
      pdist_new[i]=eprob[obsMap[obs]][i]*pdist_new[i]
      
   return normalize(pdist_new)

# Prediction.
def predictfwd(stateMap,stateIndex,obsMap,obsIndex,pdist,tprob,eprob):
   pdist_new=[0.0]*len(stateMap)
   for i in range(len(pdist_new)):

      #predict state at 't+1' based on state belief at 't'
      #p(x_t+1|e_1:t) ~ sum(p(x_t+1|x_t) * p(x_t|e_1:t)) 
      pdist_new[i]=sum([tprob[j][i]*pdist[j]
                        for j in range(len(stateMap))])

   return normalize(pdist_new)

# Likelihood.
def likelifwd(stateMap,stateIndex,obsMap,obsIndex,pdist,tprob,eprob,obs):
   pdist_new=[0.0]*len(stateMap)
   for i in range(len(pdist_new)):

      #predict state at 't' based on state belief at 't-1'
      #p(x_t|e_1:t-1) = sum(p(x_t|x_t-1) * p(x_t-1|e_1:t-1)) 
      pdist_new[i]=sum([tprob[j][i]*pdist[j]
                        for j in range(len(stateMap))])
      
      #update prediction with current observation/evidence
      #p(x_t|e_1:t) ~ p(e_t|x_t) * p(x_t|e_1:t-1) 
      pdist_new[i]=eprob[obsMap[obs]][i]*pdist_new[i]
      
   return pdist_new

# Smoothing/hindsight.
def smoothback(stateMap,stateIndex,obsMap,obsIndex,pdist,tprob,eprob,f,obs):
   pdist_new=[0.0]*len(stateMap)
   b=[0.0]*len(stateMap)
   for i in range(len(pdist_new)):

      #smooth state at time 'k' (1<=k<t) based on state belief at 't'    
      #p(e_K+1:t|x+k) = sum(p(e_k+1|x_k+1)*p(e_k+2:t|x_k+1)*p(x_k+1|x_k))
      b[i]=sum([eprob[obsMap[obs]][j] * pdist[j] * tprob[i][j]
                        for j in range(len(stateMap))])
      
      #p(x_k|e_1:t) ~ p(x_k|e_1:k) * p(e_K+1:t|x+k)  
      pdist_new[i]=f[i]*b[i]

   return normalize(pdist_new),b

# Decoding/Viterbi algorithm.
def viterbi(stateMap,stateIndex,obsMap,obsIndex,pdist,tprob,eprob,obs):
   pdist_new=[0.0]*len(stateMap)
   state_idx=[0]*len(stateMap)
   for i in range(len(pdist_new)):

      # p(x_t|x_t-1) * p(x_1,...,x_t-1|e_1:t-1)
      states=[tprob[j][i]*pdist[j] for j in range(len(stateMap))]

      # max(p(x_t|x_t-1) * p(x_1,...,x_t-1|e_1:t-1)) 
      pdist_new[i]=max(states)

      # here, we keep track of:
      # for each possible state in the present (i),
      # what is the most likely state in the past (j)
      # to have lead to this state (i) in the present
      state_idx[i]=states.index(pdist_new[i])
      
      # p(e_t|x_t) * max_val(p(x_t|x_t-1) * p(x_1,...,x_t-1|e_1:t-1))
      pdist_new[i]=eprob[obsMap[obs]][i]*pdist_new[i]
      
   return normalize(pdist_new),state_idx

Driven with:


#hmm_infer.py
from hmm import *

#load sequence of observations from disk
data=loadData('weather-test1-1000.txt')
observations = []
states=[]
for c,o in data:
  observations.append(o)
  states.append(c)

#initialize state distribution to our 
#prior state distribution
stateDist=weatherProb

print '\nPrior distribution over states (t=0):'
for i in range(len(stateDist)):
    print '   ',
    print weatherStateIndex[i],
    print '%1.3f' % stateDist[i],

obsn=10
obs_seq=observations[:obsn]
state_seq=states[:obsn]

print '\n\nUmbrella observation sequence:'
print '   ', obs_seq

####################################################
# Filtering/Monitoring:
#Compute the updated the belief state distribution
#according to the model and the current observation
####################################################
print '\nFiltering/Monitoring...'
for i in range(obsn):
    stateDist=filterfwd(weatherStateMap, \
                        weatherStateIndex, \
                        weatherObsMap, \
                        weatherObsIndex, \
                        stateDist, \
                        weatherTProb, \
                        weatherEProb, \
                        obs_seq[i])
    
    print '\nDistribution over state at t=%d:'%(i+1)
    for i in range(len(stateDist)):
        print '   ',
        print weatherStateIndex[i],
        print '%1.3f' % stateDist[i],

####################################################
# Prediction:
#Forecast state distribution 'n' time steps into
#the future
####################################################
print '\n\nPrediction...'
for i in range(15):
    stateDist=predictfwd(weatherStateMap, \
                         weatherStateIndex, \
                         weatherObsMap, \
                         weatherObsIndex, \
                         stateDist, \
                         weatherTProb, \
                         weatherEProb)
        
    print '\nDistribution over state at t+%d:'%(i+1)
    for j in range(len(stateDist)):
        print '   ',
        print weatherStateIndex[j],
        print '%1.3f' % stateDist[j],

####################################################
# Likelihood:
#Compute the likelihood of our given HMM having
#generated the observed sequence
####################################################
print '\n\nLikelihood...'

#reset our state distribution to our 
#prior state distribution
stateDist=weatherProb

for obs in obs_seq:
    stateDist=likelifwd(weatherStateMap, \
                        weatherStateIndex, \
                        weatherObsMap, \
                        weatherObsIndex, \
                        stateDist, \
                        weatherTProb, \
                        weatherEProb, \
                        obs)
    
print '\nLikelihood of observations given our HMM: ',
print '%1.6f' % sum(stateDist)

####################################################
# Smoothing/Hindsight:
#Compute distribution over past states
#given observations up to present
####################################################
print '\nSmoothing/hindsight...'
#reset our state distribution to our 
#prior state distribution
stateDist=weatherProb

fs=[[]]*obsn
for i in range(obsn):
    stateDist=filterfwd(weatherStateMap, \
                        weatherStateIndex, \
                        weatherObsMap, \
                        weatherObsIndex, \
                        stateDist, \
                        weatherTProb, \
                        weatherEProb, \
                        obs_seq[i])
    
    #record filtered state distribution beliefs
    fs[i]=stateDist
    
    print '\nDistribution over state at t=%d:'%(i+1)
    for j in range(len(stateDist)):
        print '   ',
        print weatherStateIndex[j],
        print '%1.3f' % stateDist[j],

print '\n\nt=%d\n'%obsn
print 'Umbrella observation sequence:'
print '   ', obs_seq

b=[1,1,1]
for i in range(1,obsn):
    stateDist,b=smoothback(weatherStateMap, \
                           weatherStateIndex, \
                           weatherObsMap, \
                           weatherObsIndex, \
                           b, \
                           weatherTProb, \
                           weatherEProb, \
                           fs[-(i+1)], \
                           obs_seq[-i])
        
    print '\nDistribution over state at t-%d:'%(i)
    for j in range(len(stateDist)):
        print '   ',
        print weatherStateIndex[j],
        print '%1.3f' % stateDist[j],

####################################################
# Decoding:
#Compute the most likely sequence of states given
#the observed sequence of observations
####################################################
print '\n\nDecoding...'
#reset our state distribution to our 
#prior state distribution
stateDist=weatherProb

decoded_seq=[0]*obsn
state_index=[[]]*obsn
for i in range(obsn):
    stateDist,state=viterbi(weatherStateMap, \
                            weatherStateIndex, \
                            weatherObsMap, \
                            weatherObsIndex, \
                            stateDist, \
                            weatherTProb, \
                            weatherEProb, \
                            obs_seq[i])
    
    state_index[i]=state

#backtrack to find the most likely state sequence
decoded_seq[obsn-1]=stateDist.index(max(stateDist))
for i in range(obsn-1,0,-1):
  decoded_seq[i-1]=state_index[i][decoded_seq[i]]

#turn state indices into state names
decoded_seq=[weatherStateIndex[idx]
             for idx in decoded_seq]
  
print '\nViterbi - predicted state sequence:\n   ',
print decoded_seq
print 'Viterbi - actual state sequence:\n   ',
print state_seq 
print '\nThe accuracy of the viterbi decoder is', 
print accuracy(state_seq, decoded_seq)

print '\nDecoding full sequence of observations...'
#reset our state distribution to our 
#prior state distribution
stateDist=weatherProb
obsn=len(observations)
decoded_seq=[0]*obsn
state_index=[[]]*obsn

for i in range(obsn):
    stateDist,state=viterbi(weatherStateMap, \
                            weatherStateIndex, \
                            weatherObsMap, \
                            weatherObsIndex, \
                            stateDist, \
                            weatherTProb, \
                            weatherEProb, \
                            observations[i])
    
    state_index[i]=state

#backtrack to find the most likely state sequence
decoded_seq[obsn-1]=stateDist.index(max(stateDist))
for i in range(obsn-1,0,-1):
  decoded_seq[i-1]=state_index[i][decoded_seq[i]]

#turn state indices into state names
decoded_seq=[weatherStateIndex[idx]
             for idx in decoded_seq]
   
print '\nThe accuracy of the viterbi decoder is', 
print accuracy(states, decoded_seq)

With output:


Prior distribution over states (t=0):
    sunny 0.500     rainy 0.250     foggy 0.250 

Umbrella observation sequence:
    ['no', 'no', 'no', 'yes', 'no', 'no', 'yes', 'yes', 'no', 'yes']

Filtering/Monitoring...

Distribution over state at t=1:
    sunny 0.667     rainy 0.074     foggy 0.259 
Distribution over state at t=2:
    sunny 0.728     rainy 0.042     foggy 0.231 
Distribution over state at t=3:
    sunny 0.752     rainy 0.034     foggy 0.214 
Distribution over state at t=4:
    sunny 0.282     rainy 0.424     foggy 0.294 
Distribution over state at t=5:
    sunny 0.558     rainy 0.120     foggy 0.322 
Distribution over state at t=6:
    sunny 0.679     rainy 0.055     foggy 0.265 
Distribution over state at t=7:
    sunny 0.241     rainy 0.466     foggy 0.293 
Distribution over state at t=8:
    sunny 0.082     rainy 0.722     foggy 0.197 
Distribution over state at t=9:
    sunny 0.447     rainy 0.198     foggy 0.356 
Distribution over state at t=10:
    sunny 0.142     rainy 0.600     foggy 0.258 

Prediction...

Distribution over state at t+1:
    sunny 0.285     rainy 0.445     foggy 0.270 
Distribution over state at t+2:
    sunny 0.371     rainy 0.362     foggy 0.267 
Distribution over state at t+3:
    sunny 0.423     rainy 0.316     foggy 0.262 
Distribution over state at t+4:
    sunny 0.454     rainy 0.289     foggy 0.257 
Distribution over state at t+5:
    sunny 0.472     rainy 0.273     foggy 0.255 
Distribution over state at t+6:
    sunny 0.483     rainy 0.264     foggy 0.253 
Distribution over state at t+7:
    sunny 0.490     rainy 0.258     foggy 0.252 
Distribution over state at t+8:
    sunny 0.494     rainy 0.255     foggy 0.251 
Distribution over state at t+9:
    sunny 0.496     rainy 0.253     foggy 0.251 
Distribution over state at t+10:
    sunny 0.498     rainy 0.252     foggy 0.250 
Distribution over state at t+11:
    sunny 0.499     rainy 0.251     foggy 0.250 
Distribution over state at t+12:
    sunny 0.499     rainy 0.251     foggy 0.250 
Distribution over state at t+13:
    sunny 0.500     rainy 0.250     foggy 0.250 
Distribution over state at t+14:
    sunny 0.500     rainy 0.250     foggy 0.250 
Distribution over state at t+15:
    sunny 0.500     rainy 0.250     foggy 0.250 

Likelihood...

Likelihood of observations given our HMM:  0.000654

Smoothing/hindsight...

Distribution over state at t=1:
    sunny 0.667     rainy 0.074     foggy 0.259 
Distribution over state at t=2:
    sunny 0.728     rainy 0.042     foggy 0.231 
Distribution over state at t=3:
    sunny 0.752     rainy 0.034     foggy 0.214 
Distribution over state at t=4:
    sunny 0.282     rainy 0.424     foggy 0.294 
Distribution over state at t=5:
    sunny 0.558     rainy 0.120     foggy 0.322 
Distribution over state at t=6:
    sunny 0.679     rainy 0.055     foggy 0.265 
Distribution over state at t=7:
    sunny 0.241     rainy 0.466     foggy 0.293 
Distribution over state at t=8:
    sunny 0.082     rainy 0.722     foggy 0.197 
Distribution over state at t=9:
    sunny 0.447     rainy 0.198     foggy 0.356 
Distribution over state at t=10:
    sunny 0.142     rainy 0.600     foggy 0.258 

t=10

Umbrella observation sequence:
    ['no', 'no', 'no', 'yes', 'no', 'no', 'yes', 'yes', 'no', 'yes']

Distribution over state at t-1:
    sunny 0.223     rainy 0.335     foggy 0.441 
Distribution over state at t-2:
    sunny 0.083     rainy 0.672     foggy 0.245 
Distribution over state at t-3:
    sunny 0.100     rainy 0.601     foggy 0.300 
Distribution over state at t-4:
    sunny 0.352     rainy 0.151     foggy 0.496 
Distribution over state at t-5:
    sunny 0.470     rainy 0.115     foggy 0.415 
Distribution over state at t-6:
    sunny 0.356     rainy 0.318     foggy 0.326 
Distribution over state at t-7:
    sunny 0.589     rainy 0.067     foggy 0.344 
Distribution over state at t-8:
    sunny 0.718     rainy 0.033     foggy 0.249 
Distribution over state at t-9:
    sunny 0.747     rainy 0.042     foggy 0.210 

Decoding...

Viterbi - predicted state sequence:
    ['sunny', 'sunny', 'sunny', 'sunny', 'sunny', 'sunny', 'rainy', 'rainy', 'rainy', 'rainy']
Viterbi - actual state sequence:
    ['foggy', 'foggy', 'foggy', 'rainy', 'sunny', 'foggy', 'rainy', 'rainy', 'foggy', 'rainy']

The accuracy of the viterbi decoder is 0.4

Decoding full sequence of observations...

The accuracy of the viterbi decoder is 0.636

An interesting thing to note is how predicting into the future will eventually lead back to the stationary distribution of the markov chain. Next time, we will talk about how to actually infer the parameters of our HMM from data.

Hidden Markov Model

2009-05-13T12:02:00.000-07:00

In our last post we talked about one particular application of markov chains: sampling from probability distributions. Another example was that of a random text generator (RTG), and as trivial a pursuit as it may be, an RTG is a good simple application from which to explore these concepts.

Imagine that we have a different kind of RTG, one that instead of learning word transition probabilities, learns part-of-speech (POS) transition probabilities. One categorization of POS may include: {'noun', 'verb', 'pronoun', 'preposition', 'adverb', 'conjunction', 'participle', 'article'}, but there are other variations. Imagine also that we have a POS distribution ('noun'=.3, 'verb'=.3, ... all adding up to one, of course), and a model that approximates the distribution of words given a certain POS, so that given that POS='verb', we would have that: ('eat'=.1, 'sleep'=.3, 'work'=.3, play='.3'), which are the verbs in our vocabulary. If we had that info, we could generate random text as follows:

Start with a random POS drawn from the POS distribution
Draw a random word from the given POS' word distribution
Transition to the next POS given the POS transition probability matrix
Draw a new word from the new POS' word distribution
Repeat step 3

If we generate some random text using this model and show it to someone, what they would observe (the sequence of words) would be called an observed markov chain. Unless we show them the sequence of POS tags, they would not know about it, and so that would be an unobserved (hidden) markov chain to them. These kind of models turn out to be useful in all kinds of applications, speech recognition for example, where a computer sees an observed markov chain of sound recordings and it tries to infer the what the unobserved markov chain of words that "generated those sounds" actually was. As a side note, Stanford Engineering Everywhere also includes a course on Natural Language Processing. Here's a link to Berkeley's version of that course.

A hidden markov model (HMM) is useful for a lot of things, but in essence, we learn the parameters of a HMM based on the available data (any priors and the sequence of observations) so that we can use it to calculate the most probable sequence of [latent] states given the sequence of observations (chapter 15 of AIMA has an interesting discussion on things we can do with temporal models such as HMM).

Are probabilites of sequences really useful? It turns out they are, think about a cell phone conversation you've recently had where the reception wasn't great, and the speaker on the other end said something you couldn't really understand. Imagine, for example, that you heard '?all', but you weren't sure what letter should've gone into the place of the question mark; did they say 'ball', 'call', or 'wall'? A silly way to decide what the word was would be to just use the probability of ever hearing each of those words and picking the most probable one, but then you're not taking advantage of the context of the conversation, which includes all the other words that you heard before '?all' in the conversation. Most likely, what you would do would be to "calculate" the probabilities of each word given that you were talking about a soccer game, and chose 'ball' as most likely, in other words, you would calculate probabilities of sequences. We will look into estimating the parameters of the HMM and then into calculating probabilities of sequences, but first, we examine the simpler case of doing so in observed markov chains.

One field very concerned with sequences is bioinformatics, which uses computational methods to solve problems in biology, and where common activities include mapping and analyzing DNA and protein sequences and aligning different DNA and protein sequences to compare them. For many applications, a DNA sequence can be treated as a string composed of the characters in the set {A,G,C,T}, and one simple example of something that can be done is to try and come up with a model of how DNA sequences evolve through time. For that, several substitution models have been proposed, including a number of different markov models. In order to infer the parameters of a markov chain using data, we would use a technique we have seen before, maximum likelihood, which is well described here, and which unsurprisingly yields that the estimates are the ratios present in the dataset:


#mc.py
from numpy import *

############################################################
#compute probabilities of sequences given a markov chain
############################################################
def printchain(chain):
    for node in chain:
        print node+':',
        print chain[node]
    
def nextprob(chain,node,nextnode):
    if node in chain and nextnode in chain[node].keys():
        return chain[node][nextnode]
    else:
        return 0.0

def seqprob(chain,seq):
    probseq=nextprob(chain,'b',seq[0])
    for index in range(len(seq)-1):
        probseq*=nextprob(chain,seq[index],seq[index+1])
    return probseq

chain={}
chain['b']={'A':1.0}
chain['A']={'C':0.3,'G':0.7}
chain['C']={'T':1.0}
chain['G']={'T':1.0}
chain['T']={'e':1.0}
chain['e']=None

print '\nprob from A to C:\t%f'%nextprob(chain,'A','C')
print 'prob of seq AGT:\t%f'%seqprob(chain,'AGT')
printchain(chain)

############################################################
#infer markov chain from given sequences 
############################################################

def buildchain(seqs):
    #initialize chain
    chain={'b':{},'e':None}

    #add sequence items as nodes
    for seq in seqs:
        for item in seq:
            chain.setdefault(item,{})
    for i in range(len(seqs[0])-1):
        for seq in seqs:
            chain[seq[i]].setdefault(seq[i+1],0.0)
            chain[seq[i]][seq[i+1]]+=1

    #add limit nodes
    for seq in seqs:
        chain['b'].setdefault(seq[0],0.0)
        chain['b'][seq[0]]+=1
        chain[seq[len(seq)-1]].setdefault('e',1.0)
        
    #normalize each node's transition probabilities
    for node in chain:
        if node<>'e':
            norm=1/sum(chain[node].values())
            for nextnode in chain[node]:
                chain[node][nextnode]*=norm
        
    return chain
            
seqs=['ACT','AGT']
chain=buildchain(seqs)

print '\nprob of seq ACT:\t%f'%seqprob(chain,seqs[0])
print 'prob of seq AGT:\t%f'%seqprob(chain,seqs[1])
printchain(chain)

For which we get:


prob from A to C: 0.300000
prob of seq AGT: 0.700000
A: {'C': 0.29999999999999999, 'G': 0.69999999999999996}
C: {'T': 1.0}
b: {'A': 1.0}
e: None
G: {'T': 1.0}
T: {'e': 1.0}

prob of seq ACT: 0.500000
prob of seq AGT: 0.500000
A: {'C': 0.5, 'G': 0.5}
C: {'T': 1.0}
b: {'A': 1.0}
e: None
G: {'T': 1.0}
T: {'e': 1.0}

Notice anything wrong with this picture? Watch what happens when we use the dataset below to infer a new chain and calculate probabilities of sequences:


seqs=['GTATC','GTTTC','GT-AC','GT-AC']
chain=buildchain(seqs)

print '\nprob of seq GTATC:\t%f'%seqprob(chain,seqs[0])
print 'prob of seq GTTTC:\t%f'%seqprob(chain,seqs[1])
print 'prob of seq GT-AC:\t%f'%seqprob(chain,seqs[2])
printchain(chain)

We get:


prob of seq GTATC: 0.013605
prob of seq GTTTC: 0.023324
prob of seq GT-AC: 0.190476
A: {'C': 0.666666666667, 'T': 0.333333333333}
C: {'e': 1.0}
b: {'G': 1.0}
e: None
G: {'T': 1.0}
-: {'A': 1.0}
T: {'A': 0.142857142857, 'C': 0.285714285714, '-': 0.285714285714, 'T': 0.285714285714}

When we should have gotten 25% probability for sequences GTATC and GTTTC respectively and 50% probability for sequence GT-AC.

What our markov chain model failed to capture was that given the dataset ['GTATC','GTTTC','GT-AC','GT-AC'], and assuming fixed length sequences of size 5, the probability of 'T' being followed by an item in our alphabet ['A','G','C','T','-'] changes depending on whether the 'T' we are looking at falls into the 2nd, 3rd, or 4th position in our sequence. Doing some manual calculations, we can see that we have the following transition probabilities:

when 'T' is in position 2: {'A':.25,'T':.25,'-':.5}
when 'T' is in position 3: {'T':1.0}
when 'T' is in position 4: {'C':1.0}

A hidden markov model would actually model these transitions effectively, where the "hidden" chain would be the sequence of positions. Actually, this is a very simple hidden markov model, because of the transition probabilities for the hidden chain (i.e. position 1 transitions with 100% probability into position 2, and so on), and so a few modifications to our 'buildchain' and 'seqprob' methods up there would fix our problem. In our next post, however, we will look at a more general algorithm to learn and compute with HMM's.

Markov Chain Monte Carlo II

2009-05-12T15:12:00.000-07:00

Continuing our discussion on simulation, today we discuss a class of algorithms for sampling from probability distributions based on constructing a markov chain that has the desired distribution as its equilibrium distribution: markov chain monte carlo (MCMC).

What is a markov chain's equilibrium distribution? For illustration purposes, take the following example, from Applied Linear Algebra and Matrix Analysis:

"Suppose two toothpaste companies compete for customers in a fixed market in which each customer uses either brand A or brand B. Suppose also that a market analysis shows that the buying habits of the customers fit the following pattern: each quarter, 30% of A users will switch to B, while the rest stay with A. During the same time, 40% of B users will switch to A, while the rest will stay with B."

As the python session below shows, it does not matter what the initial state distribution is (i.e. what percentage of customers uses brand A and B), after enough quarters, the state distribution will converge to an equilibrium:


#tpaste_inter.py
from numpy import *

def evolve(x,T,n):
    for i in range(n):
        x=dot(T,x)
    print str(x)+'\n' 

#transition matrix
T=array([[.7,.4],
         [.3,.6]])

#one initial state distribution
x0=array([.5,
          .5]).reshape(2,1)

evolve(x0,T,20)

#another initial state distribution
x0=array([1,
          0]).reshape(2,1)

evolve(x0,T,20)

#yet another initial state distribution
x0=array([0,
          1]).reshape(2,1)

evolve(x0,T,20)

Output:


[[ 0.57142857]
 [ 0.42857143]]

[[ 0.57142857]
 [ 0.42857143]]

[[ 0.57142857]
 [ 0.42857143]]

As shown, all the initial distributions converged to A retaining 57% of customers, with B retaining the remaining 43%. This convergence to an equilibrium distribution does not always happen in markov chains (check this out for necessary conditions for convergence), but we do need it in MCMC. PageRank is also an example of this, where the transition matrix is a link graph and the equilibrium state distribution denotes the rank of each page.

The idea in MCMC, is that we need to construct a transition matrix so that its equilibrium distribution represents the target distribution we want to sample from, so that independently of what we initialize the initial state distribution to, by applying the transition matrix enough times, we will converge to our desired distribution. One obvious benefit of this method over say rejection sampling is that in rejection sampling, we are sampling all over the distribution space (the proposal distribution's space, that is), which may generate a lot of rejections (wasted samples). The same is true with importance sampling (sampling importance resampling, that is). With MCMC, on the other hand, we are smarter about navigating the state space of the desired distribution by navigating the links on its chain. For example, imagine we are sapling from some bivariate distribution, one approach is to use use something like rejection sampling and to generate a bunch of i.i.d (x0,x1) tuples; while another approach would be to start with a random sample for 'x0', then generate 'x1' from a distribution of p(x1|x0), which gives us our first tuple (x0,x1), then we generate a new sample for 'x0' based on our last 'x1' by drawing this time from the conditional p(x0|x1), which gives us our next tuple (not i.i.d), and then we generate a new 'x1' based on this last 'x0' for the next tuple, and so on. This last approach is known as Gibbs Sampling, and is one variant of MCMC.

A more general form of MCMC is the so called Metropolis-Hastings algorithm, where rather than always taking the next link in the chain as a valid sample, we only take it if it is a much better sample than the last sample (according to some fitness function). If it is not much better, we don't take it. If they are about the same (similar fitness scores), we flip a coin to decide. Good introductory resources include Nando's paper on MCMC and chapter 11 of Bishop's book. In there you will find more details about how we go about constructing a transition matrix (the proposal conditional distribution) with the desired properties, which can be considered somewhat of an art (different choices for the proposal distribution lead to different variants of the Metropolis-Hastings algorithms, one of which is Gibbs Sampling).

Below is an example from Marsland's book, where he uses Metropolis-Hastings to sample from a mixture of two gaussians based on a proposal distribution that is a gaussian of the form N(x,b^2) for some b>0. This means that the proposal is a drawn from a Normal centered at the current value. This variant of Metropolis-Hastings is known as the "Random-walk-Metropolis-Hastings" algorithm, the reason for the name being that if we did not do the accept-reject step (evaluating the new sample by a fitness function) we would then be simulating a random walk:


#mcmc.py

# Code from Chapter 14 of Machine Learning: An Algorithmic Perspective
# by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008

# The Metropolis-Hastings algorithm
from pylab import *
from numpy import *

def p(x):
    mu1 = 3
    mu2 = 10
    v1 = 10
    v2 = 3
    return 0.3*exp(-(x-mu1)**2/v1) + 0.7* exp(-(x-mu2)**2/v2)

stepsize = 0.5
x = arange(-10,20,stepsize)
px = zeros(shape(x))
for i in range(len(x)):
    px[i] = p(x[i])
N = 5000

# random walk chain
u2 = random.rand(N)
sigma = 10
y2 = zeros(N)
#start with a random sample
y2[0] = random.normal(0,sigma)
for i in range(N-1):
    #get next sample in the chain
    y2new = y2[i] + random.normal(0,sigma)
    #apply fitness function
    alpha = min(1,p(y2new)/p(y2[i]))
    #only take this sample if better than last
    if u2[i] < alpha:
        y2[i+1] = y2new
    else:
        y2[i+1] = y2[i]

figure(1)
nbins = 30
hist(y2, bins = x)
plot(x, px*N/sum(px), color='r', linewidth=2)

show()

And for which we get:

Markov Chain Monte Carlo

2009-05-07T18:31:00.000-07:00

We move quicker now towards the reinforcement learning part of the lectures, where the concept of a markov decision process (MDP) is introduced and further discussed (starting in lecture 16). Before we get there though, we should look into a few related concepts that will allow us to build up to MDP's, particularly, simpler markov models: markov chains and hidden markov models.

To see how all these different flavors of discrete time, discrete state markov models are related, the following chart can be useful (found here):

Markov Models		Do we have control over the state transitions?
Markov Models		No	Yes
Are the states completely observable?	Yes	Markov Chain	Markov Decision Process
Are the states completely observable?	No	Hidden Markov Model	Partially Observable Markov Decision Process

So, what is a markov chain? If you think about our diagram relating probability vs. inferential statistics, there are 2 ways from which we can look at markov chains. The first one is to think about the data generating process, and if we did so, we would say that a markov chain is a system whose initial state (x0) is a probability distribution vector, and whose transition matrix (A) is composed of columns which are themselves probability distribution vectors (a stochastic matrix). This system evolves through time, which we consider to be discrete, by taking the initial state vector and "passing it through" the transition matrix, generating the next state vector (x1), and so on, for state vectors (x2), (x3), ..., (xn). After 'k' iterations, we would have seen 'k' state vectors, and we can think of this sequence of state vectors as the chain, with the state vectors being the links. This leads us to the second way of looking at markov chains, which answers the question: what can we say about the chain (the sequence of observations) that is generated? Well, the links in the chain have the so called markov property, which says that a future state vector (xk+1) is independent of all past state vectors (x0 through xk-1) given the present state vector (xk).

The markov property we mentioned characterizes a so called 1st-order markov model. You could have a 2nd-order model, where a state vector depends only on the 2 previous states, etc, until you have a model where your state vector depends on all previous states, and nth-order model. There are many sequence of events where successive observations are correlated, and so this markov property is simply a [formalism|simplification] of that fact. Obviously, the smaller the order of the markov model, the more tractable it is, and fortunately, just as the naive bayes assumption works surprisingly well in a variety of applications, so do the 1st-order models. Check the wikipedia entry on markov chains for a list of applications and for further references, including this book, whose earlier edition is available for download. Gentler introductions can be found in linear algebra books or in probability books that have a chapter or two on stochastic processes.

There are many simple examples of markov chains, and the most prevalent one throughout the web appears to be that of a random text generator; check this python script out for an example of an algorithm described in chapter 3 of The Practice of Programming. What we will do instead, is explore an application of markov chains to simulations, the Markov chain Monte Carlo algorithm.

We've very briefly talked about markov chains, but what about monte carlo methods? As wikipedia mentions, monte carlo method is a general term for basically any algorithm that uses random samplings to calculate its results, and one desired property of these algorithms is that as the number of samplings increase, so does the accuracy of the calculated results. There is another kind of randomized algorithms, the las vegas algorithms, that despite their use of random numbers, they always give the correct result (think quicksort). One illustrative example of a monte carlo method is to calculate pi. There are some good explanations here and here, the latter one also has a nice little applet of the simulation at work.

Simulation from probability distributions relies on the ability to generate random numbers from a (0,1) uniform distribution, which can be accomplished with the linear congruential method:


#lcm_imp.py
import time

m=2**32
a=1664525
c=1013904223
xi=time.time()

def seed(x):
    global xi
    xi=x

def rng():
    global xi
    xi=(a*xi+c)%m
    return xi/float(m)

Driven with:


#lcm_inter.py
from lcm_imp import *

seed(123)
for i in range(10): print rng()

Yielding:


0.283736921381
0.43513002363
0.0386512577534
0.220879904693
0.359427076299
0.590244138846
0.361280900426
0.326849908335
0.07973951241
0.647962252842

Having that, we can simulate other distributions by transforming the samples from the (0,1) uniform distribution using the inverse of the distribution function of the density we want to sample from (inverse transform method). For example, say we want to sample from an (a,b) uniform distribution where a=-1 and b=1. We know that the density of a continuous uniform distribution is f(x)=1/(a-b) in the interval a<=x<=b, and its distribution function is F(x)=(x-a)/(b-a) for some a<=x<b. The inverse of this distribution function is x=(b-a)F(x)+a and so if we plug in the the resulting samples for U(0,1) into F(x) in the inverse distribution, we will obtain results which are distributed U(a,b) for whatever 'a' and 'b' we choose:


#invtrans_inter.py
from lcm_imp import *

seed(123)

base_U_samples=[rng() for i in range(10)]

a,b=-1,1

desired_U_samples=[(b-a)*x+a for x in base_U_samples]

for item in desired_U_samples: print item

For which we get:


-0.432526157238
-0.129739952739
-0.922697484493
-0.558240190614
-0.281145847403
0.180488277692
-0.277438199148
-0.34630018333
-0.84052097518
0.295924505685

The efficacy of the inverse transform method relies on our ability to calculate the desired inverse distribution. For distributions where we can't do so, we could turn to another method, which although not very efficient in high-dimensions, it's easy to understand: rejection sampling. This technique is nicely explained in part (2/6) of Nando's lectures on Monte Carlo simulation (with demo and all). For illustrative purposes, here is an example from Stephen Marsland's recent Machine Learning: An Algorithmic Perspective book:


#rejectsamp_imp.py

# Code from Chapter 14 of Machine Learning: An Algorithmic Perspective
# by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008

# The basic rejection sampling algorithm

from pylab import *
from numpy import *

def qsample():
    return random.rand()*4.

def p(x):
    return 0.3*exp(-(x-0.3)**2) + 0.7* exp(-(x-2.)**2/0.3) 

def rejection(nsamples):
    
    M = 0.72#0.8
    samples = zeros(nsamples,dtype=float)
    count = 0
    for i in range(nsamples):
        accept = False
        while not accept:
            x = qsample()
            u = random.rand()*M
            if p(x)>u:
                accept = True
                samples[i] = x
            else: 
                count += 1
    print count   
    return samples

Driven with:


#rejectsamp_inter.py
from rejectsamp_imp import *

x = arange(0,4,0.01)
x2 = arange(-0.5,4.5,0.1)
realdata = 0.3*exp(-(x-0.3)**2) + 0.7* exp(-(x-2.)**2/0.3)
box = ones(len(x2))*0.75#0.8
box[:5] = 0
box[-5:] = 0
plot(x,realdata,'k',lw=6)
plot(x2,box,'k--',lw=6)

import time
t0=time.time()
samples = rejection(10000)
t1=time.time()
print "Time ",t1-t0

hist(samples,15,normed=1,fc='k')
xlabel('x',fontsize=24)
ylabel('p(x)',fontsize=24)
axis([-0.5,4.5,0,1])
show()

And which yields:


17857
Time  0.555863857269

Notice how it took 17875 tries to get 10000 samples from our desired distribution.

Numpy gotcha!

2009-05-07T09:39:00.000-07:00

So, in my haste to push my last post on unsupervised learning and move on to sequential data, I did not check the results as thoroughly as I should have, and so when my brother looked at them, he wisely pointed out that the separations among the clusters found by the expectation maximization application to mixture of gaussians (MoG) looked linear, which is what you'd expect from k-means but not from MoG.

I looked at it and sure enough, there were a few minor issues, which once fixed, uncovered another issue that took me a little longer to figure out. What was happening now was that all the points were being labeled as belonging to one cluster/gaussian by the MoG algorithm. I eventually narrowed it down to the cause, which was that my implementation was estimating covariance matrices with very large values, and so when it came time to find the most probable cluster for all points, the one with the largest covariance matrix was overpowering the remaining gaussians.

There was nothing wrong with the implementation per se, instead, what slipped by me was this particular instance of how numpy arrays are handled:


>>> from numpy import *
>>> a=random.ranf((5,))
>>> b=random.ranf((5,1))
>>> a*b
array([[ 0.04429077,  0.03425297,  0.63137759,  0.00791116,  0.47450932],
       [ 0.04375046,  0.03383511,  0.62367522,  0.00781465,  0.46872064],
       [ 0.02370387,  0.01833176,  0.3379054 ,  0.00423395,  0.25395146],
       [ 0.03411096,  0.02638025,  0.48626142,  0.00609285,  0.36544784],
       [ 0.02545382,  0.01968512,  0.36285147,  0.00454653,  0.27269958]])
>>> c=random.ranf((5,1))
>>> d=random.ranf((5,1))
>>> c*d
array([[ 0.18859163],
       [ 0.33695256],
       [ 0.34911791],
       [ 0.09585323],
       [ 0.20524801]])
>>> e=a*b
>>> shape(e)
(5, 5)
>>> shape(e[:,2])
(5,)
>>> e
array([[ 0.04429077,  0.03425297,  0.63137759,  0.00791116,  0.47450932],
       [ 0.04375046,  0.03383511,  0.62367522,  0.00781465,  0.46872064],
       [ 0.02370387,  0.01833176,  0.3379054 ,  0.00423395,  0.25395146],
       [ 0.03411096,  0.02638025,  0.48626142,  0.00609285,  0.36544784],
       [ 0.02545382,  0.01968512,  0.36285147,  0.00454653,  0.27269958]])
>>> e[:,2]
array([ 0.63137759,  0.62367522,  0.3379054 ,  0.48626142,  0.36285147])

A few things to note there:

The product of 2 arrays with dimensions (5,1) results in an element-wise multiplication. OK, nothing wrong here.
An array of dimensions (5,) multiplied with another one with dimensions (5,1) results in a (5,5) array. Not sure why a (5,) wouldn't default to (5,1).
Slicing the 3rd column of the array 'e' up there results in a (5,) array as opposed to a (5,1), which is what I would have expected.

Anyways, a simple 'reshape' can bring the array back to what you want, so you just gotta be aware of these things.

Also see: array or matrix, which should I use?

Lec12 - Expectation Maximization

2009-05-05T06:06:00.000-07:00

The Expectation Maximization (EM) algorithm is covered in lectures 12, 13, and 14, from its general form to a few derived instances (mixture of gaussians, mixture of naive bayes, and factor analysis).

We will illustrate unsupervised learning with implementations of mixture of gaussians (an instance of EM) and k-means. K-means approaches the problem from a geometric perspective (at each iteration, a data point is assigned to one and only one cluster, the closest one), while mixture of gaussians does so from a probabilistic perspective (a data point belongs to all clusters, with a given probability):


#em_imp.py
from numpy import *
import sys

def gen_data(sources,k,m):
    save3Dplot('sources',sources)
    kdata=[zeros((m/k,2)) for i in range(k)]
    for i in range(k):
        for j in range(2):
            mu=sources[i][j]['mu']
            sd=sources[i][j]['sd']
            kdata[i][:,j]=random.normal(mu,sd,m/k)
        saveplot('X'+str(i),kdata[i],' ')

    data=kdata[0]
    for i in range(1,k):
        data=concatenate((data,kdata[i]))

    random.shuffle(data)
    return data

def saveplot(fname,data,sep):
    m,n=data.shape
    fp=file(fname,'w')
    for i in range(m):
        for j in range(n):
            fp.write(str(data[i,j])+sep)
        fp.write('\n')
    fp.close()

def plot(name,clusters,centroids,X):
    saveplot(name+'centroids',centroids,' ')
    k=shape(centroids)[0]
    for i in range(k):
        data=X.take(where(clusters==i)[0],0)
        saveplot(name+'data'+str(i),data,' ')

def dnorm(X,muvec,covmat):
    m=shape(X)[0]
    probs=zeros((m,1))
    dt=linalg.det(covmat)
    A=1./((2*pi)*sqrt(dt))
    covinv=linalg.inv(covmat)
    for i in range(m):
        t=X[i]-muvec
        temp=-.5*dot(dot(t,covinv),transpose(t))
        probs[i,0]=A*exp(temp)
    return probs

def save3Dplot(fname,sources):
    for src in range(len(sources)):
        mu1=sources[src][0]['mu']
        sd1=sources[src][0]['sd']
        mu2=sources[src][1]['mu']
        sd2=sources[src][1]['sd']
        covmat=zeros((2,2))
        covmat[0,0]=sd1**2
        covmat[1,1]=sd2**2
        avg=array([mu1,mu2])     
        data=zeros((100,100))
        
        for i in range(100):
            for j in range(100):
                xvec=array([i,j])/10.
                data[i,j]=dnorm(xvec.reshape(1,2),avg,covmat)
        
        saveplot(fname+str(src),data,'\n')

def kmeans(X,k):
    m=shape(X)[0]
    n=shape(X)[1]
    oldcentroids=zeros((k,n))
    e=1

    #pick actual points at random to be initial centroids
    centroids=X[(random.ranf((k,1))*m).astype('int')]
    centroids=centroids.reshape(k,n)
    clusters=zeros((m,1),'int')

    #repeat until convergence
    while linalg.norm(oldcentroids-centroids)>1e-15:
        oldcentroids=centroids+0 #copy

        #compute cluster memberships
        for i in range(m):
            #repeat X[i,:] k times
            dpoints=kron(ones((k,1)),X[i,:])
            dists=sum((dpoints-centroids)**2,1).reshape(k,1)
            clusters[i,0]=argmin(dists)

        #compute new cluster centroids
        for i in range(k):
            centroids[i,:]=mean(X.take(where(clusters==i)[0],0),0)

        print 'iter: %d, norm: %s'%(e,linalg.norm(oldcentroids-centroids))
        e+=1

    #compute final cluster memberships
    for i in range(m):
        dpoints=kron(ones((k,1)),X[i,:])
        dists=sum((dpoints-centroids)**2,1).reshape(k,1)
        clusters[i,0]=argmin(dists)
    
    return clusters,centroids

def MoG(X,k,kclusters,kcentroids):
    m=shape(X)[0]
    n=shape(X)[1]
    W=zeros((m,k))
    lastll=sys.maxint
    ll=0  
    e=1
    
    #first guess parameters (using results from kmeans)
    priors=zeros((k,1))
    muvecs=zeros((k,n))
    covmats=[zeros((n,n)) for i in range(k)]
    for i in range(k):
        priors[i,0]=len(X.take(where(kclusters==i)[0],0))/float(m)
        covmats[i]=cov(transpose(X.take(where(kclusters==i)[0],0)))
        muvecs[i]=kcentroids[i]
        
    #repeat until convergence
    while abs(lastll-ll)>1e-5 and lastll>ll:
        
        #e-step: compute w's (probability of membership in each k)
        for i in range(k):
            #posterior~likelihood*prior
            W[:,i:i+1]=dnorm(X,muvecs[i,:],covmats[i])*priors[i,0]
        #normalize
        W=W/W.sum(1).reshape(m,1)

        #m-step: update parameters based on last membership
        for i in range(k):
            priors[i,0]=sum(W[:,i])/float(m)
            muvecs[i,:]=sum(W[:,i:i+1]*X,0).reshape(1,n)/sum(W[:,i])
            for j in range(n):
                for l in range(n):
                    weight=W[:,i:i+1]
                    diff1=(X[:,j]-muvecs[i,j]).reshape(m,1)
                    diff2=(X[:,l]-muvecs[i,l]).reshape(m,1)
                    covmats[i][j,l]=sum(weight*diff1*diff2)/sum(W[:,i])
            
        #compute log-likelihood
        ll1=zeros((m,1))
        ll2=W.argmax(1).reshape(m,1)
        for i in range(m):
            ll1[i,:]=dnorm(X[i].reshape(1,n),muvecs[ll2[i],:],covmats[ll2[i]])
            ll2[i,0]=priors[ll2[i]]
        ll1=(1e-15>ll1).choose(ll1,1e-15)
        ll2=(1e-15>ll2).choose(ll2,1e-15)
        lastll=ll
        ll=sum(log(ll1)+log(ll2))   
        print 'iter: %d, log-likelihood: %s'%(e,ll)
        e+=1

    #compute final membership
    clusters=W.argmax(1).reshape(m,1)
    
    return clusters,muvecs

Driven with:


#em_inter.py
from em_imp import *
import os

#build your [bivariate] gaussians
sources=[[{'mu': 7.0, 'sd': 0.7}, {'mu': 7.0, 'sd': 0.7}],
         [{'mu': 3.0, 'sd': 0.7}, {'mu': 7.0, 'sd': 0.7}],
         [{'mu': 5.0, 'sd': 1.0}, {'mu': 4.0, 'sd': 1.0}]]

k=len(sources) #number of gaussians
m=1500*k #number of data points
random.seed(123456) #for repeatability

#sample data from sources
X=gen_data(sources,k,m)

#save files for data plot
saveplot('X',X,' ')

print '\nkmeans:'

#cluster data using k-means
kclusters,kcentroids=kmeans(X,3)

#save files for kmeans plot
plot('k',kclusters,kcentroids,X)

print '\nem:'

#cluster data using MoG
MoGclusters,MoGcentroids=MoG(X,3,kclusters,kcentroids)

#save files for MoG plot
plot('MoG',MoGclusters,MoGcentroids,X)

print '\nkcentroids:'
print kcentroids

print '\nMoGCentroids:'
print MoGcentroids

#run gnuplot script
os.system('gnuplot plotinst')

Using the following gnuplot script:


#plotinst
set term gif
set output 'source_contours.gif' 
unset surface
set contour
set view 0,0,1,1                                  
splot 'sources0','sources1','sources2'            
set output 'data_color.gif'           
plot 'X0' using 2:1, 'X1' using 2:1, 'X2' using 2:1
set output 'data_nocolor.gif'                     
plot 'X' using 2:1                                            
set output 'sources3d.gif'
set hidden3d
set surface                           
set view 60,25,1,1                     
splot 'sources0','sources1','sources2'
set output 'kdata_color.gif'           
plot 'kdata0' using 2:1, 'kdata1' using 2:1, 'kdata2' using 2:1, 'kcentroids' using 2:1 pointtype 7 pointsize 2
set output 'MoGdata_color.gif'           
plot 'MoGdata0' using 2:1, 'MoGdata1' using 2:1, 'MoGdata2' using 2:1, 'MoGcentroids' using 2:1 pointtype 7 pointsize 2
set output
set term x11

We get the following output:


kmeans:
iter: 1, norm: 1.49829818106
iter: 2, norm: 0.143890700032
iter: 3, norm: 0.0288887879026
iter: 4, norm: 0.00466759639044
iter: 5, norm: 0.00179109984505
iter: 6, norm: 0.0

em:
iter: 1, log-likelihood: -165991.331942
iter: 2, log-likelihood: -165986.870227

kcentroids:
[[ 4.96865803  3.83410557]
 [ 3.03335108  6.95460482]
 [ 6.9650538   6.94975746]]

MoGCentroids:
[[ 4.97013245  3.91434377]
 [ 3.02119885  6.98255096]
 [ 6.96649683  6.97264637]]

And also the following plots:

Contours of the bivariate gaussians sources of data to be clustered:

Data generated from the gaussians:

Data generated from the gaussians (no color):

3-D view of the gaussians:

Data as clustered by the k-means algorithm:

Data as clustered by the EM-mixture of gaussians algorithm (initialized with the results of k-means):