Skip to content

Instantly share code, notes, and snippets.

@aguschin
Created December 1, 2016 06:44
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save aguschin/08594da3c71c90fe02b46ced5c234900 to your computer and use it in GitHub Desktop.
Save aguschin/08594da3c71c90fe02b46ced5c234900 to your computer and use it in GitHub Desktop.
hw1
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"from sklearn.datasets import fetch_covtype\n",
"from sklearn.preprocessing import StandardScaler\n",
"import scipy\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"import seaborn as sns"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"{'DESCR': 'Forest covertype dataset.\\n\\nA classic dataset for classification benchmarks, featuring categorical and\\nreal-valued features.\\n\\nThe dataset page is available from UCI Machine Learning Repository\\n\\n http://archive.ics.uci.edu/ml/datasets/Covertype\\n\\nCourtesy of Jock A. Blackard and Colorado State University.\\n',\n",
" 'data': array([[ 2.59600000e+03, 5.10000000e+01, 3.00000000e+00, ...,\n",
" 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
" [ 2.59000000e+03, 5.60000000e+01, 2.00000000e+00, ...,\n",
" 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
" [ 2.80400000e+03, 1.39000000e+02, 9.00000000e+00, ...,\n",
" 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
" ..., \n",
" [ 2.38600000e+03, 1.59000000e+02, 1.70000000e+01, ...,\n",
" 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
" [ 2.38400000e+03, 1.70000000e+02, 1.50000000e+01, ...,\n",
" 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
" [ 2.38300000e+03, 1.65000000e+02, 1.30000000e+01, ...,\n",
" 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]),\n",
" 'target': array([5, 5, 2, ..., 3, 3, 3], dtype=int32)}"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = fetch_covtype()\n",
"data"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X, y = data['data'][:100], data['target'][:100] == 2\n",
"X = StandardScaler().fit_transform(X)"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)"
]
},
"execution_count": 103,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sklearn.linear_model import LinearRegression\n",
"lr = LinearRegression()\n",
"lr.fit(X[:,:1], y)"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"w = np.array([lr.intercept_, lr.coef_])"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.2379, 0. , 0. , ..., 0. , 0. , 0. ],\n",
" [ 0. , 0.2379, 0. , ..., 0. , 0. , 0. ],\n",
" [ 0. , 0. , 0.2379, ..., 0. , 0. , 0. ],\n",
" ..., \n",
" [ 0. , 0. , 0. , ..., 0.2379, 0. , 0. ],\n",
" [ 0. , 0. , 0. , ..., 0. , 0.2379, 0. ],\n",
" [ 0. , 0. , 0. , ..., 0. , 0. , 0.2379]])"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cov_data = np.eye(len(y)) * y.std() ** 2\n",
"cov_data"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ypred = lr.predict(X[:,:1])"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def data_likelyhood(y, ypred, cov_data):\n",
" err = (y - ypred).reshape(-1, 1)\n",
" return 1e8 * np.exp(err.T.dot(np.linalg.inv(cov_data)).dot(err) / (-2.))[0][0]"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3.0514565633470453e-08"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_likelyhood(y, ypred, cov_data)"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def calc_error(y, ypred, cov_data):\n",
" err = np.asarray(y - ypred).reshape(-1, 1)\n",
" return err.T.dot(np.linalg.inv(cov_data)).dot(err)"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"width0 = 10\n",
"width1 = 10\n",
"num = 20\n",
"w0_array = np.linspace(w[0] - width0, w[0] + width0, num)\n",
"w1_array = np.linspace(w[1] - width1, w[1] + width1, num)\n",
"likelyhood_array = np.zeros([num, num])\n",
"for (i, w0_) in enumerate(w0_array):\n",
" for (j, w1_) in enumerate(w1_array):\n",
" pred_ = np.asarray(w0_ + w1_ * X[:,:1])[:, 0]\n",
" likelyhood_array[i, j] = calc_error(y, pred_, cov_data)"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([ -9.61 , -8.55736842, -7.50473684, -6.45210526,\n",
" -5.39947368, -4.34684211, -3.29421053, -2.24157895,\n",
" -1.18894737, -0.13631579, 0.91631579, 1.96894737,\n",
" 3.02157895, 4.07421053, 5.12684211, 6.17947368,\n",
" 7.23210526, 8.28473684, 9.33736842, 10.39 ])"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w0_array"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7f88b2c1b710>"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAw0AAAKXCAYAAADAa6ZxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X2QpVd9H/jv7Z43zegFJAFCYJmYxMdFEQwSb4mTdRBO\nUQYUL7XK2uGlAixQWhxAgmWBVYJceEvgEshQIUgCYolljQp2pVphil3C8uI1uwZrIl68pnwcRMlg\ntIoG9DKS5q1nuvePe6eqPTVzpntmnj79TH8+qi5133Pv/Z253X37/u73nOeZLC0tBQAA4Fjmek8A\nAABY3zQNAABAk6YBAABo0jQAAABNmgYAAKBJ0wAAADRpGgAAgCZNAwAA0LSp9wQOe9bP/2qXs8w9\n9ZwLe5RNkjzlrPO61b7g7LO71T53x/Zutc/ZvqVb7R1nbO5We9vWfr/qmzf3e29i0/zGfF/k4KHF\nbrUXFvrV3rf/YLfaj+1d6Fb74T0HutV+4LE93Wrft3t3t9o/eeRn3Wr/zcP3dqudJN/76z+edJ3A\nMfR6HZms38fkVNuYf1EBAIAV0zQAAABN62Z5EgAAnIjJZEOsEOpK0gAAADRpGgAAgCbLkwAAGLXJ\nxPvgQ/MIAwAATZoGAACgSdMAAAA02dMAAMCozcUhV4cmaQAAAJo0DQAAQNMJLU8qpexIcv7sy121\n1j2nbkoAALByzgg9vFU1DaWUK5O8KUlZdvFSKeX7ST5Wa73xVE4OAADob8VNQynlA0lekeT6JHcl\neWA2dF6S5yd5ZynlCbXW3z3lswQAgGOYc3K3wa0mafjNJL9Wa737iMvvTvJnpZT/M8lXk2gaAADg\nNLKatuysJP+5Mf6TJOec3HQAAID1ZjVNwzeTXFdKOfvIgVLKuUk+lOTrp2heAACwIpPJpNvHRrGa\n5UlvTnJ7kp+WUu5J8mCSSaZ7Gi5KcmeSy0/x/AAAgM5W3DTUWn+U5LmllOcmuTjTZiFJdiXZWWv9\nzgDzAwAAOlv1eRpqrTuT7BxgLgAAwDrk+FQAAEDTCZ0RGgAA1otJNs6G5F4kDQAAQJOkAQCAUXNG\n6OF5hAEAgCZJAwAAo7aRTrLWi6QBAABo0jQAAABNlicBADBqc5YnDU7SAAAANK2bpOGp51zYpe7P\nP+6JXeomyYVnn92t9hPP3tGt9uPO2tqt9llnbulWe/v2zd1qb9nW71d906Z+703MbdqY7zwtHlzq\nVvvgwcVutQ/sO9it9p49C91q73i033PL9i39nlu2zM93q72pY23oRdIAAAA0aRoAAICmdbM8CQAA\nTsTE++CD8wgDAABNkgYAAEbNGaGHJ2kAAACaNA0AAECT5UkAAIyaM0IPT9IAAAA0SRoAABi1SSQN\nQ5M0AAAATZoGAACgSdMAAAA0aRoAAIAmG6EBABi1uYn3wYfmEQYAAJokDQAAjNrEyd0GJ2kAAACa\nNA0AAECT5UkAAIzanOVJg5M0AAAATZIGAABGbRJJw9AkDQAAQJOmAQAAaNI0AAAATZoGAACgyUZo\nAABGbW7iffCheYQBAIAmSQMAAKM2cXK3wUkaAACApnWTNDzlrPO61L3w7LO71E2SCx53Zrfa5z/u\njG61zzlna7fa28/c0q32th2bu9XevK3fr/qmLfPdak/mNuY7T0uLS91qHzxwqFvthX0Hu9Xe9thC\nv9odf7+3bO73+z23QX+/Dx7q9zvGySul/OMk/yHJ8ifquSSbk1ya5GtJ9s0un8yu95pa622z2781\nyZuTXJDke0murLXeNRvbmuQjSV6WZGuSrye5otb6wGz8oiQfS/LCJI8k+Wyt9d0rnfu6aRoAAOBE\nzI1keVKt9U+S/K13bksp70ny92df3lNr/YWj3baUclmSa5K8JMmfJ3lbki+UUp5ea92b5Nokz0ny\ngiR7knwyyc1JfmN2F7cnuTPJbyV5UpIvllLuq7V+eCVz1zQAAEAHs3f/357k2Un+7nGu/qYkN9da\nd85ue12mjcNlpZTbkrw+yatrrffOxq9O8v1SygVJnprkWUkurbU+muTRUsr1s9uvqGmwpwEAgFGb\ndPzvJL0vySdrrT+ZfX12KeX2UsquUsqPSylXLbvuJUnuOvxFrXUpyXeSPC/J05Ock+Tby8Zrkr2z\n212caYqxe9n93ZWklFJ2rGSimgYAAFhjpZSnJXlFkt+fXbQ7030K1yd5cqbJwTWllNfOxs9L8uAR\nd/NAkvNnY0tHGX9w2fjRbpvZ+HFZngQAAGvvt5PcXmu9P0lqrd/OdDP0YV8updyY5HVJbplddrxo\nozV+UrGIpAEAANbe5Uk+f5zr3JPkwtnnuzJNDJY7L8n9s7HJUcbPXTZ+tNsuzcaOS9MAAABrqJTy\ny0kuSvLlZZddXkq54oirPiPJD2ef78x0f8Lh689lulfhm7PrPHjE+DOTbJndbmeSi0op5y677+cn\n+X6tdc9K5mx5EgAAozbCM0I/J8nPZkcyOuxAkg+WUn6Q6TkWXpTktUleMxu/IcmtpZRbM9378M5M\nz+nwxVrrYinl40muLqXszHQD9LVJbqu17kqyq5RyZ5IPlFLekeQpSa5Kct1KJyxpAACAtXVBkvuW\nX1Br/XySK5N8NMnDSW5M8tZa6x2z8S8leU+SzyX5WZIXJ3lprXX/7C7em2nq8N0kd8/u443LSlye\nabNwX5KvJrml1nrjSic8WVrqd+bQ5d74K/+qy0Quevzje5RN4ozQPTgj9NpzRui154zQa29fxzNC\n73n0QLfaDz+8//hXGshPH9rbrfZ9Dz16/CsN5EcPHnkAnLX1if/7o+vyifW/fM5ruj3x/W/f/vS6\nfExONUkDAADQpGkAAACabIQGAGDUTsGZmTkOSQMAANAkaQAAYNTmJt4HH5pHGAAAaNI0AAAATZoG\nAACgSdMAAAA02QgNAMCoTSYOuTo0SQMAANCkaQAAAJosTwIAYNTmLE8anKQBAABokjQAADBqk0ga\nhiZpAAAAmiQNAACMmj0Nw5M0AAAATZoGAACgSdMAAAA0rZs9DRecfXaXuk88e0eXukly/uPO6Fb7\n8Y/f1q32mY/b2q32GWf1q731zC3dam/a2u9Xfb5j7cncxnxfZGlxsVvtQ/sPdqt9sGPtzdsOdKu9\naUu/n/O5uY25jnxxcalb7QOHDnWrzca2bpoGAAA4ERMboQe3Md+GAwAAVkzTAAAANFmeBADAqDlP\nw/AkDQAAQJOkAQCAUZtE0jA0SQMAANAkaQAAYNTsaRiepAEAAGjSNAAAAE2aBgAAoEnTAAAANNkI\nDQDAqE1shB6cpAEAAGjSNAAAAE2WJwEAMGrO0zA8SQMAANAkaQAAYNQmkTQMTdIAAAA0aRoAAIAm\ny5MAABg1G6GHJ2kAAACaNA0AAECTpgEAAGiypwEAgFGb2NMwOEkDAADQpGkAAACaLE8CAGDUHHJ1\neJIGAACgSdIAAMCo2Qg9vHXTNJy7Y3uXuo87a2uXuklyzjn9ap/5uH61t5+zrVvtbR1rb9nR7zGf\nP2Nzv9pb+tWezG/MMHXp0GK32ocOLHSrvWlvv9pzm+e71Z7MbcwXS4uLS91qH1g41K32ngN9Xi/B\nxvyLCgAArNi6SRoAAOBETLIxE7e1JGkAAACaNA0AAECTpgEAAGjSNAAAAE02QgMAMGob9MjDa0rS\nAAAANEkaAAAYNWeEHp6kAQAAaJI0AAAwanOShsFJGgAAgCZNAwAA0GR5EgAAo2Yj9PAkDQAAQJOm\nAQAAaNI0AAAATZoGAACgyUZoAABGbS42Qg9N0gAAADRJGgAAGDWHXB2epAEAAGjSNAAAAE2WJwEA\nMGpzlicNTtIAAAA0SRoAABg1QcPwJA0AAECTpgEAAGjSNAAAAE2aBgAAoGndbIQ+Z/uWLnXPOrNP\n3STZ3rH2GWdt7VZ72znbutXees4Z3Wpv3tHv3z1/Rr/v99zmfk8zk/n5brV7Wjp0qFvt+YWD3Wof\n2rK/W+3Jpo35HtzS4lK32gcPLHarfda+fj/n5+xd6FZ7PXPI1eFtzGc5AABgxdZN0gAAACdiEknD\n0CQNAABAk6YBAABosjwJAIBRm9gIPThJAwAA0CRpAABg1BxydXiSBgAAoEnSAADAqAkahidpAAAA\nmjQNAABA00kvTyqlbE7yxCT31lqXTn5KAADAerKqpKGU8pFln+8opdyS5LEkP0ryWCnl90spW07t\nFAEAgJ5WmzS8McnbZp9fn+Q5Sf5ZknuSPCPJ+5IsJPnvT9H8AACgySFXh7fapmH5d+SfJ/kHtdY6\n+/ovSyn/b5JvRNMAAACnjdVuhF6+Z+GRJD88YvyeJNtOZkIAAMD6suqkoZTyc5kmDv9PklcnuXnZ\n+JVJ/vwUzQ0AAI5rEsuThrbapmFrpmnC4e/MRZk1DaWUDyZ5U5KXnqrJAQDA6aaUcnWS305yVpI/\nTfLGWutfl1IuTfL+JL+U6YGG3l9r/cyy2701yZuTXJDke0murLXeNRvbmuQjSV6W6Wv2rye5otb6\nwGz8oiQfS/LCTFcMfbbW+u6VznlVy5NqrXO11vnZ/+dqrb+ybPh/TvJLtdZvrOY+AQDgZMxNJt0+\nVquU8ttJXpnkv0jy5CTfT3JVKeWCJHdk+sL+CZmu4PlEKeXi2e0uS3JNpit9npTkC0m+UEo5Y3bX\n12Z6kKIXJPnFTF/nL18RdHuSHyd5WpJfS/KKUsqVK533SZ+n4bBa63dO1X0BAMBp6u1J3l5r/cHs\n6yuTpJTyjiS11vqp2eVfKaV8PskbMk0X3pTk5lrrztn1r8v0qKaXlVJuS/L6JK+utd47G786yfdn\nzchTkzwryaW11keTPFpKuX52+w+vZNLOCA0AAGuglHJhkr+T5LxSyl+UUn5aSvlcKeX8JJckueuI\nm9yV5Hmzz//W+Oykyt+ZjT89yTlJvr1svCbZO7vdxUnuqbXuPuK+Syllx0rmfsqSBgAA6GFEp2l4\n6uz/lye5NMl8ktuSfCLJ9kyXDy33QJLzZ5+fl+TBY4yfl+lRTo8cf3DZ+NFum9n4Y8ebuKQBAADW\nxuH25vdqrf95tpTomkxPlryUHPcwUCczflKtlaYBAADWxn2z/z+87LJ7Mn1BvznTRGC585LcP/t8\nV2N81+w+jhw/d9n40W67NBs7Lk0DAACsjb9JsjvJs5dd9neSHEjyxSTPPeL6z0vyrdnnOzPdn5Ak\nKaXMZbpX4ZuZnnD5wSPGn5lky+x2O5NcVEo5d9l9Pz/J92ute1YycXsaAAAYtclINjXUWg+VUv59\nkqtLKX+S6fkS/k2STyf5n5L8m1LK65P8YZIXJ/n1TA+hmiQ3JLm1lHJrpudoeGeSfUm+WGtdLKV8\nfHa/OzPdAH1tkttqrbuS7Cql3JnkA7OjND0lyVVJrlvp3CUNAACwdt6T5P9I8mdJ/lOSmuRtsxf3\nL0/yliQPJflQklfVWv8iSWqtX5rd9nNJfpZpU/HSWuv+2f2+N9PU4btJ7s50CdQbl9W9PNNm4b4k\nX01yS631xpVOWtIAAABrpNZ6INPG4C1HGftGpidoO9Ztb0py0zHGFo51v7PxezM9W/QJ0TQAADBq\nJ3JmZlbH8iQAAKBp3SQNO87Y3KXu9u196ibJth39am89c0u32lt2bO1We/OObd1qb+pYe377Gd1q\nz23u93M+mduY74ssLS52q724sNCt9mRuY77TuHSw5/f7ULfaC/sOdqu9fU+/57Ver5fWO0HD8Dbm\nX1QAAGDFNA0AAEDTulmeBAAAJ8JG6OFJGgAAgCZNAwAA0KRpAAAAmjQNAABAk43QAACM2iQ2Qg9N\n0gAAADRJGgAAGLWJQ64OTtIAAAA0SRoAABi1OUHD4CQNAABAk6YBAABosjwJAIBRsxF6eJIGAACg\nSdMAAAA0aRoAAIAmTQMAANBkIzQAAKNmI/TwJA0AAECTpAEAgFFzRujhSRoAAIAmSQMAAKNmT8Pw\nJA0AAECTpgEAAGiyPAkAgFGzOml4kgYAAKBJ0wAAADRpGgAAgKZ1s6dh29Y+U9myrd9DsLlj7U2d\nHu8kmT9jc8faW/vV3n5Gv9pb+/2757b0qz2Z35jviywdWuxWezK3QR/zxaVutecPLHSrvWnvxvwb\n2vO1Q6/XS+AnDwCAUZuzE3pwG/MtIQAAYMUkDQAAjNokkoahSRoAAIAmTQMAANBkeRIAAKNmH/Tw\nJA0AAECTpAEAgFFzyNXhSRoAAIAmTQMAANCkaQAAAJo0DQAAQJON0AAAjNrERujBSRoAAIAmSQMA\nAKMmaBiepAEAAGjSNAAAAE2WJwEAMGo2Qg9P0gAAADRJGgAAGLU5QcPgJA0AAECTpgEAAGjSNAAA\nAE2aBgAAoMlGaAAARs0hV4cnaQAAAJokDQAAjJqgYXiSBgAAoEnTAAAANFmeBADAqM1ZnzS4ddM0\nbN7cJ/TYtKlf2LJpy3y32vNb+33r57ds7lZ7bnO/f/fc5o7/7i1bN2TtyfzGDFOXDi32nkIXS4v9\n/t09n1t6Pqf2/FvS829oz9cOvV4vwbppGgAA4EQ45OrwtKsAAECTpgEAAGjSNAAAAE2aBgAAoMlG\naAAARs0+6OFJGgAAgCZJAwAAo+aQq8OTNAAAAE2aBgAAoMnyJAAARs3qpOFJGgAAgCZJAwAAozYn\nahicpAEAAGjSNAAAAE2aBgAAoEnTAAAANNkIDQDAqNkHPTxJAwAA0KRpAAAAmixPAgBg1CbWJw1O\n0gAAADRJGgAAGDVBw/AkDQAAQJOkAQCAUbOnYXiSBgAAoEnTAAAANGkaAACAJk0DAADQZCM0AACj\nZh/08NZN07Bpvk/oMbep30/ZZK5n7X4h06TT93pae75f7Q37mHesPdfv+71R9f1+b9Dnlo36mHf8\nG9rztUOv10vgJw8AAGhaN0kDAACciDnrkwYnaQAAAJokDQAAjNpYg4ZSyu8neVutda6U8qtJvpZk\n32x4kmQpyWtqrbfNrv/WJG9OckGS7yW5stZ612xsa5KPJHlZkq1Jvp7kilrrA7Pxi5J8LMkLkzyS\n5LO11nevdK6aBgAAWGOllGcneU2mjcFh99Raf+EY178syTVJXpLkz5O8LckXSilPr7XuTXJtkuck\neUGSPUk+meTmJL8xu4vbk9yZ5LeSPCnJF0sp99VaP7yS+VqeBAAAa6iUMklyQ5IPreJmb0pyc611\nZ611f5LrMm04LiulzCd5fZL31VrvrbU+lOTqJC8vpVxQSnlukmcleVet9dFa691Jrp/d54poGgAA\nGLXJZNLt4wRdkWRvks8ccfnZpZTbSym7Sik/LqVctWzskiR3Hf6i1rqU5DtJnpfk6UnOSfLtZeN1\nVuOSJBdnmmLsXnZ/dyUppZQdK5mwpgEAANZIKeVJSX4nyX97xNDuTPcpXJ/kyZkmB9eUUl47Gz8v\nyYNH3OaBJOfPxpaOMv7gsvGj3Taz8eOypwEAANbOh5L8+1prLaX8/OELa63fTnLpsut9uZRyY5LX\nJblldtnxoo3W+EltF5c0AADAGiilvDjJP0zyu7OLjvdC/p4kF84+35VpYrDceUnun41NjjJ+7rLx\no912aTZ2XJoGAABGbTLp97FKr0ryxCQ/KqXsSvIfk0xKKfeXUl5dSrniiOs/I8kPZ5/vzHR/QpKk\nlDKX6V6Fb86u8+AR489MsmV2u51JLiqlnLvsvp+f5Pu11j0rmbjlSQAAsDauSvKvl339c0n+NMkv\nZ7qh+cZSyg8yPcfCi5K8NtPDsibToy3dWkq5NdO9D+/M9JwOX6y1LpZSPp7k6lLKzkw3QF+b5LZa\n664ku0opdyb5QCnlHUmeMpvLdSuduKYBAADWQK314SQPH/66lLI5yVKt9f9L8vlSypVJPpppM3Ff\nkrfWWu+Y3fZLpZT3JPlckidkes6Fl84Ov5ok701yZpLvJplP8keZngjusMuTfGJ2vw8nuaHWeuNK\n565pAABg1E7i0Kdd1Vr/OtMX+Ie//mSmJ2U71vVvSnLTMcYWkrxl9nG08XszPVv0CbGnAQAAaJI0\nAAAwaiMNGkZF0gAAADRpGgAAgCbLkwAAGLWxboQeE0kDAADQpGkAAACaNA0AAECTpgEAAGiyERoA\ngFGzD3p4kgYAAKBJ0gAAwKg55OrwJA0AAECTpAEAgFETNAxP0gAAADRpGgAAgCbLkwAAGLU565MG\nJ2kAAACaJA0AAIyaoGF4kgYAAKBJ0wAAADRpGgAAgCZNAwAA0GQjNAAAozaxE3pwkgYAAKBJ0gAA\nwKgJGoa36qahlHJxkjckeV6S82cX35/kW0luqrX+xambHgAA0NuqlieVUl6Z5E+SPD7J55L8j7OP\n25I8Ocm3Sim/caonCQAA9LPapOE9Sf5ZrfUrRxsspbw0yYeS3HGyEwMAgJWYzFmfNLTVboS+KNOk\n4Vi+kuRpJzwbAABg3Vlt0/CXSf5FY/xVs+sAAMCamEz6fWwUJ7I86Y5SyluS3JXkwSSTJOcluSTT\nlOHlp3KCAABAX6tKGmqtX03y95LcmuTMJM9O8stJtia5Jcnfq7V+4xTPEQAA6GjVh1yttd6X6WZn\nAABgAzjlZ4Qupew51fcJAAD0M8QZoTfQlhAAAHqbbKQdyZ2sqmkopXzmVN8nAACwvq32Bf6lSf4q\nyd0DzAUAAFiHVts0/MskH07yslrrI0e7QinlN096VgAAsEJWJw1vtYdc/VKSP8i0eTgW3zYAADiN\nnMghV687zvgZJz4dAABYHRuhh3fKD7kKAACcXhzpCACAURM0DE/SAAAANGkaAACAJk0DAADQpGkA\nAACa1s1G6IOHFrvUXTy41KVukiwt9qzd5/FOkqVO3+tp7UP9am/Yx7xf7Y1qo36/+/6OdXxu2aiP\nece/oT1fO/R6vbTu2Qk9OEkDAADQpGkAAACa1s3yJAAAOBHOCD08SQMAANAkaQAAYNQEDcOTNAAA\nAE2aBgAAoMnyJAAARm0yZ33S0CQNAABAk6YBAABo0jQAAABN9jQAADBqDrk6PEkDAADQpGkAAACa\nLE8CAGDUJtYnDU7SAAAANEkaAAAYNUHD8CQNAABAk6YBAABosjwJAIBRsxF6eJIGAACgSdMAAAA0\naRoAAIAmexoAABg1WxqGJ2kAAACaNA0AAECT5UkAAIyaQ64OT9IAAAA0SRoAABg3b4MPzkMMAAA0\nrZukYWFhsUvdgwf71E2SgwcOdat9aP/BfrUPLHSrPb/Q79+9uNDv3z2Z25jvD0zmN+a/e+lQv+e1\nxQP7+9Xu+Du22PG5pedzas+/JT3/hvZ87dDr9RKsm6YBAABOhI3Qw9uYb8MBAAArpmkAAACaNA0A\nAECTpgEAAGiyERoAgFGzD3p4kgYAAKBJ0gAAwKg55OrwJA0AAECTpAEAgFETNAxP0gAAADRpGgAA\ngCbLkwAAGDfrkwYnaQAAAJo0DQAAQJOmAQAAaNI0AAAATTZCAwAwapM5G6GHJmkAAACaJA0AAIya\nI64OT9IAAAA0SRoAAGCNlFJ+OcmHkjw3yd4kf5zkrbXW+0splyZ5f5JfSvKjJO+vtX5m2W3fmuTN\nSS5I8r0kV9Za75qNbU3ykSQvS7I1ydeTXFFrfWA2flGSjyV5YZJHkny21vrulc5b0gAAwKhNJpNu\nH6tRStmS5EtJvprkCUmemeRJSW4opVyQ5I5MX9g/IcmVST5RSrl4dtvLklyT5NWz23whyRdKKWfM\n7v7aJM9J8oIkv5jp6/ybl5W/PcmPkzwtya8leUUp5cqVzl3TAAAAa2N7kv8hyQdqrQu11p9l+mL+\nmUlelaTWWj9Vaz1Qa/1Kks8necPstm9KcnOtdWetdX+S65IsJbmslDKf5PVJ3ldrvbfW+lCSq5O8\nvJRyQSnluUmeleRdtdZHa613J7l+dp8romkAAGDUJpN+H6tRa32o1voHtdbFJCmllCSvTfLZJJck\nueuIm9yV5Hmzz//WeK11Kcl3ZuNPT3JOkm8vG6+ZLn+6JMnFSe6pte4+4r5LKWXHSuZuTwMAAKyh\n2f6C/5RkPsnHk/xOkv890+VDyz2Q5PzZ5+clefAY4+dlmjocOf7gsvGj3Taz8ceON2dJAwAArKFa\n649qrVuTlNnHp2dDx8suTmb8pA5Mq2kAAIAOZnsLrk7yL5IcyDQRWO68JPfPPt/VGN+VaVNw5Pi5\ny8aPdtul2dhxaRoAAGANlFJeVEr5yyMuXpp9/Fmmh2Fd7nlJvjX7fGem+xMO39dcpnsVvpnkh5ku\nP1o+/swkW2a325nkolLKucvu+/lJvl9r3bOSudvTAADAuI3nlND/McnZpZTfy3Qfw5mZHkb1/0py\nQ5J3lFJen+QPk7w4ya9negjVzMZvLaXcmuk5Gt6ZZF+SL9ZaF0spH09ydSllZ6YboK9NclutdVeS\nXaWUO5N8oJTyjiRPSXJVpkdgWpF10zTs23+wS90D+/rUTZKFjrUPdnq8k2TT3oVutQ9t2d+t9mRu\nNE9op9TS4mK32pO5jRmm9nzMFxc6/n7v2duv9t5+zy2HOj6n9vxb0vNvaM/XDr1eL3Fq1Fp3l1L+\naZKPZros6NFMz9nw39Raf1pKeXmSf5vk3yW5J8mraq1/Mbvtl0op70nyuUzP43BnkpfODr+aJO/N\ntAn5bqYbrP8o0xPBHXZ5kk8kuS/Jw0luqLXeuNK5r5umAQAATsSY3pibNQEvOsbYNzI9QduxbntT\nkpuOMbaQ5C2zj6ON35vp2aJPyMZ8Gw4AAFgxTQMAANBkeRIAAKM2nn3Q4yVpAAAAmiQNAACMm6hh\ncJIGAACgSdMAAAA0aRoAAIAmTQMAANBkIzQAAKNmH/TwJA0AAECTpAEAgFGbzIkahiZpAAAAmjQN\nAABAk+VJAACM2sRO6MFJGgAAgCZJAwAA4yZoGJykAQAAaNI0AAAATZoGAACgSdMAAAA02QgNAMCo\nOeTq8CQNAABAk6QBAIBRkzQMT9IAAAA0aRoAAIAmy5MAABg3b4MPzkMMAAA0SRoAABg1G6GHt26a\nhsf2LnSpu2dPn7pJsu2xfrU3bzvQrfbc5vlutSebNma4trS41K323OZ+TzOT+X4/az0tHTrUrfbi\nwsFutQ/UIRE+AAAOdElEQVTt3d+t9sJj+7rVPvBYv3/3/kf7/S3Z1/FvaM/XDr1eL8HGfAUFAACs\nmKYBAABo0jQAAABN62ZPAwAAnAgboYcnaQAAAJo0DQAAQJPlSQAAjJvVSYOTNAAAAE2SBgAARm0y\nJ2oYmqQBAABokjQAADBuDrk6OEkDAADQpGkAAACaNA0AAECTpgEAAGiyERoAgFGzD3p4kgYAAKBJ\n0wAAADRZngQAwKhNrE8anKQBAABokjQAADBuc5KGoUkaAACAJk0DAADQZHkSAACjZiP08CQNAABA\nk6YBAABo0jQAAABN9jQAADButjQMTtIAAAA0aRoAAICmdbM86eE9B7rU3fHo5i51k2Tbtn4P/6Yt\n/frFyQY9a+PSwcVutecPLPSrvaXf79hkfmO+L7J0qN/P2qGOP2uH9varfeCx/d1q73t4X7faex/p\n9+/e82if1w1J8kjH2r1eL613Drk6vI35FxUAAFixdZM0AADAidioqxjWkqQBAABo0jQAAABNlicB\nADBuNkIPTtIAAAA0SRoAABg1h1wdnqQBAABo0jQAAABNq1qeVEr5r2qtty37+rVJXpfkwiQ/TPKx\nWusdp3SGAABAV6tNGj59+JNSylVJPpzkT5N8KElN8ulZIwEAAJwmVrsRevkuk7cn+ee11i8fvqCU\n8kdJbkxyy8lPDQAAVsA+6MGtNmlYWvb5tiRfO2L8a0kuOKkZAQAA68rJHHL1T5K8MMk3ll32T5L8\n+GQmBAAAqzGZEzUMbbVNw7ZSyg9nn5+d5HFJLk2SUsoVST6Y5MpTNz0AAKC31TYNLzri693LPv9p\nklfWWj9/clMCAADWk1U1DbXWP26M/a8nPx0AAFglZ4Qe3Ck/uVspZc+pvk8AAKCfk9kIfSxaPQAA\n1sxE0jC41Z4R+jOn+j4BAID1bbUv8C9N8ldJ7h5gLgAAwDq02qbhXyb5cJKX1VofOdoVSim/edKz\nAgAA1o1VbYSutX4pyR9k2jwci0VlAABwGln1/oNa63XHGT/jxKcDAACr5IzQgzvlh1wFAABOL450\nBADAqDnk6vAkDQAAQJOmAQAAaLI8CQCAcbM6aXCSBgAAoGndJA0PPLanS93tW/o9BFs2z3erPbdB\nD022tLjUrfbiwqFutTft7fdzPr+1X+3J3MZ8X2RpcbFb7UP7D3arfbBj7f2PHuhWe+8j+7vVfvSh\nfrUffrhf7Yc6Pua9Xi+tdzZCD29j/kUFAABWTNMAAAA0aRoAAIAmTQMAANC0bjZCAwDACdmgB3hZ\nS5IGAACgSdMAAAA0WZ4EAMCoOU/D8DQNAACwRkopL0nyqSRfrbW+ctnlv5rka0n2zS6aJFlK8ppa\n622z67w1yZuTXJDke0murLXeNRvbmuQjSV6WZGuSrye5otb6wGz8oiQfS/LCJI8k+Wyt9d0rnbem\nAQCAcRtJ0lBKeWeS1yf5q2Nc5Z5a6y8c47aXJbkmyUuS/HmStyX5Qinl6bXWvUmuTfKcJC9IsifJ\nJ5PcnOQ3Zndxe5I7k/xWkicl+WIp5b5a64dXMnd7GgAAYG3sTfL8JHefwG3flOTmWuvOWuv+JNdl\nmkRcVkqZz7QZeV+t9d5a60NJrk7y8lLKBaWU5yZ5VpJ31VofrbXeneT62X2uiKYBAIBRm0wm3T5W\no9b60VrrI42rnF1Kub2UsquU8uNSylXLxi5Jctey+1pK8p0kz0vy9CTnJPn2svGaaZNySZKLM00x\ndi+7v7uSlFLKjpXMXdMAAAD97c50n8L1SZ6caXJwTSnltbPx85I8eMRtHkhy/mxs6SjjDy4bP9pt\nMxs/LnsaAACgs1rrt5NcuuyiL5dSbkzyuiS3zC47XrTRGj+pjR+SBgAAWJ/uSXLh7PNdmSYGy52X\n5P7Z2OQo4+cuGz/abZdmY8elaQAAgM5KKZeXUq444uJnJPnh7POdme5POHz9uUz3Knxzdp0Hjxh/\nZpIts9vtTHJRKeXcZff9/CTfr7XuWcn8LE8CAGDc5sZxyNXjOJDkg6WUH2R6joUXJXltktfMxm9I\ncmsp5dZM9z68M9NzOnyx1rpYSvl4kqtLKTsz3QB9bZLbaq27kuwqpdyZ5AOllHckeUqSqzI9AtOK\naBoAAGANlFL2ZrokaPPs61ckWaq1bq+1fr6UcmWSjyb5uST3JXlrrfWOJKm1fqmU8p4kn0vyhEzP\nufDS2eFXk+S9Sc5M8t0k80n+KNMTwR12eZJPzO734SQ31FpvXOncNQ0AALAGaq1nHGf8k5melO1Y\n4zcluekYYwtJ3jL7ONr4vZmeLfqEaBoAABi11Z4vgdWzERoAAGiSNAAAMG6ShsFJGgAAgCZJAwAA\nozY5PQ65uq5JGgAAgCZNAwAA0KRpAAAAmjQNAABAk43QAACMm0OuDk7SAAAANK2bpOG+3bu71N0y\nP9+lbpLMbdDDgy0uLnWrffDAYrfaC/sOdqu9eVu/X/VNW/r9jm3UQ/Atdf0dO9Stds/fsX2PLXSr\nvefRA91qP/zw/m61f/rQ3m6179/9WLfavV4vwbppGgAA4ERMLE8anOVJAABAk6QBAIBxkzQMTtIA\nAAA0aRoAAIAmy5MAABi1jXq0vLUkaQAAAJo0DQAAQJOmAQAAaLKnAQCAcXPI1cFJGgAAgCZNAwAA\n0GR5EgAA42Z50uAkDQAAQJOkAQCAUZtIGgYnaQAAAJo0DQAAQJPlSQAAjNuc5UlDkzQAAABNmgYA\nAKBJ0wAAADRpGgAAgCYboQEAGLXJxPvgQ/MIAwAATZIGAADGzRmhBydpAAAAmiQNAACM2kTSMDhJ\nAwAA0KRpAAAAmixPAgBg3OYsTxraumkafvLIz7rU3TQ/36Vub4uLS91qH1g41K32WfsOdqu9fc/m\nbrW3bOv3q75pU79Ac27Txvwjsniw3+/3wYOL3Wof6Pj7vWfPQrfajzx6oFvthx7Z3632/bsf61b7\n3t27u9Xu9XoJLE8CAACaNA0AAECTpgEAAGhaN3saAADgRDhPw/AkDQAAQJOkAQCAcZM0DE7SAAAA\nNEkaAAAYt4n3wYfmEQYAAJo0DQAAQJPlSQAAjNpkzkbooUkaAACAJk0DAADQpGkAAACaNA0AAECT\njdAAAIybM0IPTtIAAAA0SRoAABi1iaRhcJIGAACgSdMAAAA0WZ4EAMC4TbwPPjSPMAAA0CRpAABg\n1CZzNkIPTdIAAAA0aRoAAIAmTQMAANCkaQAAAJpshAYAYNycEXpwkgYAAKBJ0gAAwKhNJA2DWzdN\nw988fG/vKay5g4cOdat9oGPtPQe2d6t9zt6FbrV3nLG5W+1tW/v9qm/e3C/Q3DS/McPUg4cWu9Ve\nWOhXe9/+g91qP9bxueXhPQe61X7gsT3dat+3e3e32j955Gfdam/E10usDxvzLyoAALBi6yZpAACA\nEzLxPvjQPMIAAECTpAEAgHGbsxF6aJIGAACgSdMAAAA0aRoAAIAmTQMAANBkIzQAAKPmjNDDkzQA\nAABNkgYAAMbNyd0G5xEGAACaNA0AAECT5UkAAIyajdDDkzQAAABNkgYAAMbNRujBeYQBAIAmTQMA\nANCkaQAAAJo0DQAAQJON0AAAjNpkziFXhyZpAAAAmjQNAABAk+VJAACMmzNCD07SAAAANEkaAAAY\ntYkzQg/OIwwAADRJGgAAGDd7GgYnaQAAAJomS0tLvecAAACsY5IGAACgSdMAAAA0aRoAAIAmTQMA\nANCkaQAAAJo0DQAAQJOmAQAAaNI0AAAATZoGAACgSdMAAAA0aRoAAICmTb0ncDJKKRcl+ViSFyZ5\nJMlna63v7jsrTkellMUk+5MsJZnM/v+JWuvbuk6M0SulvCTJp5J8tdb6yiPGLk3y/iS/lORHSd5f\na/3M2s+S08GxftZKKb+a5GtJ9s0uOvwc95pa621rPlFgXRp105Dk9iR3JvmtJE9K8sVSyn211g/3\nnRanoaUkv1hr/XHviXD6KKW8M8nrk/zVUcYuSHJHkn+V5NYk/zjJ50spf1lrvWtNJ8rotX7WZu6p\ntf7CGk4JGJnRLk8qpTw3ybOSvKvW+mit9e4k1yd5U9+ZcZqazD7gVNqb5PlJ7j7K2KuS1Frrp2qt\nB2qtX0ny+SRvWMsJctpo/awBHNeYk4aLM31nZPeyy+5KUkopO2qtj3WaF6ev3yul/MMkZyX5X5K8\n3c8ZJ6PW+tEkKaUcbfiSTJ/TlrsryX898LQ4DR3nZy1Jzi6l3J5porUvyfW11t9fo+kBIzDapCHJ\neUkePOKyB2b/P3+N58Lp70+T/IckfzfJP8h0H82/6zojTnfHeo7z/MaptjvJ9zJN65+c6TKma0op\nr+05KWB9GXPSkFguwhqptf7K8i9LKe/KdH35G2utC73mxWnPcxyDq7V+O8mlyy76cinlxiSvS3JL\nl0kB686Yk4Zdmb4Tt9x5mW5Y3bX202GDuSfJfJIndp4Hp69jPcfd32EubDz3JLmw9ySA9WPMTcPO\nJBeVUs5ddtnzk3y/1rqn05w4DZVSnl1K+eARFz8j00Ow3tthSmwMOzPd17Dc85J8q8NcOI2VUi4v\npVxxxMXPSPLDHvMB1qfRLk+qtX6nlHJnkg+UUt6R5ClJrkpyXd+ZcRq6P8mbSin3J/lwkqcleV+S\nm2qtSz0nxmntD5P8Tinl9bPPX5zk15O8oOusOB0dSPLBUsoPknw9yYuSvDbJazrOCVhnJktL433N\nU0q5MMknkvyTJA8nuaHW+rtdJ8VpqZTyj5L8XpK/n+mRRW5J8q9rrQd6zotxK6XszXRJ5ebZRQeT\nLNVat8/G/1GSf5vpyd3uSfLuWusdHabKyK3gZ+0NSf67JD+X5L4kv1trvaXDVIF1atRNAwAAMLwx\n72kAAADWgKYBAABo0jQAAABNmgYAAKBJ0wAAADRpGgAAgCZNAwAA0KRpAAAAmjQNAABAk6YBAABo\n0jQAAABN/z8LFbTxE5IgnQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f88c7e42210>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(10, 8))\n",
"sns.heatmap(likelyhood_array, xticklabels=5, yticklabels=5)"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"cov_params = np.eye(len(w))\n",
"\n",
"width0 = 10\n",
"width1 = 10\n",
"num = 101\n",
"w0_array = np.linspace(w[0] - width0, w[0] + width0, num)\n",
"w1_array = np.linspace(w[1] - width1, w[1] + width1, num)\n",
"evidence_array = np.zeros([num, num])\n",
"for (i, w0_) in enumerate(w0_array):\n",
" for (j, w1_) in enumerate(w1_array):\n",
" pred_ = np.asarray(w0_ + w1_ * X[:,:1])[:, 0]\n",
" evidence_array[i, j] = data_likelyhood(y, pred_, cov_data) * \\\n",
" data_likelyhood(np.array([w0_, w1_]).reshape(-1, 1), w.reshape(-1, 1), cov_params)"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"Правдоподобие модели без учёта нормировки"
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3.0541290086635438"
]
},
"execution_count": 87,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"evidence_array.sum()"
]
},
{
"cell_type": "code",
"execution_count": 119,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"likelyhood_array = []\n",
"width = 10\n",
"num = 11\n",
"\n",
"for n_features in range(1,X.shape[1]):\n",
" lr = LinearRegression()\n",
" lr.fit(X[:,:n_features], y)\n",
" predicted = lr.predict(X[:,:n_features])\n",
" likelyhood_array.append(data_likelyhood(y, predicted, cov_data))"
]
},
{
"cell_type": "code",
"execution_count": 122,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.text.Text at 0x7f88b2a8d950>"
]
},
"execution_count": 122,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA08AAAF0CAYAAAD7HAJdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XuU5GV95/F3dfdMzzDDIMxwh4kX4OsFLzCAJihJRKOL\nIS4bs5uom2w4K7pmFYzxiDFHs8luoiEak00EghswRgnuyh4V3YCG1QTBhJEgN/misKMiDIzOyFy7\n+lb7x+9XMzVN93RVdfVUU7/365ymuur3VP2eoZ+pmk8/z/P91RqNBpIkSZKkAxvqdwckSZIk6anA\n8CRJkiRJbTA8SZIkSVIbDE+SJEmS1AbDkyRJkiS1wfAkSZIkSW0wPEmSJElSGwxPkiRJktQGw5Mk\nSZIktcHwJEnqqYi4JiIePcDxqyPikZb7myLiUws8509HxHRE/Fx5//0RMRURy8v7X4mImxdyjoWK\nwnRE/Go/+yFJ6p7hSZLUa43yay5vB54/o32vztt0GXBsZo73+BySpAob6XcHJEnVkpk7Fumlay3n\n2A3sXqTzSJIqyvAkSTqoIuIa4FWZeewcx48AbgO+D/yrzJyIiBcDvwucBqwGvgFcmpm3zfEavwu8\nD1jRMvtERLwS+DBwMvAIcElmfq7l+E8Cvw+cRfEZ+S3gjzLzupY2a4APAr8ArAMeA64HfrsMbUTE\nSHme1wOjwFeBP233/5EkaWly2Z4k6WCbc1lfRKwAbgCeAH6hDE4nA39P8Zn1KuDFwMPAlyLilA7O\n8XTgrcC/B86kCE9/ExGHlud+Tnme7cDLKILaPwLXRsT5La9zA/DzwFuAAH6zfM2/bmnzPuDNwHuB\nFwJ/A3xkrj+3JOmpwZknSdKSEBE14FrgcODs5iwORTiZBH4xM3eWbS8ENpXH3tLmKY4GLszMbeVr\n/Cnwt8DzgK8Dl1As9ft3mTlRPueSiHg58Dbg8+XM1EuBX8rMz5dtNkXEeuCyiDg+M38A/Drwmcy8\nsmzzUEQcC/xxR/9TJElLijNPkqSl4s8pZpVemZlbWx4/C/inZnACyMw68DXg9A5e/9vN4FTaQrFP\n6tDy/hnA7S3BqenWlvOcSTF7dMssbWrAaRFxGHA8cMccbSRJT1HOPEmSloLzgFXAGMWeplZrgOdH\nxMxCE8sp9hu1a2YBieYSumagWQN8e5bnbWdfwGrePjFLm+bxZpudM9osVqEMSdJB4syTJGkp2Aa8\nCEjguub1mVqO/SPwAor9Q82v51DsTeqVHwOHzfL4YeUxWm5ntjus5fiu8vtDZrR52kI7KEnqL8OT\nJGkpuC0z7wV+BXgm8Cctx74OPBt4ODMfan5RfIZtbmm30GIM/wScOSO4AZwN/HNLmxpwzow2LwOm\ngG+USwMfp1iC2OqcHvRRktRHLtuTJC2GoYg4epbHxw70pMz8dkS8DfiriPhyZv5vihLf/4Gi6t0f\nAj8CXkkRsC6l2CsFC99P9GfleT5VljqfBi6mqKj31rJ/GyPiZuBDEbGLopT52cBvA9dk5uPla30C\neFtZ2OKrFEHqDQvsnySpz5x5kiQthiMpSoHP/LqaJ8++7FdWPDOvAa4DroqIEzPzQeCnKfZE3Qzc\nT1H97h2Z+eczXudAZjveet4HgHMpqv3dBtwOnAq8JjP/oeU5/5qiXPnHyr78HkXAa6369zvANcBl\nwDeBNwIXztM/SdISV2s0XEEgSZIkSfPpeNleeS2LjwIvoagcdF1mXjpH27dTLHU4BriL4krud5TH\nVgAfAH6R4reJtwO/Wa557+g8kiRJkrTYulm2dz3wfYortb8CuCAiLpnZqLwa+/spliocTbHE4YaI\nWFk2+SOKdeIvobgexveA/93peSRJkiTpYOgoPEXEGRSlYt+dmTvLdegfBi6apflFwNWZubG8mOFl\nFGvLzy+P/xj4rcz8QWbuAT4CPCsijunwPJIkSZK06DqdeTod2JSZ21seuwOIiFg1o+0GWq6unpkN\n4E6Kq7OTme/LzK+2tF9PUYVpa4fnkSRJkqRF12l4WktxscJWW8vbdW22ndmOiDicolLRZZk53uF5\nJEmSJGnRdXOdp06uozFv24g4Fvg/wDeA/9LleZ6k0Wg0arWFXvJDkiRJ0gDoSTDoNDxtoZgVarWW\nYi/Tljbb3t28ExHPAr4MfB64uFza1+l5ZlWr1di+fQ9TU9PtNNcAGh4eYs2alY6DCnMMyDEgx4Ac\nA2qOgV7oNDxtBNZHxBGZ2VxGdxZwX2bunqXtBoqrrBMRQxR7mT5W3l8L3Ah8LDP/2wLOM6epqWkm\nJ/1LUnWOAzkG5BjQYo6BRqPBxtzCIz/ctSivr4UZGqqxcsUy9oxNMD29uNc3PeaIQzjrOUfh6qfB\n1VF4ysw7I+J24AMR8U6KEuPvoKikR0TcD1yYmbcClwPXRsS1FNd4ehdFQYgvlC/3AeDrswSnec8j\nSZK0VHzt7s381Re/1e9uaIk4/NBRTjnxaf3uhhZJN3ueXgdcBWwGngAuz8wrymMnA6sBMvPGiHgP\n8GngSIqL4J5Xli0H+HVgMiJ+kWI5Xq28fVNmfnKe80iSJPXdoz/axd98KQEYGa6xfGS4zz3Sk9SK\n7RyNRqP4l+YiOmbtIRy3zsLQg6zWaCzyKOqfxrZtu1ymUWEjI0McfvgqHAfV5RiQY0CLOQYmJqf4\nr3/9Db7/+E6Gh2q891c38PRj1vT0HFo43wdUjoGerKXstFS5JEmSgE/f/CDff3wnAL/0M88yOEkV\nYHiSJEnq0L88sIW/v+NhAF7wrLW88swT+9wjSQeD4UmSJKkDW7eP7S0Qcdjq5Vz4mudYXU2qCMOT\nJElSm6amp/nLz93LrrFJasBFP/9c1hyyvN/dknSQGJ4kSZLa9PmvbeKBh58A4DU/9RM85+lH9LlH\nkg4mw5MkSVIb8nvb+PytmwA46YTDeO1Ln9HfDkk66AxPkiRJ89i5Z4K//Px9NBpwyOgIF53/XIaH\n/GeUVDX+rZckSTqARqPBX33hW2zbUQfg1897NusOW9nnXknqB8OTJEnSAXz5Gw9z53d+CMDPnnY8\nG+KoPvdIUr8YniRJkubw3c07+J//9zsAnHDkKv7dy0/qc48k9ZPhSZIkaRZj45Nc8bl7mZxqsHxk\niDe/9lSWLxvud7ck9ZHhSZIkaRafvOkBHtu6G4DXv/IUjl+3qs89ktRvhidJkqQZbrtnM1+7ZzMA\nZz77KF72gmP73CNJS4HhSZIkqcVjW3fz1zclAOsOW8GvvfrZ1Gq1PvdK0lIw0u8OSJJUNXd++4d8\n4eub2Llnst9dGXg1YHi4xtRUg0abz9m5e5z6+BTDQzXe/NrnccgK/7kkqeC7gSRJB8m2HXU+9aUH\n+MYDW/rdFbXhgnOeybOOO6zf3ZC0hBieJElaZNPTDW6+42Gu/4eHGBufAuDwQ0d50UnriqkRLZqh\nWo3R0RHq9UmmG+3OPcExRxzCuRtOWMSeSXoqMjxJkrSIvvfYDj7+d/fz/x7dAUCtBuduOIELXvZM\nVo76MbzYRkaGOPzwVWzbtovJyel+d0fSU5zv2pIkLYKx8Uk+e8v/40u3P7x3xmP90av5tVc/m2cc\nu6bPvZMkdcPwJElSj935nR/yyZuSH22vAzC6bJgLXvYMzj3jBIaHLHQrSU9VhidJknpk24461375\nATbmvoIQLzppHW945SmsPWxFH3smSeoFw5MkSQs0Pd3g//7LD/jMVx/cWxDiaauX84ZXnsLppxzp\nNYIkaUAYniRJKk1OTfM3Nz3Aoz/a1dHztu+e4LGtu4GieN7LN5zAvznHghCSNGh8V5ckqXTPQ1v5\nh28+0vXzTzyqKAjxzOMsCCFJg8jwJElSacee8b3fn3byOoaH2lxuV6sRJz6NnzntOAtCSNIAMzxJ\nklRq7lcC+I0Lns9Qu+FJklQJ/npMkqRSvQxPy5cNGZwkSU9ieJIkqdSceVqxbLjPPZEkLUWGJ0mS\nSs2ZpxXLXdUuSXoyw5MkSaWx8UkARpc78yRJejLDkyRJpbGJ5syT4UmS9GSGJ0mSSs1le848SZJm\nY3iSJKk05p4nSdIBGJ4kSSpZbU+SdCCGJ0mSShaMkCQdiOFJkqRS3YIRkqQDMDxJklTat+fJ8CRJ\nejLDkyRJwNT0NBOT04AFIyRJszM8SZIE1Men934/asEISdIsDE+SJLGvWAS4bE+SNDvDkyRJ7CsW\nAYYnSdLsDE+SJLGvWARYqlySNDvDkyRJ7B+eLBghSZqN4UmSJPbf8+TMkyRpNoYnSZKA+rh7niRJ\nB2Z4kiQJGGstGGGpcknSLAxPkiQBY/UiPA3Vaiwb8eNRkvRkfjpIksS+UuWjy4ep1Wp97o0kaSky\nPEmSxL6CEe53kiTNxfAkSRL7CkYYniRJczE8SZLEvus8jVosQpI0B8OTJEnsC0/OPEmS5mJ4kiSJ\nfQUjViwf6XNPJElLleFJkiQsGCFJmp/hSZIkWvY8GZ4kSXMwPEmShAUjJEnzMzxJkoSlyiVJ8zM8\nSZIqr9FotFTbs2CEJGl2hidJUuVNTk0z3WgAzjxJkuZmeJIkVV5z1gksGCFJmpvhSZJUea3hyZkn\nSdJcDE+SpMqrt4Ynq+1JkuZgeJIkVd7+y/YsGCFJmp3hSZJUeWMTk3u/d9meJGkuhidJUuXV3fMk\nSWqD4UmSVHkWjJAktaPjhd0RsR74KPASYAdwXWZeOkfbtwNvBY4B7gIuycw7Wo6fBPwtcFxmHjfj\nudNAHWgAtfL2qsy8uNM+S5J0IJYqlyS1o5tdsdcDtwO/DBwNfDEiNmfmR1obRcT5wPuBVwF3AxcD\nN0TEszJzT0T8LPAJ4FZgv+BUagCnZOb3u+ijJEltGxsv9jwtGxlieMhFGZKk2XX0CRERZwAvAN6d\nmTsz80Hgw8BFszS/CLg6MzdmZh24jCIQnV8ePwI4F/jCHKerlV+SJC2q+kQx8zRqmXJJ0gF0+uu1\n04FNmbm95bE7gIiIVTPabiiPAZCZDeBO4Mzy/mcyM+c53wcj4rsRsTUirpzlHJIkLdhYvQhP7neS\nJB1Ip+FpLbBtxmNby9t1bbad2W4utwE3AScBP0mxx+ov2u6pJEltGpswPEmS5tfNnqdOltJ1vewu\nM89uvRsR7wY+FxFvysyJdl5jeNh161XW/Pk7DqrLMaB2x8D45DQAK0ZHGBlxvAwS3wfkGFAvf/ad\nhqctFDNKrdZS7GXa0mbbuzs8Z9MmYBg4CvhBO09Ys2Zll6fSIHEcyDGg+cbA1HQDgEMPWc7hh7tC\nfBD5PiDHgHqh0/C0EVgfEUdkZnO53lnAfZm5e5a2Gygq6hERQxR7pj4230ki4kXAGzPzt1oefi5F\n6fJH2u3s9u17mJqabre5Bszw8BBr1qx0HFSYY0DtjoGdu8aL9kOwbduug9U9HQS+D8gxoOYY6IWO\nwlNm3hkRtwMfiIh3AscD76CopEdE3A9cmJm3ApcD10bEtRTXeHoXMMaTq+vNtrTvceCiiHgc+Ajw\ndOD3gCvLwhNtmZqaZnLSvyRV5ziQY0DzjYE99aJU+fKRYcfKgPJ9QI4B9UI3e55eB1wFbAaeAC7P\nzCvKYycDqwEy88aIeA/waeBIimtDnVeWLScibgTOoShaMRIReyiW//1cZt4SEecBHwR+hyJ0XVN+\nL0lST1kwQpLUjo7DU2Y+ArxmjmPDM+5fCVw5R9tXzXOeW4CzD9RGkqReGBsvr/NkeJIkHYBlRyRJ\nlVcfb848dbMgQ5JUFYYnSVKlTTca1JvL9pY58yRJmpvhSZJUac1ZJ3DPkyTpwAxPkqRKa846gXue\nJEkHZniSJFXamDNPkqQ2GZ4kSZW2/7I9C0ZIkuZmeJIkVdrY+OTe70ctGCFJOgDDkySp0vZbtjdq\neJIkzc3wJEmqtNaCEZYqlyQdiOFJklRpY+55kiS1yfAkSaq0ZniqAcuW+bEoSZqbnxKSpEprFoxY\nvnyYoVqtz72RJC1lhidJUqU1S5V7jSdJ0nwMT5KkSmsu27NYhCRpPoYnSVKl7Q1PFouQJM3D8CRJ\nqrRmqfJRl+1JkuZheJIkVVqzYIR7niRJ8zE8SZIqzYIRkqR2GZ4kSZXW3PM0asEISdI8DE+SpEqz\nYIQkqV2GJ0lSpTX3PFkwQpI0H8OTJKnSmtX2VhqeJEnzMDxJkiprcmqayakG4MyTJGl+hidJUmU1\n9zuB1fYkSfMzPEmSKqveEp5Gl1kwQpJ0YIYnSVJlNYtFgDNPkqT5GZ4kSZU1NuGyPUlS+wxPkqTK\nat3zZMEISdJ8DE+SpMqqWzBCktQBw5MkqbL2D08WjJAkHZjhSZJUWRaMkCR1wvAkSaqsZsGI4aEa\nI8N+JEqSDsxPCklSZY3Vi/DkrJMkqR2GJ0lSZdUnDE+SpPYZniRJldXc8zRqsQhJUhsMT5Kkympe\n58mZJ0lSOwxPkqTKapYqH11meJIkzc/wJEmqLGeeJEmdMDxJkiprzIIRkqQOGJ4kSZXVnHmyYIQk\nqR2GJ0lSZdXLanvOPEmS2mF4kiRV1t7rPFkwQpLUBsOTJKmSGo2GBSMkSR0xPEmSKml8cppGo/h+\n1PAkSWqD4UmSVEnNWScwPEmS2mN4kiRVUrNYBMAKq+1JktpgeJIkVVLrzJMFIyRJ7TA8SZIqab/w\nNGp4kiTNz/AkSaqkZplygFFnniRJbTA8SZIqab+ZJ/c8SZLaYHiSJFXS2H4FI5x5kiTNz/AkSaqk\n/UqVu2xPktQGw5MkqZLqZXhavmyIoaFan3sjSXoqMDxJkiqpOfNkmXJJUrsMT5KkSmrOPFksQpLU\nLsOTJKmSxiaKghGjFouQJLXJ8CRJqqS9y/YMT5KkNhmeJEmV1Fy258yTJKldhidJUiVZMEKS1CnD\nkySpksYsGCFJ6pDhSZJUSWPjFoyQJHXG8CRJqqT6hAUjJEmdMTxJkiqpbrU9SVKHDE+SpMqZmp5m\nfHIacM+TJKl9hidJUuXUx6f3fj9qtT1JUpsMT5KkymkWiwCX7UmS2md4kiRVTrNYBBieJEnt63ih\nd0SsBz4KvATYAVyXmZfO0fbtwFuBY4C7gEsy846W4ycBfwscl5nHdXseSZI60bzGE1iqXJLUvm5m\nnq4Hvg88HXgFcEFEXDKzUUScD7wfeCNwNHADcENErCyP/yzwFeChhZxHkqROtYYnC0ZIktrVUXiK\niDOAFwDvzsydmfkg8GHgolmaXwRcnZkbM7MOXAY0gPPL40cA5wJfWOB5JEnqSN2ZJ0lSFzqdeTod\n2JSZ21seuwOIiFg1o+2G8hgAmdkA7gTOLO9/JjOzB+eRJKkjFoyQJHWj07UKa4FtMx7bWt6uA3a1\n0XZdj88zp+Fh62FUWfPn7zioLseA5hoDE1P7SpWvXrmMkRHHyKDyfUCOAfXyZ9/NQu/aIrXt5XMB\nWLNm5UJfQgPAcSDHgGaOgdpwMds0VIOjjjyUWm3BHzla4nwfkGNAvdBpeNpCMSvUai3FXqYtbba9\nu8fnmdP27XuYavntoqpleHiINWtWOg4qzDGgucbAtif2AEWxiB//eHe/uqeDwPcBOQbUHAO90Gl4\n2gisj4gjMrO5jO4s4L7MnPnps5Fi39MnACJiiGIv08d6fJ45TU1NMznpX5KqcxzIMaCZY2D32ARQ\nFItwbFSD7wNyDKgXOloAmJl3ArcDH4iIQyPi2cA7KK7HRETcHxE/VTa/HPjViHhxWZ78d4Axnlxd\n70lrJeY7jyRJC9GstmexCElSJ7rZ8/Q64CpgM/AEcHlmXlEeOxlYDZCZN0bEe4BPA0dShKHzyrLl\nRMSNwDkUAW4kIvZQLMv7ucy8ZZ7zSJLUtbGJIjyNLjM8SZLa13F4ysxHgNfMcWx4xv0rgSvnaPuq\nbs8jSdJCjNWdeZIkdc6ajZKkyqlPNMNTNwswJElVZXiSJFVO8yK5o848SZI6YHiSJFXOmAUjJEld\nMDxJkiqnGZ4sGCFJ6oThSZJUOZYqlyR1w/AkSaqURqNhwQhJUlcMT5KkSpmcmmZqugE48yRJ6ozh\nSZJUKc39TmC1PUlSZwxPkqRKaQ1PKywYIUnqgOFJklQp9dbw5MyTJKkDhidJUqXsv2zPghGSpPYZ\nniRJlTI2Mbn3e2eeJEmdMDxJkirFZXuSpG4ZniRJlTJmeJIkdcnwJEmqlNbwtNxqe5KkDhieJEmV\nMjZe7HkaGR5iZNiPQUlS+/zUkCRVSn2imHlyyZ4kqVOGJ0lSpTSX7RmeJEmdMjxJkirF8CRJ6pbh\nSZJUKc1S5aOGJ0lShwxPkqRK2TvzZKU9SVKHDE+SpEqpl9X2Viwf6XNPJElPNYYnSVKljLlsT5LU\nJcOTJKlSxixVLknqkuFJklQpFoyQJHXL8CRJqpR9pcrd8yRJ6ozhSZJUGdONBvUJq+1JkrpjeJIk\nVUZzyR64bE+S1DnDkySpMpqzTmDBCElS5wxPkqTKGBs3PEmSumd4kiRVRn2/8GTBCElSZwxPkqTK\nGBuf3Pv9qAUjJEkdMjxJkipjv2V7o4YnSVJnDE+SpMrYr2CEM0+SpA4ZniRJlTFmqXJJ0gIYniRJ\nldEMTzVguTNPkqQOGZ4kSZXRLBixfPkwQ7Van3sjSXqqMTxJkiqjWarcazxJkrpheJIkVcZYWTDC\nYhGSpG4YniRJlTFWb848eYFcSVLnDE+SpMpoliq30p4kqRuGJ0lSZTQLRrjnSZLUDcOTJKkyLBgh\nSVoIw5MkqTKa13katWCEJKkLhidJUmWMjVswQpLUPcOTJKkyLBghSVoIw5MkqTKaBSNWGp4kSV0w\nPEmSKmFyaprJqQbgzJMkqTuGJ0lSJTT3O4EFIyRJ3TE8SZIqod4SniwYIUnqhuFJklQJzf1O4HWe\nJEndMTxJkiphbKJ15snwJEnqnOFJklQJrcv2LBghSeqG4UmSVAlj4848SZIWxvAkSaoEC0ZIkhbK\n8CRJqoTWghGWKpckdcPwJEmqhGbBiOGhGstG/PiTJHXOTw9JUiWM1Yvw5H4nSVK3DE+SpEqoTxie\nJEkLY3iSJFVCs9reqMUiJEldMjxJkiqhWTDCmSdJUrcMT5KkSmiWKrfSniSpW4YnSVIlNJftOfMk\nSeqW4UmSVAljFoyQJC2Q4UmSVAkWjJAkLZThSZJUCXULRkiSFqjjX79FxHrgo8BLgB3AdZl56Rxt\n3w68FTgGuAu4JDPvKI+NAn8KvAYYBb4CvCUzt5bHp4E60ABq5e1VmXlxp32WJGnvdZ4sGCFJ6lI3\nM0/XA98Hng68ArggIi6Z2SgizgfeD7wROBq4AbghIlaWTf4AOA14MXBK2ZerW16iAZySmYdk5sry\n1uAkSepYo9FoWbZneJIkdaej8BQRZwAvAN6dmTsz80Hgw8BFszS/CLg6MzdmZh24jCIQnR8Rw8CF\nwO9l5iOZ+WPgvcDPR8Qx5fNr5ZckSQsyPjlNo1F877I9SVK3Op15Oh3YlJnbWx67A4iIWDWj7Yby\nGACZ2QDuBM4EngUcBvxLy/EE9pTPa/pgRHw3IrZGxJWznEOSpHmN1Sf3fu/MkySpW52Gp7XAthmP\nbS1v17XZdl15rDHL8W0tr3MbcBNwEvCTFHus/qLD/kqStHfJHsAKq+1JkrrUzSdIJ0vp5ms75/HM\nPLv1bkS8G/hcRLwpMyfaOfnwsMUEq6z583ccVJdjQM2f/cRUY+9jq1aMMDLimKgK3wfkGFAvf/ad\nhqctFLNGrZqzSFvabHt3eaxW3t/dcvwI4PE5zr0JGAaOAn7QTmfXrFk5fyMNPMeBHAMaGtm3VO+o\ndYdy+OGuAq8a3wfkGFAvdBqeNgLrI+KIZklx4CzgvszcPUvbDcAnACJiiGLP1FXAQxRL9DZQVO4j\nIk4FlgMbI+JFwBsz87daXu+5FKXLH2m3s9u372FqarqzP6EGxvDwEGvWrHQcVJhjQM0xsPXH+z6i\nxuvjbNu2q4+90sHk+4AcA2qOgV7oKDxl5p0RcTvwgYh4J3A88A6KSnpExP3AhZl5K3A5cG1EXEtx\njad3AWPAFzNzOiL+EnhvRGykKBTxB8BnMnNLRCwDLoqIx4GPUJRF/z3gyrLwRFumpqaZnPQvSdU5\nDuQY0O49+1Z7jwwNOR4qyPcBOQbUC90sAHwdRWjaDNwMXJOZV5THTgZWA2TmjcB7gE8DPwLOBc4r\ny5YDvA/4OvBN4EHgCeBN5XMfAc4DXgv8ELgF+CLw7i76K0mquP0LRlhtT5LUnY4LRpTB5jVzHBue\ncf9K4Mo52k4Abyu/Zjt+C3D2bMckSerE2HhLqfJlhidJUncsOyJJGnjNmafly4YYGvL665Kk7hie\nJEkDr16GpxXOOkmSFsDwJEkaeHvKZXteIFeStBCGJ0nSwGvOPI1aLEKStACGJ0nSwBszPEmSesDw\nJEkaeM3wZJlySdJCGJ4kSQOvWarcghGSpIUwPEmSBt6+mScLRkiSumd4kiQNPAtGSJJ6wfAkSRp4\n+0qVG54kSd0zPEmSBl7dghGSpB4wPEmSBtrU1DTjk9MAjFowQpK0AIYnSdJAaxaLAAtGSJIWxvAk\nSRpoe+qTe7932Z4kaSEMT5KkgWZ4kiT1iuFJkjTQmhfIBUuVS5IWxvAkSRpo+888uedJktQ9w5Mk\naaCN1fcVjHDmSZK0EIYnSdJA290682SpcknSAhieJEkDbcyCEZKkHjE8SZIGWnPPU60Gy0b82JMk\ndc9PEUnSQGvOPK1YPkKtVutzbyRJT2WGJ0nSQNszXhSMcMmeJGmhDE+SpIG2Z+/Mk+FJkrQwhidJ\n0kBrLtsbtdKeJGmBDE+SpIHmzJMkqVcMT5KkgbanpWCEJEkLYXiSJA20ZngadeZJkrRAhidJ0kAb\nG3fZniSpNwxPkqSBtqdelCq3YIQkaaEMT5KkgWbBCElSrxieJEkDq9Fo7C1VbsEISdJCGZ4kSQNr\nYmqaqekGYMEISdLCGZ4kSQOrPj6193uX7UmSFsrwJEkaWGP1lvBkwQhJ0gIZniRJA2tswpknSVLv\nGJ4kSQOreY0ngFELRkiSFsjwJEkaWPst23PmSZK0QIYnSdLActmeJKmXDE+SpIHVvMYTWKpckrRw\nhidJ0sAwj7n0AAALOUlEQVSqt8w8jVptT5K0QIYnSdLA2lPueVo2PMTIsB95kqSF8ZNEkjSw6mW1\nvRWjzjpJkhbO8CRJGljNghEu2ZMk9YLhSZI0sJqlyq20J0nqBcOTJGlgNWeeVox6gVxJ0sIZniRJ\nA6tZqnyFy/YkST1geJIkDaz6eHPmyfAkSVo4w5MkaWDtKxjhsj1J0sIZniRJA6u5bG+lM0+SpB4w\nPEmSBtbemSer7UmSesDwJEkaWHtLlVswQpLUA4YnSdJAmm40qFuqXJLUQ4YnSdJAalbaA2eeJEm9\nYXiSJA2k5qwTOPMkSeoNw5MkaSDtN/NkwQhJUg/4qzhJ0pI3MTnNTbd/j02P7mj7ObvLMuVgeJIk\n9YbhSZK0pD2+bTeXf/Zevru5/eA00+pDlvewR5KkqjI8SZKWrH/+1mNc83/uZ6xcgnfculWsXtH+\nR1etVuOFcRTHr1vF5OT0YnVTklQRhidJ0pIzPjHFtX//bb565yMADNVq/JuffiavfvF6hmq1tl9n\nZGSIww9fxbZtuxarq5KkCjE8SZKWlEd+uIsrPnsPD28pAs/aNaO8+RdO5aQTDutzzyRJVWd4kiQt\nGV+7+1E+cVMyPlEssTvt5HX8+nnPYfXKZX3umSRJhidJ0hIwNj7JJ258gNvu3QzA8FCNf/vyk3jF\nhhOodbBMT5KkxWR4kiT11fce28EVn72XzVt3A3DU01by5tc+j2ccu6bPPZMkaX+GJ0lSXzQaDb5y\n5yNc++VvMzlVLNM76zlH8WuvfjYrR/14kiQtPX46SdKAajQa7Ng9wXSj0e+uPMnk5DSf/sqDbLz/\ncQCWjQzx+leczDkvPM5lepKkJcvwJEkD5sc769x6z2ZuuevRvUvhlrJj1x7Cf3rtqZxw1Op+d0WS\npAPqODxFxHrgo8BLgB3AdZl56Rxt3w68FTgGuAu4JDPvKI+NAn8KvAYYBb4CvCUzt3Z6Hkmqusmp\nab75nR9yy12PcvdDW5fkbNNsXvr8Y3nDK09hdPlwv7siSdK8upl5uh64Hfhl4GjgixGxOTM/0too\nIs4H3g+8CrgbuBi4ISKelZl7gD8ATgNeDOwGPgZcDby2k/NIUpU9/PhObrn7UW69ZzM790zsd+yU\nEw7jxc89mlVLtMz3usNW8szjLAohSXrq6Cg8RcQZwAuAl2fmTmBnRHyYIhjNDDUXAVdn5sbyuZeV\n7c6PiM8AFwJvzMxHyuPvBe6LiGOAEzo4jyRVyu6xCf7pvsf4x7seZdPmHfsde9rq5Zz9/GN56fOP\n5egjDulTDyVJGkydzjydDmzKzO0tj90BRESsysxdLY9vAK5t3snMRkTcCZwJ3AkcBvxLy/GMiD3l\n847v4Dw6SBqNBg2ABjRo0GhQfhWPNxr7HoMG0/sdK/4z3Wh5rfJ1aMB08eDe128uOSraNF+gM8Mj\nQ2yvT7F9+x6mJqcX+KfXU9GgjYFtO+p87Z7NfCO37K1OB8U1kU475Uhe+vxjOfUZRzA0ZMEFSZIW\nQ6fhaS2wbcZjW8vbdcCuNtquK481Zjm+reV4u+eZ1a/9lxsZq0/O12w/5T/dF12/tiLMPO+T/rx7\ng8rs4UjS0rL+6NWc88Lj+MlTj+HQQ5b3uztL0vDw0H63qh7HgBwD6uXPvps9T538SnO+tgc6vqBf\nnX78/a/yV6+SJADWrFnZ7y6ozxwDcgyoFzqNYVsoZoVaNWeRtrTZ9vHyWG2W40e0HG/3PJIkSZK0\n6DoNTxuB9RFxRMtjZwH3ZebMi4lspNi/BEBEDFHsmfo68BDFsrzW46cCy8vndXIeSZIkSVp0tUaH\nG3Ai4lbgHuCdFIUdvgBclplXRMT9wIWZeWtEvIqiYMS/orjG07soKuxFZtYj4g+BVwAXAHsoypTv\nzsxfnu88C/wzS5IkSVLHutk99TqKMLMZuBm4piXQnAysBsjMG4H3AJ8GfgScC5yXmfWy7fsoZqG+\nCTwIPAG8qc3zSJIkSdJB1fHMkyRJkiRVkTUbJUmSJKkNhidJkiRJaoPhSZIkSZLaYHiSJEmSpDYY\nniRJkiSpDYYnSZIkSWrDSL870GsRsR74KPASYAdwXWZe2t9eabGVF2X+OHBzZr5+xrGXA38IPBv4\nHvCHmfmpg99LLaby7/5HgHOACeDvgIszc7tjoBoi4oXAh4AzKC6+/lXg7Zn5uGOgeiLiTyjeA4bK\n+46BCoiIaaAONIBaeXtVZl7sGKiOiHgv8BvAocBtwJsy87u9GAODOPN0PfB94OnAK4ALIuKSvvZI\niyoi3kXxj+YHZjl2DPBZikB9JHAJcFVEnH5QO6mD4fPAVuBEYAPwPOCPHQPVEBHLgRspLqp+JHAq\ncDRwuWOgeiLiRcC/p/iHMxFxLI6BqmgAp2TmIZm5sry92PeB6oiI3wBeT/HL1GOB+4B39GoMDNTM\nU0ScAbwAeHlm7gR2RsSHgYsp/nGtwbQHOAv4M2B0xrE3AJmZHy/v/31EfA74j8BbD14XtZgi4jDg\nduA9mbkH2BMRHwfehmOgKg4Bfhu4JjOngR9FxPXAf8YxUCkRUQMup5iF/K/lw46B6qiVXzM5Bqrj\nN4HfzMzvlPcvAYiId9KDMTBQ4Qk4HdiUmdtbHrsDiIhYlZm7+tQvLaLM/HOAiJjt8AaKMdDqDuDf\nLnK3dBBl5hMUb36tTgR+gGOgEjLzx8BfNe9H8YbwH4DrcAxUzVsofqn2KfaFp9NxDFTJByPip4A1\nFO8B78T3gUqIiOOAZwBrI+JeihUIN1OEo56MgUFbtrcW2Dbjsa3l7bqD3BctDXONCcfDACtnof8z\n8N9wDFRKRKyPiDpwL/BPwO/iGKiMiDia4mf+n2YccgxUx23ATcBJFPvfX0KxTMsxUA0nlLevA15O\nsSLtROAqejQGBi08wexTtao2x0SFRMTZFHtf3p2ZN5cPOwYqIjO/l5mjQJRfnygPOQaq4UPA/8jM\nnOWYY6ACMvPszLw6MyfKcXApxf6XERwDVdD8GX8wMx/LzEeA9wO/wL4iIgsyaOFpC0WqbLWW4n/W\nloPfHS0Bc42Jx/vQFy2yiDgf+AJFhbW/KB92DFRQZj4IvBf4FWAcx8DAi4hzgZ8Cfr98qPUfSb4P\nVNcmYBiYxjFQBZvL2ydaHttE8X6wjB6MgUELTxuB9RFxRMtjZwH3ZebuPvVJ/bWRYo1rqzMplvNo\ngJTr268BfjEzP9lyyDFQARHxsxFx/4yHG+XXP1OUL2/lGBg8bwCOAr4XEVuAbwC1iHgcuBvHwMCL\niBdFxB/PePi5wBjwRRwDVfAwsB14Uctjz6D4JVpPxkCt0WgspINLTkTcCtxDsTnweIrfQl+WmVf0\ntWNadBFxNTDaep2niDgS+DZF5ZVPAucCnwZenJn39qWj6rmIGAbuAv4kMz8245hjoAIiYg1wP8Uy\nvd8FVlNc+20l8EvAd3AMDLSy6uaqlodOpNj/cjzFkq27cQwMtLJYwP0UhUI+QnHZmuuBLwF/gJ8F\nlRARH6JYpvdqimu+Xg98i6Ii64LHwCCGp+MoNoX9DMWU3eWZ+fsHfJKe0iJiD8Vvl5eVD00Cjcw8\npDz+UuC/U1wQbRNwaWZ+tg9d1SIpf8ZfpbgwYvOiiM3bAH4Cx8DAi4jnAX9O8ZvEnRQVlt6ZmY/6\nPlA9EfETwEOZOVzedwxUQPlz/iDwfIoZp2uA38nMccdANZTX/fsQ+/a6/S/gbZm5uxdjYODCkyRJ\nkiQthkHb8yRJkiRJi8LwJEmSJEltMDxJkiRJUhsMT5IkSZLUBsOTJEmSJLXB8CRJkiRJbTA8SZIk\nSVIbDE+SJEmS1AbDkyRJkiS1wfAkSZIkSW0wPEmSJElSG/4/tUwbzI4AhL0AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f88b2817390>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(10, 4))\n",
"plt.plot(likelyhood_array)\n",
"plt.title('Likelihood')"
]
},
{
"cell_type": "code",
"execution_count": 127,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"# features: 1, evidence: 3.051457\n",
"# features: 2, evidence: 3.054281\n",
"# features: 3, evidence: 18.526575\n"
]
},
{
"ename": "KeyboardInterrupt",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m\u001b[0m",
"\u001b[0;31mKeyboardInterrupt\u001b[0mTraceback (most recent call last)",
"\u001b[0;32m<ipython-input-127-b881e626daea>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0mw_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw_\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[0mpred_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mw_\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw_\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 25\u001b[0;31m \u001b[0mevidence\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mdata_likelyhood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlabels\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpred_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcov_data\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mdata_likelyhood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw_\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mweights_opt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcov_params\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 26\u001b[0m \u001b[0;32mprint\u001b[0m \u001b[0;34m'# features: %d, evidence: %f'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mn_features\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mevidence\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 27\u001b[0m \u001b[0mevidence_array\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mevidence\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-51-204906cd80f1>\u001b[0m in \u001b[0;36mdata_likelyhood\u001b[0;34m(y, ypred, cov_data)\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mdata_likelyhood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mypred\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcov_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0merr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mypred\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0;36m1e8\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcov_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2.\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m/home/a/miniconda2/lib/python2.7/site-packages/numpy/linalg/linalg.pyc\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m 524\u001b[0m \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'D->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 525\u001b[0m \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 526\u001b[0;31m \u001b[0mainv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 527\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 528\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mKeyboardInterrupt\u001b[0m: "
]
}
],
"source": [
"from itertools import product\n",
"\n",
"likelyhood_array = []\n",
"evidence_array = []\n",
"width = 5\n",
"num = 11\n",
"labels = y\n",
"\n",
"for n_features in range(1,X.shape[1]):\n",
" cov_params = np.eye(n_features + 1)\n",
" evidence = 0\n",
" data = X[:,:n_features]\n",
" lr.fit(data, labels)\n",
" weights_opt = np.append(np.array([lr.intercept_]), lr.coef_)\n",
" \n",
" w_array = np.zeros([len(weights_opt), num])\n",
" for i in range(len(weights_opt)):\n",
" w_array[i, :] = np.linspace(weights_opt[i] - width0, weights_opt[i] + width0, num)\n",
" \n",
" predicted = lr.predict(data)\n",
" likelyhood_array.append(data_likelyhood(labels, predicted, cov_data))\n",
" for w_ in product(*[x for x in w_array]):\n",
" w_ = np.array(w_)\n",
" pred_ = w_[0] + np.asarray(data).dot(w_[1:])\n",
" evidence += data_likelyhood(labels, pred_, cov_data) * \\\n",
" data_likelyhood(w_.reshape(-1, 1), weights_opt.reshape(-1, 1), cov_params)\n",
" print '# features: %d, evidence: %f' % (n_features, evidence)\n",
" evidence_array.append(evidence)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"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.12"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment