Famafrench bayes

pyfolio/bayesian.py

@@ -19,6 +19,7 @@
 import scipy as sp
 from scipy import stats
 import seaborn as sns
+import theano.tensor as tt
 
 import matplotlib.pyplot as plt
 
@@ -34,14 +35,16 @@ def model_returns_t_alpha_beta(data, bmark, samples=2000):
     return sets. Usually, these will be algorithm returns and
     benchmark returns (e.g. S&P500). The data is assumed to be T
     distributed and thus is robust to outliers and takes tail events
-    into account.
+    into account.  If a pandas.DataFrame is passed as a benchmark, then
+    multiple linear regression is used to estimate alpha and beta.
 
     Parameters
     ----------
     returns : pandas.Series
         Series of simple returns of an algorithm or stock.
-    bmark : pandas.Series
-        Series of simple returns of a benchmark like the S&P500.
+    bmark : pandas.DataFrame
+        DataFrame of benchmark returns (e.g., S&P500) or risk factors (e.g.,
+        Fama-French SMB, HML, and UMD).
         If bmark has more recent returns than returns_train, these dates
         will be treated as missing values and predictions will be
         generated for them taking market correlations into account.
@@ -55,10 +58,15 @@ def model_returns_t_alpha_beta(data, bmark, samples=2000):
         of the posterior.
     """
 
-    if len(data) != len(bmark):
-        # pad missing data
+    bmark = bmark.dropna()
+
+    if data.shape[0] != bmark.shape[0]:
         data = pd.Series(data, index=bmark.index)
 
+    if bmark.ndim > 1:
+        Nbmark = bmark.shape[1]
+    else:
+        Nbmark = 1
     data_no_missing = data.dropna()
 
     with pm.Model():
@@ -69,16 +77,19 @@ def model_returns_t_alpha_beta(data, bmark, samples=2000):
         nu = pm.Exponential('nu_minus_two', 1. / 10., testval=.3)
 
         # alpha and beta
-        beta_init, alpha_init = sp.stats.linregress(
-            bmark.loc[data_no_missing.index],
-            data_no_missing)[:2]
-
-        alpha_reg = pm.Normal('alpha', mu=0, sd=.1, testval=alpha_init)
-        beta_reg = pm.Normal('beta', mu=0, sd=1, testval=beta_init)
-
+        X = bmark.loc[data_no_missing.index]
+        X['ones'] = np.ones(len(X))
+        y = data_no_missing
+        alphabeta_init = np.linalg.lstsq(X, y)[0]#[:2]
+
+        alpha_reg = pm.Normal('alpha', mu=0, sd=.1, testval=alphabeta_init[-1])
+        beta_reg = pm.Normal('beta', mu=0, sd=1, 
+                testval=alphabeta_init[:-1], shape=Nbmark)
+        bmark_theano = tt.as_tensor_variable(bmark.ix[data_no_missing.index].T)
+        mu_reg = alpha_reg + tt.dot(beta_reg, bmark_theano)
         pm.T('returns',
              nu=nu + 2,
-             mu=alpha_reg + beta_reg * bmark,
+             mu=mu_reg,
              sd=sigma,
              observed=data)
         start = pm.find_MAP(fmin=sp.optimize.fmin_powell)
@@ -528,7 +539,7 @@ def run_model(model, returns_train, returns_test=None,
         Out-of-sample returns. Datetimes in returns_test will be added to
         returns_train as missing values and predictions will be generated
         for them.
-    bmark : pd.Series (optional)
+    bmark : pd.Series or pd.DataFrame (optional)
         Only used for alpha_beta to estimate regression coefficients.
         If bmark has more recent returns than returns_train, these dates
         will be treated as missing values and predictions will be

pyfolio/examples/linear algebra tests.ipynb

@@ -0,0 +1,195 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "x = np.array([0, 1, 2, 3])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "y = np.array([-1, 0.2, 0.9, 2.1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "A = np.vstack([x, np.ones(len(x)), np.ones(len(x))]).T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(4, 3)"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 0.,  1.,  1.],\n",
+       "       [ 1.,  1.,  1.],\n",
+       "       [ 2.,  1.,  1.],\n",
+       "       [ 3.,  1.,  1.]])"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "mcc = np.linalg.lstsq(A, y)[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ 1.   , -0.475, -0.475])"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "mcc"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH6ZJREFUeJzt3XtwVPXdx/H3F4xILGAUAQmXyGCnVR4vVGhCFVdsHWsV\nb5HaKuoz9sGxajsOdtSqTUTE1pHWVm3Fjtbbg9RGR8FLBaMx2gYMVXxEBAoY0HhpEQrCogb4Pn9k\nicmy2Wz2bPb6ec1k3LPnt+f8Dke+OXz2d37H3B0RESkcvTLdARERSS8VfhGRAqPCLyJSYFT4RUQK\njAq/iEiBUeEXESkwgQq/mQ03s5fM7G0zW25mP4nRJmRmW8zsjcjPDUH2KSIiwewT8PMtwFXuvszM\nvgL8w8wWufs7Ue1edvfJAfclIiIpEOiK390/cvdlkdfbgHeAoTGaWpD9iIhI6qQs4zezMuAYYEnU\nKgcmmNmbZvasmR2eqn2KiEj3BY16AIjEPDXATyNX/u29Dgx397CZfRd4EvhqKvYrIiLdZ0Hn6jGz\nIuBp4Dl3vyOB9u8C33D3TVHva9IgEZEkuHu34vSgo3oMuA9Y0VnRN7PBkXaY2Xhaf9lsitXW3fP2\np6qqKuN90LHp+HR8+feTjKBRz7eAC4D/M7M3Iu/9HBgRKeRzgErgMjPbCYSB8wLuU0REAghU+N39\nVbr4V4O73w3cHWQ/IiKSOrpzN01CoVCmu9Bj8vnYQMeX6/L9+JIR+MvdVDEzz5a+iIjkCjPDu/nl\nbkqGc4pIZkXGT0ieS9XFsQq/SJ7Qv5jzWyp/uSvjFxEpMCr8IiIFRoVfRKTAKOMXyVPhcJj6+sXU\n1DSwYcMORozoS2VlBRMnllNcXJy2bSTr1ltvZd26dfzxj39Maduu9OrVizVr1jBq1Kgu21ZXV7N2\n7VoefvjhwPtNJw3nFMkDkSF9bctTpvyCxsY+NDdX0NJSDhQDYYqKFlNa2sC4cZ/z2GMz4m4zFdvY\n44EHHmD27NmsW7eO/v37c9ZZZ3HrrbcyYMCApI+5p3Sn8N90002sWbMmocJ/8cUXM3z4cG6++eak\n+hV9jqPeT99cPSKSfcLhMI2NfWhqup6Wlkm0FmyAYlpaJtHUdD2NjfsSDod7dBt7zJ49m2uvvZbZ\ns2ezdetWFi9ezPr16/nOd75DS0tLzM/s2rWrW8cs3aPCL5Jn6usX09xcEbdNc3MF9fWLe3QbAFu3\nbqW6upq77rqLk08+md69ezNy5Egee+wxmpqaeOSRR4DWyKSyspKpU6cyYMAAHnjgAaqrq5k6dWrb\nth566CFGjhzJwIEDmTlzJmVlZbz44ottn9/TtqmpiV69erW1P/jgg5k1a1bbdl577TUqKiooKSlh\n6NChXHnllZ3+Aor27rvvcsIJJ9C/f39OPvlkNm7c2GH9ueeeyyGHHMIBBxzACSecwIoVKwC49957\nmTt3Lrfddhv9+vXjjDPOAOCXv/wlo0ePpn///hxxxBE8+eSTCfUjKBV+kTxTU9MQiWY619JSQU1N\nQ49uA+Dvf/87n332GWeffXaH9/fff39OPfVUFi1a1Pbe/PnzOffcc9myZQvnn39+h3HrK1as4PLL\nL+fRRx/lww8/ZMuWLXzwwQdt62ONcf/b3/7G6tWrqa2tZcaMGaxatQqAffbZh9/+9rd88sknNDQ0\nUFtby+9///u4x7HHD3/4Q8aNG8cnn3zCjTfeyIMPPthh39/73vdYs2YN//73vxk7diznn38+ANOm\nTeP888/nmmuu4dNPP+Wpp54CYPTo0bz66qts3bqVqqoqLrjgAj766KOE+hKECr9IntmwYQdfRjOd\nKY6067ltAGzcuJGBAwfSq9fepWbIkCEdrpgnTJjA5Mmtj+beb7/9OuTZNTU1TJ48mQkTJlBUVMSM\nGTM6FNxY2XdVVRV9+vThyCOP5KijjmLZsmUAjB07lvHjx9OrVy9GjhzJtGnTePnll7s4VtiwYQNL\nly7l5ptvpqioiOOPP57TTz+9w74vvvhi9t9/f4qKiqiqquLNN9/k008/7bSflZWVDBkyBIApU6Zw\n2GGH8dprr3XZl6BU+EXyzIgRfWmdAT2ecKRdz20DYODAgWzcuJHdu3fvte7DDz/k4IMPblseNmxY\np9v54IMPOqzv27cvBx10UNx97ymoAMXFxWzfvh2A1atXc9ppp3HIIYcwYMAArr/+ej755JO429rT\nh5KSEvr2/fKYR44c2fZ6165dXHvttYwePZoBAwZw6KGHAuwVB7X30EMPccwxx1BSUkJJSQnLly9P\nqC9BqfCL5JnKygqKiuJn70VFDVRWdp7hp2IbABUVFfTp04fHH3+8w/vbtm3jr3/9KyeddFLbe/Gm\nJBg6dCjvv/9+2/KOHTuSLpCXXXYZhx9+OGvWrGHLli3ccsstMX8xRTvkkEPYvHlzhy+0169f39bv\nuXPnMn/+fGpra9myZQvvvvsu8OVVfvTxrV+/nmnTpnH33XezadMmNm/ezJgxY9Iy9YYKv0iemTix\nnNLS+Nl7aWkDEyd2nuGnYhsAAwYMoKqqiiuvvJLnn3+elpYWmpqamDJlCsOHD+/w5W0855xzDgsW\nLKChoYEvvviC6urqpAvktm3b6NevH8XFxaxcuZI//OEPCX1u5MiRHHvssVRVVdHS0sKrr77K008/\n3WG7ffr04cADD2T79u38/Oc/7/D5wYMHs27durbl7du3Y2YMHDiQ3bt386c//Ynly5cndUzdpcIv\nkmeKi4sZN+5zyspmUlRUy5eRTZiiolrKymYybtwXcW/ASsU29vjZz37GrFmzuPrqqxkwYADl5eWM\nHDmS2tpaioqKgNar4egr4vbvHXHEEdx5552cd955DB06lH79+jFo0CD69OkT8/Px/vVw++23M3fu\nXPr378+0adM477zzEv7s3LlzWbJkCQceeCAzZszgoosualt34YUXMnLkSEpLSxkzZgwVFRUdtnXJ\nJZewYsUKSkpKOPvsszn88MOZPn06FRUVDBkyhOXLl3Pcccd1+ecJrcNt//rXF/nRj25JqH003cAl\nkgdi3dyT63fuxrNt2zZKSkpYs2ZNh5w9n+05x3vfWLd/t2/gUuEXyQOd3dWZTxYsWMBJJ52EuzN9\n+nQaGxv5xz/+kelupY2ZsX37do444jc0NV3ffo3u3BWR/DR//nxKS0spLS1l7dq1zJs3L9NdSrtE\nbqxLhK74RfJAIVzxFzoz45JLZnLffVfR8R4LXfGLiOStxG6s65oKv4hIjkjsxrquBSr8ZjbczF4y\ns7fNbLmZ/aSTdr8zs3+a2ZtmdkyQfYqIFKroG+vG8FZS2wl6xd8CXOXuRwDlwOVm9vX2DczsVGC0\nux8GTAMSu1tCREQ62HNj3VdZxQJO4wnO7vpDMQQq/O7+kbsvi7zeBrwDDI1qNhl4MNJmCXCAmQ0O\nsl8RkUJU/Pnn3FW0gIZeY3mlVyljSG5Ct5Rl/GZWBhwDLIlaVQq81275faDz2ZhEpCD069ePpqam\nHtt+9Hz+7dXV1TF8+PC25TFjxlBfX99jfUmZr32N7006mv3Wvs2Jz3yfqZckNp10tJQ8c9fMvgLU\nAD+NXPnv1SRqOea4s+rq6rbXoVCIUCiUiu6JSAaVlZXxr3/9i969ewOtwxJXr17dYbriWI8lLCsr\n4/7772fSpElJ7Tfe1AvR0jVHTmCLFlG3aRN1DzwAQJwJTeMKXPjNrAh4HHjE3WM9PqYZGN5ueVjk\nvb20L/wikh/MjKeffrrbBVz3JsRw5JGEoMNF8U033dTtzQQd1WPAfcAKd7+jk2bzgQsj7cuB/7j7\nx0H2KyK5r1evXqxdu3avxxJOnjyZCy+8kA0bNnD66afTr18/br/9dgAWL17MhAkTKCkp4eijj+7w\nAJWuHosYT/RjHKdMmcJFF11E//79GTNmTIepIT744APOOeccBg0axKhRo7jzzjtT9CeSRu6e9A9w\nHLAbWAa8Efn5LnApcGm7dncBa4A3gbGdbMtFJDnZ/PenrKzMX3jhhb3eNzNfu3atu7tffPHFfuON\nN+71udra2rbl999/3w866CB/7rnn3N190aJFftBBB/nGjRvd3b28vNynT5/uX3zxhdfX13u/fv18\n6tSpMfv00ksv+bBhw2Luq6qqyvfbbz9/7rnnfPfu3X7dddd5eXm5u7vv2rXLx44d6zfffLO3tLT4\nunXrfNSoUf78888n+8eTsM7OceT9btXuQFGPu79KAv9qcPcrguxHRALqRt4dVxLRi7tz5plnss8+\nreXmxBNP5IknnojZLp5HHnmEU089lVNOOQWAb3/72xx77LE888wzhEIhli5dyosvvtjpYxG74/jj\nj2/bzwUXXMAdd7QGGo2NjWzcuJEbbrgBgEMPPZQf/ehHzJs3j5NPPjmpfWVCSr7cFZEsl8Gs3Mx4\n6qmnkv6Sdo/169fzl7/8hQULFrS9t3PnTiZNmtTpYxHfe++9WJvq0uDBX444Ly4u5rPPPmP37t2s\nX7++bV977Nq1i4kTJya1n0xR4ReRjIs1Aif6vREjRjB16lTuvffevdquX7++7bGIe54TsH79+raR\nRKkyfPhwDj30UFavXp3S7aab5uoRkYyLfizhnvfWrl3btnzBBRewYMECFi5cyK5du/jss8+oq6uj\nubm5y8cipsr48ePp168ft912Gzt27GDXrl0sX76cpUuXpnxfPUmFX0QyIt5jCQGuu+46Zs6cSUlJ\nCb/+9a8ZNmwYTz31FLNmzWLQoEGMGDGC2bNntz0oPd5jEbvaf/T7sR4DCdC7d2+efvppli1bxqhR\nozj44IOZNm0aW7duTfrPIRM0H79IHtCY9/zX2TmOvK/5+EVEpHMq/CIiBUaFX0SkwKjwi4gUGBV+\nEZECo8IvIlJgdOeuSJ7ozvzzUthU+EXygMbwp8GqVTB9OqxcCbNnw+TJqZv8Ls0U9YiIxLN5M1x1\nFXzrWxAKwdtvwxln5GzRBxV+EZHYdu6Eu++Gr30NduyAFSvg6quhT59M9ywwRT0iItEWLmy9yh88\nGBYtgiOPzHSPUkqFX0RkjzzK8eNR1CMikoc5fjwq/CJSuPI4x49HUY+IFKY8z/HjUeEXkcJSIDl+\nPIp6RKQwFFiOH48Kv4jktwLN8eMJHPWY2f3A94B/uft/xVgfAp4C9jxJ+XF3nxl0vyIiXSrgHD+e\nVGT8fwLuBB6K0+Zld5+cgn2JiHRNOX5cgaMed38F2NxFM/2Ji0jPU46fkHRk/A5MMLM3zexZMzs8\nDfsUkUKiHL9b0jGc83VguLuHzey7wJPAV2M1rK6ubnsdCoUIhUJp6J6I5LQCy/Hr6uqoq6sLtA1L\nxTzeZlYGLIj15W6Mtu8C33D3TVHvu+YUF5GEKccHWh/A4+7dOvAej3rMbLBFHg1kZuNp/WWzqYuP\niYjEphw/sFQM53wUOAEYaGbvAVVAEYC7zwEqgcvMbCcQBs4Luk8RKUA7d8KcOTBjBpx1VmuOP2hQ\npnuVk1IS9aSCoh4R6VT7HP+OO/I+x++OZKIezdUjItlLOX6P0JQNIpJ9lOP3KBV+EckeGo+fFop6\nRCQ7FNh4/ExS4ReRzFKOn3aKekQkM5TjZ4wKv4ikl3L8jFPUIyLpoxw/K6jwi0jPU46fVRT1iEjP\nUY6flVT4RST1lONnNUU9IpJayvGzngq/iKSGcvycoahHRIJRjp9zVPhFJDnK8XOWoh4R6T7l+DlN\nhV9EEqccPy8o6hGRrinHzysq/CLSOeX4eUlRj4jEphw/b6nwi0hHyvHznqIeEWmlHL9gqPCLFDrl\n+AUncOE3s/vN7GMzeytOm9+Z2T/N7E0zOyboPkUkRRYuhKOOgscfb83x77kHBg3KdK+kh6Xiiv9P\nwCmdrTSzU4HR7n4YMA34Qwr2KSJBrFoFp50GP/4x3HIL1Nbqy9sCErjwu/srwOY4TSYDD0baLgEO\nMLPBQfcrIkmIleOfeaZy/AKTjoy/FHiv3fL7wLA07FdE9lCOL+2kazhn9OWEx2pUXV3d9joUChEK\nhXquRyKFQuPx80pdXR11dXWBtmHuMWtw9zZiVgYscPf/irHuHqDO3edFllcCJ7j7x1HtPBV9EZEI\njccvCGaGu3frxKbjin8+cAUwz8zKgf9EF30RSaHNm2HGDHj4Ybj22tYRO1GRTjgcpr5+MTU1DWzY\nsIMRI/pSWVnBxInlFBcXZ6jjki6Br/jN7FHgBGAg8DFQBRQBuPucSJu7aB35sx34b3d/PcZ2dMUv\nEsTOnTBnTmvRP+us1v/GGJo5ZcovaGzsQ3NzBS0t5UAxEKaoaDGlpQ2MG/c5jz02I+3dl+Rk5Irf\n3X+QQJsrgu5HROJIMMcPh8M0Nvahqen6qDXFtLRMoqlpEjCTcDisK/88pjt3RXJZ+/H4s2Z1OR6/\nvn4xzc0VcTfZ3FxBff3iVPdUsogKv0guSnJenZqahki807mWlgpqahpS2FnJNir8Irkk4Hj8DRt2\n0Jrpx1McaSf5StMyi+SKFIzHHzGiLxAmfvEPR9pJvtIVv0i262aOH09lZQVFRfHz+6KiBior438P\nILlNhV8kW/XA/PgTJ5ZTWho/vy8tbWDixPjfA0huU+EXyTY9OK9OcXEx48Z9TlnZTIqKammNfaB1\nHH8tZWUzGTfuCw3lzHMpmbIhFXQDlwgdc/w77uixeXV0527+SOYGLhV+kWygeXUkSckUfkU9Ipmk\n59xKBqjwi2SC5seXDNI4fpF00/z4kmEq/CLpohxfsoSiHpGephxfsowKv0hPUY4vWUpRj0hPUI4v\nWUyFXySVlONLDlDUI5IKyvElh6jwiwShHF9ykKIekWQpx5ccpcIv0l3K8SXHKeoRSZRyfMkTKvwi\nXVGOL3kmcOE3s1PMbKWZ/dPMromxPmRmW8zsjcjPDUH3KZI2CxfCUUfB44+35vj33AODBmW6VyKB\nBMr4zaw3cBfwbaAZaDSz+e7+TlTTl919cpB9iaSVcnzJY0Gv+McDa9y9yd1bgHnAGTHa6W+M5Abl\n+FIAghb+UuC9dsvvR95rz4EJZvammT1rZocH3KdI6inHlwISdDhnIs9KfB0Y7u5hM/su8CTw1VgN\nq6ur216HQiFCoVDA7okkQOPxJYfU1dVRV1cXaBuBnrlrZuVAtbufElm+Dtjt7r+K85l3gW+4+6ao\n9/XMXUkv5fiSBzLxzN2lwGFmVmZm+wLfB+ZHdWqwWevfJjMbT+svm017b0okTZTjS4ELFPW4+04z\nuwJ4HugN3Ofu75jZpZH1c4BK4DIz2wmEgfMC9lkkOTt3wpw5MGMGnHVWa46voZlSgAJFPamkqEd6\n1J4cf8gQ+M1vlONL3kgm6tFcPZLf9uT4q1bB7bcrxxdBUzZIvtqT4x93HJx4IixfrhxfJEKFX/JL\n9Hj8t99uveLXeHyRNop6JH+0z/E1Hl+kUyr8kvuU44t0i6IeyV3K8UWSosIvuUc5vkgginoktyjH\nFwlMhV9yg3J8kZRR4ZekhMNh6usXU1PTwIYNOxgxoi+VlRVMnFhOcXFx6na0eXPrFAuPPALXXtv6\nJCxFOiKBaMoG6bYpU35BY2MfmpsraGkpB4qBMEVFiyktbWDcuM957LEZwXYSPa/OjBmaV0ckBk3Z\nID0uHA7T2NiHpqbro9YU09IyiaamScBMwuFw8lf+yvFFepRG9Ui31Ncvprm5Im6b5uYK6usXd3/j\nq1bBaafB5ZfDrFnwwgsq+iI9QIVfuqWmpiES73SupaWCmpqGxDeq8fgiaaXCL92yYcMOWjP9eIoj\n7bqg8fgiGaGMX7plxIi+tD5PJ17xD0faxaEcXyRjdMUv3VJZWUFRUfz8vqiogcrKTr4HUI4vknEq\n/NItEyeWU1oaP78vLW1g4sSo7wGU44tkDRV+6Zbi4mLGjfucsrKZFBXV0hr7QOs4/lrKymYybtwX\nXw7lVI4vknV0A5ckJaE7d/WcW5Eel8wNXCr8knqaV0ckbZIp/Ip6JHWU44vkBBV+CU45vkhOCVz4\nzewUM1tpZv80s2s6afO7yPo3zeyYoPuULLJwIRx1FDzxROt4/Hvu0WRqIlku0A1cZtYbuAv4NtAM\nNJrZfHd/p12bU4HR7n6YmX0T+AMQ/55/yX6rVsHVV8PKlcrxRXJM0Cv+8cAad29y9xZgHnBGVJvJ\nwIMA7r4EOMDMBgfcr2RK+xw/FFKOL5KDghb+UuC9dsvvR97rqs2wgPuVdNu5E37/e+X4Inkg6Fw9\niY6/jL4cjPm56urqttehUIhQKJRUpyTFNK+OSNaoq6ujrq4u0DYCjeM3s3Kg2t1PiSxfB+x291+1\na3MPUOfu8yLLK4ET3P3jqG1pHH+2UY4vkvUyMY5/KXCYmZWZ2b7A94H5UW3mAxdGOlgO/Ce66EuW\nUY4vktcCFX533wlcATwPrAD+7O7vmNmlZnZppM2zwDozWwPMAX4csM/SU5TjixQETdkgrTSvjkhO\n0sPWpfuU44sUHE3ZUKiU44sULBX+QqMcX6TgKeopJBqPLyKo8BcG5fgi0o6innymHF9EYlDhz0fK\n8UUkDkU9+UY5voh0QYU/XyjHF5EEKerJdcrxRaSbVPhzlXJ8EUmSop5cpBxfRAJQ4c8lyvFFJAUU\n9eQC5fgikkIq/NlMOb6I9ABFPdlKOb6I9BAV/myjHF9EepiinmyhHF9E0kSFP9OU44tIminqySTl\n+CKSASr8maAcX0QySFFPOinHF5EsoMKfDsrxRSSLJB31mNmBwJ+BkUATMMXd/xOjXROwFdgFtLj7\n+GT3mZOU44tIljF3T+6DZrcBG939NjO7Bihx92tjtHsX+Ia7b+pie55sX7KScnwRSQMzw927VVyC\nRD2TgQcjrx8EzozTtnAqnnJ8EclyQQr/YHf/OPL6Y2BwJ+0ceMHMlprZ/wTYX3ZTji8iOSJuxm9m\ni4AhMVZd337B3d3MOstpvuXuH5rZwcAiM1vp7q/EalhdXd32OhQKEQqF4nUveyjHF5E0qauro66u\nLtA2gmT8K4GQu39kZocAL7n717r4TBWwzd1nx1iXexm/cnwRybB0Z/zzgYsiry8CnozRoWIz6xd5\nvT9wMvBWgH1mB+X4IpLDghT+XwLfMbPVwKTIMmY21MyeibQZArxiZsuAJcDT7r4wSIczSjm+iOSB\npKOeVMv6qKd9jv+b3yjHF5GskEzUo7l6uqIcX0TyjKZs6IxyfBHJUyr80ZTji0ieU9TTnsbji0gB\nUOEH5fgiUlAKO+pRji8iBagwC79yfBEpYIUX9SjHF5ECVziFXzm+iAhQCFGPcnwRkQ7yt/ArxxcR\niSk/ox7l+CIincqvwq8cX0SkS/kR9SjHFxFJWG4XfuX4IiLdlrtRj3J8EZGk5F7hV44vIhJI7kQ9\nyvFFRFIi+wu/cnwRkZTK7qhHOb6ISMplZ+FXji8i0mOyK+pRji8i0uOSLvxmdq6ZvW1mu8xsbJx2\np5jZSjP7p5ldE2+bW4cOZ8PqNYQbG5Xji4j0kCBX/G8BZwH1nTUws97AXcApwOHAD8zs6521P+6z\nWkYvuoojTvxfpkz5RYCuZZ+6urpMd6HH5POxgY4v1+X78SUj6cLv7ivdfXUXzcYDa9y9yd1bgHnA\nGZ01fotv0tIyiaam62ls3JdwOJxs97JOPv/Pl8/HBjq+XJfvx5eMns74S4H32i2/H3mvS83NFdTX\nL+6RTomIFLK4hd/MFpnZWzF+Tk9w+55sx1paKqipaUj24yIi0glzT7o2t27A7CVguru/HmNdOVDt\n7qdElq8Ddrv7r2K0DdYREZEC5e7dGvqYqnH8ne10KXCYmZUBHwDfB34Qq2F3Oy4iIskJMpzzLDN7\nDygHnjGz5yLvDzWzZwDcfSdwBfA8sAL4s7u/E7zbIiKSrMBRj4iI5Ja03rmbyM1cZva7yPo3zeyY\ndPYvqK6Oz8xCZrbFzN6I/NyQiX4mw8zuN7OPzeytOG1y+dzFPb5cPncAZjbczF6K3HS53Mx+0km7\nnDyHiRxfrp5DM9vPzJaY2TIzW2Fmt3bSLvFz5+5p+QF6A2uAMqAIWAZ8ParNqcCzkdffBBanq39p\nOr4QMD/TfU3y+I4HjgHe6mR9zp67BI8vZ89dpP9DgKMjr78CrMqzv3+JHF/OnkOgOPLffYDFwHFB\nzl06r/gTuZlrMvAggLsvAQ4ws8Fp7GMQid6slpNfYrv7K8DmOE1y+dwlcnyQo+cOwN0/cvdlkdfb\ngHeAoVHNcvYcJnh8kKPn0N333M26L60XmZuimnTr3KWz8CdyM1esNsN6uF+pksjxOTAh8k+xZ83s\n8LT1rufl8rlLRN6cu8gou2OAJVGr8uIcxjm+nD2HZtbLzJYBHwMvufuKqCbdOnfpnJY50W+Ro38j\n58q3z4n083VguLuHzey7wJPAV3u2W2mVq+cuEXlx7szsK0AN8NPIlfFeTaKWc+ocdnF8OXsO3X03\ncLSZDQCeN7OQu9dFNUv43KXzir8ZGN5ueTitv5XitRkWeS8XdHl87v7pnn+yuftzQJGZHZi+Lvao\nXD53XcqHc2dmRcDjwCPu/mSMJjl9Drs6vnw4h+6+BXgGODZqVbfOXToLf9vNXGa2L603c82PajMf\nuBDa7vr9j7t/nMY+BtHl8ZnZYLPWhwuY2Xhah9NGZ3W5KpfPXZdy/dxF+n4fsMLd7+ikWc6ew0SO\nL1fPoZkNNLMDIq/7At8B3ohq1q1zl7aox913mtmem7l6A/e5+ztmdmlk/Rx3f9bMTjWzNcB24L/T\n1b+gEjk+oBK4zMx2AmHgvIx1uJvM7FHgBGBg5Ma9KlpHL+X8uYOuj48cPncR3wIuAP7PzPYUjZ8D\nIyAvzmGXx0funsNDgAfNrBetF+sPu3ttkNqpG7hERApMdj16UUREepwKv4hIgVHhFxEpMCr8IiIF\nRoVfRKTAqPCLiBQYFX4RkQKjwi8iUmD+H7zWL1hRa4t/AAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x1062b1710>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(x, y, 'o', label='Original data', markersize=10)\n",
+    "plt.plot(x, m*x + c, 'r', label='Fitted line')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 2",
+   "language": "python",
+   "name": "python2"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}