{
"metadata": {
"name": "Homework2"
},
"nbformat": 3,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"source": [
"Homework #2: Putting it all together (developing and fitting a model)"
]
},
{
"cell_type": "markdown",
"source": [
"Ok, so in this homework, you are going to put together the pieces introduced in the last few lectures and actually fit a couple models to some data. At the end, I hope you'll end up with some useful code you can use in your own research.",
"",
"If you remember from our lecture, the following figure provides a high-level \"overview\" of what we will accomplish in this homework. ",
"",
"",
"The basic idea is that we have our model, which is represented in the box labeled \"Cognitive Model\". The model takes as input a set of parameters (e.g., decay) and outputs predicted quantities (any number of predictions might come out of the model depending on what it is designed to do). Then, we want some way of relating the model predictions to the empirical data. Our preferred way of doing this is through a secondary probabilistic \"model\" which will give us a likelihood score for the data given the current model and parameters (for the reasons discussed in class and in the readinings). ",
"",
"Then, the whole system is wrapped up in a search algorithm (`fmin`) which tries to find the parameters which maximize the likelihood of the data (or or minimize the negative log of the likelihood scores).",
"",
"So far in class we've written code that does almost all of these steps. The main goal of this homework is to pull it all together."
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 0.1: Load libraries"
]
},
{
"cell_type": "markdown",
"source": [
"You are going to need to load various libraries in developing your code for this homework. **It is a good idea to put them up top in this file so when you come back to work on them it will be the first cell you execute.** Thus, this isn't really a \"step\" as much as a organizational technique you'll want for later. I've got you started by importing the `random()`, `randint()` and `normalvariate()` functions which generated pseudo-random numbers. Hard to lose points on this one!"
]
},
{
"cell_type": "heading",
"level": 3,
"source": [
"load some libs"
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"# load your libraries here using:",
"# import ",
"# or",
"# from import ",
"",
"from random import random, randint, normalvariate"
],
"language": "python",
"outputs": [],
"prompt_number": 125
},
{
"cell_type": "heading",
"level": 3,
"source": [
"define some constants/global variables"
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"# we'll also define some constants that will be useful later on",
"",
"# the human experiment data",
"humandata = array([0.31, 0.47, 0.51, 0.7, 0.72, 0.8]) # estimated from Figure 2.7 in the Lewandosky book",
"conditions = array([1.64, 1.92, 2.33, 2.64, 3.73, 3.94])"
],
"language": "python",
"outputs": [],
"prompt_number": 91
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 1.0: Create a Python function which includes your code for the phonological loop model."
]
},
{
"cell_type": "markdown",
"source": [
"The first step is to kind of \"wrap up\" our phonological loop code into a easy to use/manipulated function. Look back at homework #1 if you've forgotten how to write a function in Python. Some of you will have done this step already in class. ",
"",
"Basically, you just want to take the code you wrote and refactor it so the code is in a function. The arguments to the function should be the \"free\" parameters of the model.",
"",
"In particular, this is the design spec you are going for:",
"",
"1. make a function called `phonoloop()` that takes as arguments the mean and standard deviation of the population you are sampling the decay rate from.",
"2. the function should return a array of values which are the probability of correct recall for the items on the list.",
"",
"To help get you started I will give you the basic template for this function:"
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"def phonoloop(meandecay, stddecay):",
" # put your code here and return something other than '0'!",
" return 0"
],
"language": "python",
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"Of course, you don't want to literally `return 0`, but you want to return the vector of values we used in our code (it was probably called somethign like `pCor/nReps`).",
"",
"At the end you should be able to call phonoloop in a cell by itself like this:",
"",
"`phonoloop(0.5, 0.2)`",
"",
"and get something like this as output:",
"",
"`array([ 0.606, 0.664, 0.84 , 0.86 , 0.974, 0.988])`",
"",
"If you can do this you have created the three interconnected boxes shown here:",
"",
"",
"",
"",
"Nice job!"
]
},
{
"cell_type": "heading",
"level": 3,
"source": [
"Step 1.1"
]
},
{
"cell_type": "markdown",
"source": [
"Ok, so now I'd like you to make a small change to your code. The code we developed in class randomly draws a decay rate for each person using the `normalvariate` function. However, it is possible that the decay rate might go below zero (which is non-sensical... what is a negative decay rate?). Thus, please adjust your code to not allow decay rates smaller than zero or larger than 1.0. This can be done using a hard cutoff (although another approach would be to draw from something like the Beta distribution which is bounded between zero and one)."
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 2.0: Plot the output of the model compared with the human data"
]
},
{
"cell_type": "markdown",
"source": [
"Before we get deeper into the model, please make sure you have the necessary code to plot the human data along side you model. We defined a vector above (step 0.1) with the human data in it. For reference, here is the simple plotting function I stepped through in class:"
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"def simplelistplot(datax, datay, axrange, ylab, xlab, joined=False, data2x=None, data2y=None):",
" fig=plt.figure()",
" ax=fig.add_subplot(111)",
" if joined==True:",
" ax.plot(datax,datay,'ko-',antialiased=True,markersize=3,linewidth=1)",
" else:",
" ax.plot(datax,datay,'ko',antialiased=True,markersize=3)",
" if data2x!=None and data2y!=None:",
" ax.plot(data2x, data2y, 'ro', antialiased=True, markersize=5)",
" plt.axis(axrange)",
" plt.ylabel(ylab)",
" plt.xlabel(xlab)",
" plt.show()"
],
"language": "python",
"outputs": [],
"prompt_number": 93
},
{
"cell_type": "markdown",
"source": [
"Try to understand what this function is doing, then write some code which, for a given set of parameters, will plot your models predictions alongside the human data. Your plot should looks something like this:",
"",
""
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 3.0: Create some fake subject data"
]
},
{
"cell_type": "markdown",
"source": [
"While the Lewandowsky text gives us the average proportion correct for each speech rate in this hypothetical experiment, we don't have access to the individual data from subjects. However, we can make up some fake data for the purposes of exploring the consequences of fitting to individuals or to the group. ",
"",
"To do this we will have to make a few assumptions. First, lets assume that the participants in this particular experiment came into the lab and were tested for their memory with six different lists each have different length words. Thus, the design is within-subject (rather than each participant contributing a single data point). In addition, we will make the completely **unreasonable** assumption that people's performance on the different tasks is independent. This is probably not true in an experiment like this... some people have really good memory, others don't. However, for simplicity lets just assume that each individuals performance on list A is independent from their performance on list B.",
"",
"Given this we can generate some \"fake\" data in this experiment which conforms to the overall means reported in the textbook.",
"",
"Our basic assumption is that for each person, the proportion correct in the overall population determines the probability of correctly recalling one of the list words. So in the condition where the speech rate was 1.64 we know that the probability of recalling each of the 5 list items is 0.31. From this we can generate data reflecting the number of items a subject might remember by flipping a weighted coin:",
"",
"",
"",
"To help you out I've provided this example function which samples N binomial trials with probability `p`:"
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"def sample(p, ntimes):",
" mysum = 0.0",
" for i in range(int(ntimes)):",
" if random() < p:",
" mysum += 1",
" return mysum, ntimes"
],
"language": "python",
"outputs": [],
"prompt_number": 55
},
{
"cell_type": "markdown",
"source": [
"We'll do this for 30 subjects in all 6 conditions. This data will then be used later when we explore the consequences of fitting to individual data or to the group mean.",
"",
"Your goal in this step is to generate a matrix called `fakehumandata` that looks something like this:",
"",
"",
"",
"(except it should have 30 rows rather thatn 5). Each row is a subject in the experiment, each entry in a row is a tuple (set of two numbers) which are the number of items correctly remember on the list and the number of items on the list. So for example `(1.0, 5)` means 1 out of 5 items were correctly recalled by this subject.",
"",
"Next create a second matrix called `fakehumandatap` which show the proportions corresponding to the entries in `fakehumandata`. For example it would look like this:",
"",
"",
"",
"This is basically just the same matrix of data with the number of correct recalls divided by the number of items on the list. This format will be needed later for fitting using standard \"goodness of fit\" measures such as RMSE.",
"",
"Finally, let's make a \"group\" data vector that has the counts of the number of times each list item was remembered in the group as a whole. If there were 30 subjects in our sample and 5 items per list, this would mean 30*5 = 150 possibly recall events. Of these, we expect that `p` will be correctly remembered. This is kind of a group analog to our `fakehumandata` matrix. You can do this by simply summing the columns of the `fakehumandatamatrix`. The end result should look something like this:",
"",
"",
"",
"Let's call is `fakehumandatagroup`.",
"",
"All of this is good practice using python to solve an actual problem. If you run into problems feel free to email me or contact a friend in class."
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 4.0: Compute a measure of discrepancy or \"goodness of fit\""
]
},
{
"cell_type": "markdown",
"source": [
"Write a function that computes the root mean square error between the human data in the experiment (the probability of correct recall we copied from the textbook figures) and your model output (i.e., `phonoloop()`). We did this in class but it is good practice. Remember the RMSE is defined as:",
"",
"$$rmse = \\sqrt{ \\frac{\\sum_i^N{(x_i-p_i)^2}}{N} }$$",
"",
"where $x_i$ is the $i$th human data point and $p_i$ is your model's prediction for the $i$th data point. For sake of clarity lets write the function like",
"`rmse(humandata, modeldata)`."
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 5.0: Compute a measure of likelihood"
]
},
{
"cell_type": "markdown",
"source": [
"In class, we went over some code in Python which can compute a binomial likelihood. Step 3.0 should have convinced you that the binomial might be a nice probability model for relating the output of our model (probability of recall) to the human data. In this step, you should write a function similar to the one you developed in Step 4.0 but will compute the log likelihood of a vector of human data given a vector of probabilities given by your model. ",
"",
"Specficially, you want a function that you can call like this:",
"",
"`ll( fakehumandata[0], phonoloop(0.5, 0.1) )`",
"",
"and will return the log likelihood (`ll`) of the data in `fakehumandata[0]` (the first \"fake\" subject we generated) given the probabilities the model makes for a mean decay rate of 0.5 and a standard deviation of 0.1. The log function in python is... `log()`!",
"",
"Remember that to compute the likelihood of multiple data points you multiply the probabilities together. However, if you first take the log of the probability you ADD them together!.",
"",
"One thing to keep in mind... What happens when `p=1.0`? The likelihood of anything other than a perfect outcome (all list items remembered) will be 0.0. But what is the log(0.0)? -inf to python. Thus, it is typically to add a small constant (e.g., 0.00000000001) to all probabilities before taking the `log` so that if they go to zero they won't generate a numerical problem."
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 6.0: Adjust parameters by hand to find parameters that provide a good fit according to either measure"
]
},
{
"cell_type": "markdown",
"source": [
"Ok, congratulations! You have all the pieces you need now to fit you model. In fact, all the pieces in this figure are now in place:",
"",
"",
"",
"Before we try an automated search proceedure, try to play with the model parameters by hand to find the best fit. Perform the fit in two ways, both to the \"group\" data (rather than individuals). ",
"",
"First, try to find the best parameters to **minimize** the RMSE between the model predictions and the overall probabilities of correct recall for the group (the number we estimated from the figure in the book).",
"",
"Next, try to find the best parameters to **minimize** the **NEGATIVE log likelihood** of the model prediction and the aggregated frequency of correct recall data we generated in step 3.0.",
"",
"Are the results of the best fit you found the same or different between the RMSE and ML approach? ",
"",
"A couple hints: since decay rate is a random parameter in the model, we may have to run the model many times over and over again to get stable estimate. The code we developed in class did this. You may want to make sure you adjust your `phonoloop` function accordingly to ensure stable estimates of the model outputs. Otherwise, you'll tend to chase your own tail as the fit numbder change when you repeat the fit simple due to the randomness inherent in the model."
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 7.0: Use `fmin` to automatically search for good parameters"
]
},
{
"cell_type": "markdown",
"source": [
"Ok, maybe you did a good job at finding parameters, but it's still a hassle. Let's let an automatic search algorithm help us! This was always my favorite part of modeling because it's like you give a kind of complex problem solving task to the computer and it does a better job than you might be able to do! It's like AI!",
"",
"When you finish step 7.0 you will have completed this entire loop:",
"",
"",
"",
"But how do we do it? ",
"First, import the `fmin` (and perhaps `brute`) function from `scipy.optimize` (place this up in step 0.1 if you like to stay organized).",
"",
"`from scipy.optimize import fmin, brute`",
"",
"Next, we want to write a little \"wrapper function\" that will go around our `ll` or `rmse` function and place the arguments in the correct order for `fmin`. In particular, write a function called `fitmodel` that takes two arguments: `params` (a vector of the parameters we are optimizing) and `data` which will be the data we are trying to fit (in an assumed vector format in case there are multiple bits of data we are fitting). In particular,",
"",
"`fitmodelll(params, data)`",
"",
"should return the negative log likelihood of `data` using the vector of parameters in `params`. `params[0]` could be the mean decay rate and `params[1]` could be the standard deviation. Basically, the names of these variables are being lost.. it is just the order the appear in the `params` vector matters.",
"",
"Write a second function called `fitmodelrmse()` that works identically but instead calls the rmse code you wrote above."
]
},
{
"cell_type": "markdown",
"source": [
"Ok now you can use `fmin`! Here is the syntax:",
"",
"`res=fmin(evalf, iparams, args = [data], maxfun = 1000, maxiter = 500, full_output=1)`",
"",
"`evalf` is the evaluation function you are using (either `fitmodelll` or `fitmodelrmse`). ",
"",
"`iparams` is the initial parameters you start with. You might want to embed this in a `for`-loop and set the intial parameters randomly over and over (if the best fit parameters always end of the same, you can be more confident you are not in a local minima). ",
"",
"`args` is the additional arguments you want to provide to `evalf` (i.e., the second argument which we called `data`). It should be a list of additional arguments which is why i wrote `[data]`. This part can be a bit tricky sometimes. For example, when i tried to fit the `fakehumandatagroup` vector I had to set `args = [[data]]` otherwise python wasn't treating my vector of numbers as a entire piece of data. It's something to just be aware of.",
"",
"`maxfun`, `maxiter`, and `full_output` should be relatively self-explanatory given our class discussion but you can always google `scipy fmin` to learn more about the options to the function.",
"",
"If you hook everything up right, `fmin` should start searching the parameter space using the Nelder-Mead simplex algorithm trying to minimize either the RMSE or negative log likelihood (depending on which function you pass as `evalf`). Using the syntax I described above `res=fmin(...)`, at the end the variable `res` will contain all we want to know about the fit. In particular, a vector that looks like this",
"",
"`(array([ 0.8393479, 0.14426697]), 0.18626325456192419, 37, 122, 0)`",
"",
"The first `array` is the best fit parameters (in the order you interpretted them in `fitmodelll()` or `fitmodelrmse()`, next is the value of the best evaluation (i.e., the lowest RMSE or negative log likelihood score). Next is the total number of function evaluation, iterations, etc... (these are typically less important and can be ignored).",
"",
"You should repeat the fit many times (say 30? 100?) and record the best fit overall.",
"",
"Do the values that `fmin` finds match what you found by hand?",
"",
"After you find the best fit parameters, use the code you developed in step 2.0 to plot the resulting data. Make sure it looks reasonable!"
]
},
{
"cell_type": "heading",
"level": 1,
"source": [
"Step 8.0: Compare fits to individuals versus the group "
]
},
{
"cell_type": "markdown",
"source": [
"In the final step of this homework, I would like you to explore the consequences of fitting to individuals versus the group. We have data structures from step 3.0 that define various individual versus group data. There are basically four optimaization you can do",
"",
"- Fit to the group data using RMSE",
"- Fit to each individual using RMSE",
"- Fit to the group data using ML",
"- Fit to each individual using ML",
"",
"When I say \"fit to each individual\" i mean run fmin once for each person in the experiment and report the average best fit parameter values you found.",
"When I say \"fit to the group data\" I mean run fmin once to find the parameters that best capture the average of everyone.",
"",
"Summarize what you found. Do the parameter estimates from the individuals look like typical distribution (Gaussian?).",
"",
"To help you explore, here is some example code that will plot a histogram of a set of numbers:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# http://matplotlib.sourceforge.net/examples/api/histogram_demo.html",
"import numpy as np",
"import matplotlib.pyplot as plt",
"import matplotlib.mlab as mlab",
"",
"mu, sigma = 100, 15",
"x = mu + sigma * np.random.randn(10000)",
"",
"fit = plt.figure()",
"ax = fit.add_subplot(111)",
"ax.set_ylabel(\"probability\")",
"n, bins, patches = ax.hist(x, 50, normed=1, facecolor='green', alpha=0.75)",
"plt.show()"
],
"language": "python",
"outputs": [
{
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD9CAYAAACbSYGGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9UVHX+P/DnWIluebAJf+1B0EmOM5jFCMNFQ5x2iyFN\nodVOTZttip1Z1KTVIZddtrDOJqzbSlFHZzNsi6jTxjlfzRLCthFOGzNwPG4bP7JIyMoIIRESW8T3\n9w8O99PIjINwLzD4fJwz53DvvO/1/TqlT+59v+/7aoQQAkREREM0bqQ7QEREYwMDhYiIFMFAISIi\nRTBQiIhIEQwUIiJSBAOFiIgUoWqglJeXw2AwICIiAvn5+V7bZGZmQqfTITo6GvX19QCAc+fOQZIk\nREVFIS4uDjt37pTbd3R0IDk5GWFhYUhJSUFnZ6eaJRAR0QCpGijp6elwOBw4dOgQXnjhBZw6dcrj\ne7fbjYqKClRXV8Nut8NutwMAJkyYgA8++ABHjx7F4cOH8dJLL+Hzzz8HAOzatQthYWH47LPPEBoa\nit27d6tZAhERDZBqgdLe3g4ASEhIQHh4OBITE+FyuTzauFwurFq1ClqtFlarFXV1dfJ3P/vZzwAA\nnZ2dOH/+PIKCggD0hlBqaiqCgoKwdu3afuckIqKRoVqgVFVVQa/Xy9uRkZGorKz0aON2uxEZGSlv\nT5kyBQ0NDQCAnp4e3HLLLZg2bRo2btyImTNn9juvXq+H2+1WqwQiIroMV4/kHy6EwMUrv2g0GgDA\nVVddhf/85z9obGzE0qVLceutt8JoNPZr70vfeYiIaOCGshqXalcoJpNJHmQHgJqaGsTFxXm0kSQJ\ntbW18nZLSwt0Op1Hm1mzZmHp0qXylYjJZJJvjdXV1cFkMvnsQ19gjbXPE088MeJ9YH2sj/WNvc9Q\nqRYowcHBAHpnejU2NqKsrAySJHm0kSQJxcXFaG1tRVFREQwGAwDg1KlTOH36NACgtbUV7733Hlas\nWCEfU1BQgK6uLhQUFPQLKSIiGhmq3vLKy8uDzWZDd3c3Nm3ahJCQEDgcDgCAzWZDbGws4uPjERMT\nA61Wi8LCQgDAyZMn8Zvf/AY9PT2YPn067HY7ZsyYAQBIS0vDAw88gLlz52LBggXIzc1VswQiIhog\njVDiOmcU0mg0ilzCjUZOpxNms3mku6Ea1hfYWF/gGuq/mwwUIiICMPR/N7n0ChERKYKBQkREimCg\nEBGRIhgoRESkCAYKEREpgoFCRESKYKAQEZEiGChERKQIBgoRESmCgUJERIpgoBARkSIYKEREpAgG\nChERKYKBQkREimCgEBGRIhgoRESkCAYKEREpgoFCRESKuHqkO0A0lpgtZjS3NffbP007Dc5S5/B3\niGgYMVCIFNTc1ozpG6f32//t89+OQG+IhhdveRERkSIYKEREpAgGChERKYKBQkREimCgEBGRIjjL\ni+gy+ZoaDABNJ5owHf1neRFdCRgoRJfJ19RgAPjC/oXX/U2NTTCYDP32n/zqJGaEzui3n8+tUCBS\n9ZZXeXk5DAYDIiIikJ+f77VNZmYmdDodoqOjUV9fDwA4ceIEbrvtNsybNw9msxlFRUVy++zsbISG\nhsJoNMJoNKKkpETNEogU0YMeTN84vd+n63yX1/2+roCIRjNVr1DS09PhcDgQHh4Oi8UCq9WKkJAQ\n+Xu3242KigpUV1ejtLQUdrsdBw4cwDXXXIOdO3ciKioKp06dQmxsLFasWIHrrrsOGo0GmzdvxubN\nm9XsOhERXSbVrlDa29sBAAkJCQgPD0diYiJcLpdHG5fLhVWrVkGr1cJqtaKurg4AMH36dERFRQEA\nQkJCMG/ePFRVVcnHCSHU6jYREQ2SaoFSVVUFvV4vb0dGRqKystKjjdvtRmRkpLw9ZcoUNDQ0eLT5\n/PPPUVNTg9jYWHlffn4+4uLikJubi46ODpUqIBo5fWMu3j5mi3mku0fk1YgOygsh+l1taDQa+eeO\njg7ce++92LlzJ6699loAQFpaGh5//HGcOXMGGRkZcDgcsNvtXs+fnZ0t/2w2m2E2mxWvgUgNfWMu\n3nBdMFKK0+mE0+lU7HyqBYrJZEJGRoa8XVNTg6SkJI82kiShtrYWFosFANDS0gKdTgcA6O7uxsqV\nK7F69WokJyfLx0ydOhUAEBwcjA0bNmD9+vUDChSiscLXjDHODKPLdfEv2tu2bRvS+VS75RUcHAyg\nd6ZXY2MjysrKIEmSRxtJklBcXIzW1lYUFRXBYOj9SyKEQGpqKm666SY8+uijHsecPHkSAHD+/HkU\nFRVh6dKlapVANCr5mjHGmWE00lS95ZWXlwebzYbu7m5s2rQJISEhcDgcAACbzYbY2FjEx8cjJiYG\nWq0WhYWFAIAPP/wQhYWFuPnmm2E0GgEA27dvR1JSErZu3YqjR49i/PjxSEhIQFpampolEBHRAKka\nKEuWLJFnbvWx2Wwe2zk5OcjJyfHYFx8fjwsXLng95yuvvKJsJ4mISBFcy4uIiBTBQCEiIkUwUIiI\nSBFcHJLIB1+rCnNFYSLvGChEPvhaVdjXisJEVzre8iIiIkXwCoVojOAT9DTSGChEY4Sv9b+49hcN\nFwYKXfE4+E6kDAYKXfE4+E6kDA7KExGRIhgoRESkCAYKEREpgoFCRESKYKAQEZEiGChERKQIBgoR\nESmCgUJERIpgoBARkSIYKEREpAgGChERKYKBQkREimCgEBGRIrjaMF0RfC1RD3CZeiKlMFDoiuBr\niXqAy9QTKYW3vIiISBEMFCIiUgQDhYiIFKFqoJSXl8NgMCAiIgL5+fle22RmZkKn0yE6Ohr19fUA\ngBMnTuC2227DvHnzYDabUVRUJLfv6OhAcnIywsLCkJKSgs7OTjVLICKiAVI1UNLT0+FwOHDo0CG8\n8MILOHXqlMf3brcbFRUVqK6uht1uh91uBwBcc8012LlzJ2pqavDWW28hKytLDo5du3YhLCwMn332\nGUJDQ7F79241SyAiogFSLVDa29sBAAkJCQgPD0diYiJcLpdHG5fLhVWrVkGr1cJqtaKurg4AMH36\ndERFRQEAQkJCMG/ePFRVVQHoDaHU1FQEBQVh7dq1/c5JREQjQ7VAqaqqgl6vl7cjIyNRWVnp0cbt\ndiMyMlLenjJlChoaGjzafP7556ipqUFsbGy/8+r1erjdbrVKICKiyzCiz6EIISCE8Nin0Wjknzs6\nOnDvvfdi586duPbaa+VjBio7O1v+2Ww2w2w2D6m/RERjidPphNPpVOx8qgWKyWRCRkaGvF1TU4Ok\npCSPNpIkoba2FhaLBQDQ0tICnU4HAOju7sbKlSuxevVqJCcne5y3rq4ORqMRdXV1MJlMPvvw00Ah\nIiJPF/+ivW3btiGdT7VbXsHBwQB6Z3o1NjairKwMkiR5tJEkCcXFxWhtbUVRUREMBgOA3quQ1NRU\n3HTTTXj00Uf7HVNQUICuri4UFBQgLi5OrRKIiOgyqHrLKy8vDzabDd3d3di0aRNCQkLgcDgAADab\nDbGxsYiPj0dMTAy0Wi0KCwsBAB9++CEKCwtx8803w2g0AgC2b9+OpKQkpKWl4YEHHsDcuXOxYMEC\n5ObmqlkCERENkKqBsmTJEnnmVh+bzeaxnZOTg5ycHI998fHxuHDhgtdzTpo0Cfv27VO2o0RENGRc\nHJJojGtqbILBZPD63TTtNDhLncPbIRqzGChEY1wPenyutPzt898Oc29oLONaXkREpAgGChERKYKB\nQkREimCgEBGRIhgoRESkCL+B8qtf/QrvvPOOz+dCiIiIgAEESlpaGl577TXMmTMHv//97/Hpp58O\nR7+IiCjA+A2UO+64A0VFRThy5AhmzZqFX/7yl1i0aBFee+214egfEREFiAGNobS2tuLll1/Gnj17\nsGDBAmzatAmHDx9GSkqK2v0jIqIA4fdJ+bvvvhv19fVYvXo13n77bcyYMQMAcN99911y6XgiGv18\nLcvCJVloMPwGysMPP4ylS5d67Pvxxx8RFBQkv5aXiAKTr2VZuCQLDYbfW15//OMf++1buHChKp0h\nIqLA5fMK5eTJk/jmm2/Q1dWFI0eOQAgBjUaD7777DkFBQcPZRyIiCgA+A6W0tBT/+Mc/8PXXX2PL\nli3y/vDwcDz11FPD0jkiIgocPgPloYcewkMPPYTi4mKsXLlyOPtEREQByGegvPrqq1i9ejUaGxvx\nt7/9Td7fd+tr8+bNw9JBIiIKDD4D5ezZswCAjo4OaDQaeX9foBCNRmaLGc1tzf32N51ownR4f8kU\nESnDZ6D0vfs9Ozt7uPpCNGTNbc1ep8F+Yf9iBHpDdGXxGSiPPPKIz4M0Gg2ee+45VTpERESByWeg\nREdHQ6PRQAjR7zve8iIiootdcpYXERHRQPkMlPT0dDz77LNYvnx5v+80Gg3279+vaseIiCiw+AyU\nBx98EAA8Hmrsw1teRER0sUuOoQCA2WwGAHz99dfQaDT4+c9/PiwdIyKiwOJ3tWGXy4WHH35YHpwf\nN24cXnzxRcTGxqreOSIiChx+A2Xr1q3Ys2ePHCBVVVXIyMiA0+lUu29ERBRA/C5f39HRAYPh/17A\nYzAY0NHRMaCTl5eXw2AwICIiAvn5+V7bZGZmQqfTITo6GvX19fL+tWvXYtq0aZg/f75H++zsbISG\nhsJoNMJoNKKkpGRAfSEiInX5vEIpLi4G0DuGsmzZMtx9990QQmDfvn1YsmTJgE6enp4Oh8OB8PBw\nWCwWWK1WhISEyN+73W5UVFSguroapaWlsNvtOHDgAABgzZo1eOSRR+TJAX361hHjWmJERKOLz0B5\n++235dlcOp0OH3/8MQBg9uzZ+P777/2euL29HQCQkJAAAEhMTITL5cKyZcvkNi6XC6tWrYJWq4XV\nakVWVpb83eLFi9HY2Oj13N4etqQrC9fsIhp9fAbKyy+/PKQTV1VVQa/Xy9uRkZGorKz0CBS3243V\nq1fL21OmTEFDQwNuvPHGS547Pz8f//znP3H33Xdj/fr1mDRp0pD6SoGHa3YRjT5+B+X/97//4YMP\nPkBpaSm+//57+aqloKBgyH+4EKLf1Ya/Z1zS0tLw+OOP48yZM8jIyIDD4YDdbvfa9qcLW5rNZnkK\nNBERAU6nU9EJVn4DJSsrC2fPnsW7776L9PR0vP7667jtttv8nthkMiEjI0PerqmpQVJSkkcbSZJQ\nW1sLi8UCAGhpaYFOp7vkeadOnQoACA4OxoYNG7B+/foBBQoREXm6+Bftbdu2Del8fmd5vf/++8jP\nz8fEiRORnp6OgwcP4v333/d74uDgYAC9M70aGxtRVlYGSZI82kiShOLiYrS2tqKoqMhjNpkvJ0+e\nBACcP38eRUVFWLp0qd9jiOjyNDU2wWAy9PuYLeaR7hqNYn6vUK666ipoNBoYjUaUlZUhIiJCfvmW\nP3l5ebDZbOju7samTZsQEhICh8MBoPd9K7GxsYiPj0dMTAy0Wi0KCwvlY61WKw4fPozW1lbMnDkT\nTz75JNasWYOtW7fi6NGjGD9+PBISEpCWljbI0onIlx70eB2j+vb5b0egNxQo/AbKww8/jLa2Njz6\n6KOw2+345ptv8NRTTw3o5EuWLEFdXZ3Hvr4Xd/XJyclBTk5Ov2Nff/11r+d85ZVXBvRnExHR8BpQ\noACAVqvl0/FEROST30A5ffo0XnrpJfmJ9DvvvBOpqanyGAkREREwgEH57OxsNDU1Yfv27di+fTua\nmprwxBNPDEffiIgogPi9QikpKUFNTQ2uuuoqAIDRaMS8efNU7xgREQUWv1coK1euxHPPPYe2tja0\ntbXh+eefx8qVK4ejb0REFEB8XqFcd9118lPrP/zwg/zwoBAC1157Lf785z8PTw+JiCgg+AyUzs7O\n4ewHEREFOL9jKABQV1eH/fv3Q6PRYMWKFR6LPhIREQEDGEPZs2cPHnroIYwb19t0zZo12LNnj+od\nIyKiwOL3CmXv3r0oKSnB9ddfD6D3Qcdly5Zh3bp1qneOiIgCh98rlMmTJ6O1tVXebmtrw+TJk1Xt\nFBERBR6/VyibN29GUlKSvBJwfX29vMAjEV1Z+lYh9maadhqcpc7h7RCNKpcMlAsXLqC1tRXHjh1D\nZWUlNBoNJEmSx1OI1OTrNb8AX/U7UnytQgxwJWLyEyjjxo1DTk4O7rnnHixatGi4+kQEwPdrfgG+\n6pdoNPJ7qZGSkoLHHnsMn3zyify0fFtb23D0jYiIAojfMZSCggJoNBq89dZbHvuPHz+uWqeIiCjw\n+A2U+vp6vPPOOzh48CA0Gg2SkpJw1113DUffiIgogPgNlNzcXHz88cewWq0AgDfeeAOffPIJl7An\nIiIPfgPljTfewNGjRxEUFAQAWL58OaKiohgoRETkwe+g/KJFi/Duu+/K2wcPHsTChQtV7RQREQUe\nv1coH330Efbu3Su/8re9vR16vR7z58+HRqPBxx9/rHoniYho9PMbKD+9OiEiIvLFb6DMmjVrGLpB\nRESBjmuoEBGRIhgoRESkiAG9sZFITb4WgeQCkESBhYFCI87XIpBcAJIosDBQiEgRvt6VwvekXDlU\nHUMpLy+HwWBAREQE8vPzvbbJzMyETqdDdHQ06uvr5f1r167FtGnTMH/+fI/2HR0dSE5ORlhYGFJS\nUtDZ2almCUQ0QH3vSrn44+udNjT2qBoo6enpcDgcOHToEF544QWcOnXK43u3242KigpUV1fDbrfD\nbrfL361ZswYlJSX9zrlr1y6EhYXhs88+Q2hoKHbv3q1mCURENECqBUp7ezsAICEhAeHh4UhMTITL\n5fJo43K5sGrVKmi1WlitVtTV1cnfLV68GNdff32/87rdbqSmpiIoKAhr167td04iIhoZqgVKVVUV\n9Hq9vB0ZGYnKykqPNm63G5GRkfL2lClT0NDQMODz6vV6uN1uBXtNRESDNaKD8kIICCE89mk0Gr/H\nDFR2drb8s9lshtlsvpzuERGNaU6nE06nU7HzqRYoJpMJGRkZ8nZNTQ2SkpI82kiShNraWlgsFgBA\nS0sLdDqd3/PW1dXBaDSirq4OJpPJZ9ufBgoREXm6+Bftbdu2Del8qt3y6luduLy8HI2NjSgrK4Mk\nSR5tJElCcXExWltbUVRUBIOh/5TDi0mShIKCAnR1daGgoABxcXGq9J+IiC6PqrO88vLyYLPZcPvt\nt2P9+vUICQmBw+GAw+EAAMTGxiI+Ph4xMTF45plnsGPHDvlYq9WKRYsW4dixY5g5cyb27t0LAEhL\nS8OXX36JuXPn4uuvv8Zvf/tbNUsgIqIBUnUMZcmSJR4ztwDAZrN5bOfk5CAnJ6ffsa+//rrXc06a\nNAn79u1TrpNERKQILg5JRESKYKAQEZEiuJYXDRuuKkw0tjFQaNhwVWGisY2BQkSq4irEVw4GChGp\nqm8V4ot9+/y3I9AbUhMH5YmISBEMFCIiUgQDhYiIFMFAISIiRTBQiIhIEQwUIiJSBAOFiIgUwUAh\nIiJFMFCIiEgRDBQiIlIEA4WIiBTBQCEiIkUwUIiISBEMFCIiUgQDhYiIFMH3oRDRiPD14i2AL98K\nVAwUIhoRvl68BfDlW4GKt7yIiEgRDBQiIlIEb3mR4swWM5rbmvvtbzrRhOnwfouDiAIfA4UU19zW\n7PXe+Bf2L0agN0Q0XHjLi4iIFMFAISIiRagaKOXl5TAYDIiIiEB+fr7XNpmZmdDpdIiOjkZ9fb3f\nY7OzsxEaGgqj0Qij0YiSkhI1SyAfzBYzDCaD10/TiaaR7h4RjQBVx1DS09PhcDgQHh4Oi8UCq9WK\nkJAQ+Xu3242KigpUV1ejtLQUdrsdBw4c8Hrs/fffjxtuuAEajQabN2/G5s2b1ew6+eFrnATgWAnR\nlUq1K5T29nYAQEJCAsLDw5GYmAiXy+XRxuVyYdWqVdBqtbBarairq/N5bGVlpXycEEKtbhMR0SCp\nFihVVVXQ6/XydmRkpEcoAL1XKJGRkfL2lClT0NDQ4PfY/Px8xMXFITc3Fx0dHWqVQEREl2FEpw0L\nIfpdbWg0mksek5aWhscffxxnzpxBRkYGHA4H7Ha717bZ2dnyz2azGWazeahdJiIaM5xOJ5xOp2Ln\nUy1QTCYTMjIy5O2amhokJSV5tJEkCbW1tbBYLACAlpYW6HQ6aLVan8dOnToVABAcHIwNGzZg/fr1\nAwoUIiLydPEv2tu2bRvS+VS75RUcHAygd7ZWY2MjysrKIEmSRxtJklBcXIzW1lYUFRXBYOhdeXTy\n5Mk+jz158iQA4Pz58ygqKsLSpUvVKoGIiC6Dqre88vLyYLPZ0N3djU2bNiEkJAQOhwMAYLPZEBsb\ni/j4eMTExECr1aKwsPCSxwLA1q1bcfToUYwfPx4JCQlIS0tTswQiIhogVQNlyZIl8sytPjabzWM7\nJycHOTk5AzoWAF555RVlO0lEo46vd6XwPSmjG9fyIqJRx9e7UvielNGNS68QEZEiGChERKQI3vKi\nS+K7TWg04djK6MZAoUviu01oNOHYyujGW15ERKQIBgoRESmCgUJERIpgoBARkSIYKEREpAjO8iIA\nnB5MREPHQCEAnB5MgY3Pp4wODBQiCnh8PmV04BgKEREpgoFCRESKYKAQEZEiGChERKQIDspfQXxN\nDQY4PZiIho6BcgXxNTUY4PRgIho6BgoRjVm+nk8B+IyKGhgoYxCfeifq5ev5FIDPqKiBgTIG8al3\nIv/4dL3yGChEdEXi0/XK47RhIiJSBAOFiIgUwVteAYyD70TK49jK4DFQRjl/DyNKuVK//Rx8Jxo8\njq0MHgNllOPDiEQUKFQdQykvL4fBYEBERATy8/O9tsnMzIROp0N0dDTq6+v9HtvR0YHk5GSEhYUh\nJSUFnZ2dapYwKp2uPz3SXVAV6wtsY70+p9PZb5/ZYobBZOj3MVvMw96/kaTqFUp6ejocDgfCw8Nh\nsVhgtVoREhIif+92u1FRUYHq6mqUlpbCbrfjwIEDXo+9//77ccMNN2DXrl0ICwvDm2++iS1btmD3\n7t2w2+1qljHqnK4/jcn6ySPdDdWwvsA2VuvrG1s59c0phPw8xPM7H7efXXbXFTUeo1qgtLe3AwAS\nEhIAAImJiXC5XFi2bJncxuVyYdWqVdBqtbBarcjKyvJ5bGVlJZYtWwa3242srCwEBQVh7dq12L59\nu1olDNqlxj1OfnUSM0JnDHg/B9iJRoe+sZVz/+8cpqd4/p30dfv5ShuPUS1QqqqqoNfr5e3IyEg5\nFPq43W6sXr1a3p4yZQoaGhpw/Phxn8f+9Lx6vR5ut1uxPvsKAl+/TVxqlpW331aA3v/xfD3Fzqfb\nia4Mw7XG2OX+mzZkQiVlZWXivvvuk7d37dolsrKyPNr8+te/FiUlJfK2JEmioaHB67F/+tOfhBBC\nzJw5U3R1dQkhhPjhhx9EWFiY1z8fAD/88MMPP5f5GQrVrlBMJhMyMjLk7ZqaGiQlJXm0kSQJtbW1\nsFgsAICWlhbodDpotVqfx5pMJtTV1cFoNKKurg4mk8nrn9+bKURENFxUm+UVHBwMoHe2VmNjI8rK\nyiBJnreBJElCcXExWltbUVRUBIOh9xJw8uTJPo+VJAkFBQXo6upCQUEB4uLi1CqBiIgug6qzvPLy\n8mCz2dDd3Y1NmzYhJCQEDocDAGCz2RAbG4v4+HjExMRAq9WisLDwkscCQFpaGh544AHMnTsXCxYs\nQG5urpolEBHRQA3phtkocf78eREVFSXuuusuIYQQZ86cEStWrBAzZ84UycnJoqOjY4R7OHidnZ3i\nwQcfFBEREcJgMIjKysoxVd/f//53sXDhQrFgwQKRnp4uhAjs/35r1qwRU6dOFTfddJO871L1PPvs\ns2LOnDnCYDCIioqKkejyZfFWn91uF3q9XhiNRpGeni7Onj0rfxdI9Xmrrc9f//pXodFoRGtrq7wv\nkGoTwnd9BQUFQq/Xi8jISPHYY4/J+wdT35gIlGeeeUbcf//9Yvny5UIIIXJzc8XGjRvFuXPnxIYN\nG8SOHTtGuIeDt2XLFpGVlSW6urpEd3e3OH369Jipr7W1VcyaNUt0dnaKnp4eceedd4qSkpKArq+8\nvFwcOXLE4y+tr3qam5vF3LlzRVNTk3A6ncJoNI5UtwfMW33vvfee6OnpET09PWLdunViz549QojA\nq89bbUII8eWXXwqLxSJmzZolB0qg1SaE9/r++9//iri4OHHs2DEhhBDfffedEGLw9QX8asNfffUV\n3n33Xaxbt04eiHe73UhNTZWfVXG5XCPcy8E7dOgQ/vCHP2DChAm4+uqrERwcPGbqmzhxIoQQaG9v\nR1dXF86ePYvJkycHdH2LFy/G9ddf77HPVz0ulwtJSUkICwvDkiVLIIRAR0fHSHR7wLzVd8cdd2Dc\nuHEYN24cLBYLDh8+DCDw6vNWGwBs3rwZf/nLXzz2BVptgPf6Dh48iNTUVERERADofXQDGHx9AR8o\nv/vd77Bjxw6MG/d/paj5rMpw+uqrr3Du3DmkpaVBkiTk5uaiq6trzNQ3ceJE7Nq1C7NmzcL06dNx\n6623QpKkMVNfH1/1uFwueSIKAMydOzfga33xxRexfPlyAL1BGuj17du3D6Ghobj55ps99o+F2gDg\nvffewyeffIKYmBisW7cOtbW1AAZfX0AHyoEDBzB16lQYjUaPacJijEwZPnfuHI4dO4aVK1fC6XSi\npqYGb7755pipr6WlBWlpaaitrUVjYyM++ugjHDhwYMzU1+dy6tFoNCr2RF1PPvkkJk2ahHvuuQeA\n97oDqb6zZ8/i6aefxrZt2+R9fTUFem19zp07h7a2NlRUVCA5ORkbN24EMPj6AjpQ/v3vf2P//v2Y\nPXs2rFYr/vWvf2H16tXysyoALvmsymg3Z84czJ07F8uXL8fEiRNhtVpRUlIyZupzu92Ii4vDnDlz\ncMMNN+Cee+5BRUXFmKmvj696+p7D6lNfXx+wtb788ssoLS31mKkZ6PU1NDSgsbERt9xyC2bPno2v\nvvoK0dHRaG5uDvja+sTFxeHee+/FxIkTsXz5ctTX1+PcuXODri+gA+Xpp5/GiRMncPz4cbzxxhv4\nxS9+gVdMLSfsAAABeklEQVRffXVMPasSEREBl8uFCxcu4J133sHtt98+ZupbvHgxqqur0dbWhh9/\n/BEHDx5EYmLimKmvj696YmNjUVpaii+//BJOpxPjxo3DpEmTRri3l6+kpAQ7duzA/v37MWHCBHl/\noNc3f/58NDc34/jx4zh+/DhCQ0Nx5MgRTJs2LeBr67Nw4UIcPHgQQgi4XC7ceOONmDBhwuDrG/rc\ngdHB6XTKs7wCedrpxT799FMhSZK45ZZbxJYtW0RnZ+eYqm/v3r0iISFBxMTEiKysLNHT0xPQ9d13\n331ixowZYvz48SI0NFQUFBRcsp68vDxx4403CoPBIMrLy0ew5wPTV98111wjQkNDxUsvvSTmzJkj\nwsLCRFRUlIiKihJpaWly+0Cqz9t/u5+aPXu2x7ThQKpNCO/1nT9/XthsNqHX60VKSopwu91y+8HU\npxFijN2wJiKiERHQt7yIiGj0YKAQEZEiGChERKQIBgoRESmCgUJERIpgoBARkSL+P2Qt64KEDI7S\nAAAAAElFTkSuQmCC\n"
}
],
"prompt_number": 210
},
{
"cell_type": "code",
"collapsed": true,
"input": [],
"language": "python",
"outputs": []
}
]
}
]
}