{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\tHiroshi TAKEMOTO\n", "\t(take.pwave@gmail.com)\n", "\t\n", "\t

Sageでグラフを再現してみよう：データ解析のための統計モデリング入門第４章

\n", "\t

\n", "\t\tこの企画は、雑誌や教科書にでているグラフをSageで再現し、\n", "\t\tグラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。\n", "\t

\n", "\t

\n", "\t\t前回に続き、データ解析のための統計モデリング入門\n", "\t\t（以下、久保本と書きます）\n", "\t\tの第４章の例題をSageを使って再現してみます。\n", "\t

\n", "\t

\n", "\t\t4章のメインは、久保本図4.9（以下に引用）のバイアス補正とその分布にあると思います。\t\t\n", "\t

\n", "\t

\n", "\t\t\n", "\t

\n", "\t

\n", "\t\t数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。\n", "\t\tこの機会にSageを分析に活用してみてはいかがでしょう。\n", "\t

\n", "\n", "\n", "\n", "\t

前準備

\n", "\t

\n", "\t\t最初に必要なライブラリーやパッケージをロードしておきます。jupyter用に新しくなったRUtil.pyも使います。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%HTML\n", "" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# RとPandasのデータフレームを相互に変換する関数を読み込む\n", "# Rの必要なライブラリ\n", "r('library(ggplot2)')\n", "r('library(jsonlite)')\n", "\n", "# python用のパッケージ\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt \n", "import seaborn as sns\n", "%matplotlib inline\n", "\n", "# jupyter用のdisplayメソッド\n", "from IPython.display import display, Latex, HTML, Math, JSON\n", "# sageユーティリティ\n", "load('script/sage_util.py')\n", "# Rユーティリティ\n", "load('script/RUtil.py')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

あてはまりの良さとモデルの複雑さ

\n", "\t

\n", "\t\t「あてはまりの良さとモデルの複雑さ」を説明するために、３章のデータを使ってk=1とk=7の結果の比較が示されています。\n", "\t\tこのような図を簡単に再現できれば、実際の分析の時にも役立つと思います。\n", "\t

\n", "\t

\n", "\t\t3章のデータをも一度読み込み、statsmodelsのglmを使ってk=1とその結果をプロットします。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 3章のデータをもう一度読み込む\n", "d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv')\n", "d_C = d[d['f'] == 'C']\n", "d_T = d[d['f'] == 'T']" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# statsmodelsを使ってglmを計算します\n", "import statsmodels.formula.api as smf\n", "import statsmodels.api as sm\n", "from scipy.stats.stats import pearsonr" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# k=1でのglm回帰を実行\n", "fit_1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "# λの予測値をλ_1にセット\n", "d['lam_1'] = fit_1.predict()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAESCAYAAAD67L7dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+clGW9//HXLAsIArsDbOoIDZh2dczI7GEetE0kHxgd\nS82sHl+tEKMOoB7j9P2SnY6ldvScPamkKH2PPyJ5nLIsfyQJlEZ97YTFkQxKvSxl96wOKTC7C8qv\nhZ3vHzOLs+PO7vy672tmrvfzH2fnx3V/ruu+/ezNPTPvjaRSKURExB8NrgsQEZFwqfGLiHhGjV9E\nxDNq/CIinlHjFxHxjBq/iIhnGoPegDHmJOAh4GZr7R3GmEbgu8DxwC7g49banqDrEBGRtEDP+I0x\nY4Fbgcey7l4AvGqtPQ34AdAaZA0iIjJQ0Gf8+4C5wJez7vsIcA2AtfaugLcvIiI5Am381to+YL8x\nJvvuacCHjTH/DmwDFllru4OsQ0RE3uDizd0I8Ky19izgT8BXHNQgIuKtwN/cHcRfgf+Xub0O+PpQ\nTz548FCqsXFE0DWJiNSbSL4HXDT+NaSv+68E3gvYoZ7c1bUnhJLC1dIynu3bd7suwxnNX/P3ef4Q\nzhq0tIzP+1igjd8YcwpwExAHeo0xHwf+F3CrMeYyYDfw2SBrEBGRgYJ+c3cTcNYgD30iyO2KiEh+\n+uauiIhn1PhFRDyjxi8i4hk1fhERz6jxi4h4Ro2/RGvWrOb2279V8XHvv/8+Zs36W/bt21fxsUVE\nQI2/LJG834srzdq1P6WrK0lLy1sqO7CISBYX39ytG6kU3HbbLTz77J/o7T3AeeddyLnnnscNN1xL\nc3MUa5+ju7uLiy/+LI8++hN27eph+fL/AAb/Rt2ZZ85mzJgx/Pzna8OdiIh4peYb/w9/8Rc2Pvdq\nRcc89R1v4ROzjy/oubFYjCuu+CL79+/nk588n3PPPQ+AxsZGvvWtO7juun/mT3/azLJld3D99dew\nadN/E4+fO+hYY8aMqdgcpHw9ySRPLF1CU0c7PfE4rW230BSd6LoskbLVfON3raenh4UL59PYOJKe\nnjfSpU888Z0ATJo0mXh8GgATJ07itddec1GmlOCJpUuY9/ADRIDU05tYSYRz71zpuiyRstV84//E\n7OMLPjuvNGufpa+vj9tvv4uGhgbmzDnz8GMjRowY9HYqlSpg5Aq/eSAlaepoP7wnIpmfReqB3twt\nw7Zt2zjqqKNoaGjg17/+FX19hzh48GAFRi7kl4MErSceP7wnUkBP5l9uIrWu5s/4XfrAB2axZcsf\nuOKKL9DaOovTT2/lppv+dcBzsj/5M9yngO699x42bvwtyeROvvSlK3nnO9/FwoVXBFC5FKK17RZW\nEslc459Ga9vNrksSqYhIYZce3Nm+fXd1F1gC3/PINX/N3+f5Q2h5/FX1h1i89sorf2XJkkUcPNgH\npK/5RyIRTj75FObP/7zj6kTEB2r8ITvqqKNZtWqV92c8IuKO3twVEfGMGr+IiGfU+EVEPBN44zfG\nnGSM+YsxZlHO/ecYY/qC3r6IiAwU6Ju7xpixwK3AYzn3jwa+DCSC3H6Qli9fhrXPkkzuZO/evUyZ\nMpUJEybwjW+0uS5NRGRIQX+qZx8wl3STz/YVYDnw7wFvPzCXX34VkM7l37r1BRYt+gfHFYmIFCbQ\nSz3W2j5r7f7s+4wxbwdmWGt/TIihNH9Y/zg/v+0Wnvvdk2FtUqRoPckkqxfM44k5s1i94LP0dCVd\nlyR1yMXn+G8GQs0h+OWK5byr7V84+/XX+X1zlN/9Sxvvu+iTYZYgUhAlgkoYQm38xpgYYID/NMZE\ngGOMMeuttWfle000OpbGxhH5Hi7IyIfu5x2vvw7Ae7q76Pjx92lZ9Lmyxuw3fvwRjBkzipaWwf+4\nSj7FPr/eaP6Dz39yonNAIujkRGddrlU9zqlYLtcgzMYfsdYmgBP67zDGbB2q6QN0de0pe8MHDg38\n8NC+g30V++bs7t372Lv3QFHj+Z5Vovnnn/+O2BRSbEyf8QM7YlPrbq183/8QWlZP3seC/lTPKcBN\nQBzoNcZcCHzMWtv/F0tCCWAbdck8nv7G13j3rl1saGmh+TPzw9isSNGUCCphCLTxW2s3AXnP6K21\nxwW5/X5nzLsM+64ZfG/TU5xwRivvzfx1LJFq0xSdqGv6EjhvQtrMe0/FvPfUio87d+7gfz9XRKRa\nKbJBRMQzavwiIp5R4xcR8Ywav4iIZ9T4RUQ8o8YvIuIZbz7OGYSXXurk1ltvoru7m76+Pk46aQaL\nF/8DI0eOdF2aV3qSSZ5YuiTzpac4rW230BSd6Los8UwtHYdq/CXq6+vjn/7p/7BkyVLe/e6TAVi2\n7JusXHkXCxYsdFydXxRsJtWglo5Dbxr/+vW/549//CunnRbnfe87sezxNm78LdOmTTvc9AEWLbqS\nhgZdPQtbU0f7gGCzpo52h9WIr2rpOPSiS61Y8Tjz5x/F9dd/gksuGcX9928oe8yOjnaOP94MuG/U\nqFE0Nnrzu7Rq9MTjh0OfUkBPfJrDasRXtXQcetGlfvSjQ7z+erpJd3e/h/vu+xEXXVTemJFIhL6+\nQxWoTsqlYDOpBrV0HHrR+COR1JA/lyIen8aPf/yDAff19vbS2fk/HHfc28oeXwqnYDOpBrV0HHpx\nqeeSS45gwoSngRQtLRv4zGfKf6f91FNP45VXXuE3v/k1kH6zd8WKW1m//rFhXiki4lYklQolEr9k\n27fvrkiBTz31LJs2tXPGGYYTT6xMGnQyuZN/+7dvkEzupLFxJKeeehrz539+2Nf5/ocoNH/N3+f5\nQ2h/iCXv3zT3pvFXE98PfM1f8/d5/uC+8XtxqUdERN6gxi8i4hk1fhERzwT+cU5jzEnAQ8DN1to7\njDFTgXuAkcAB4BJr7atB1yEiImmBnvEbY8YCtwLZn3G8Hvi2tXYW6V8I/xhkDSIiMlDQZ/z7gLnA\nl7PuW5i5H2A78J6Aa5AqlUx2s3Tpejo6JhCP99DWNptotNl1WVJhuamVF9xzF+l/8JfP5TFUy8dv\noI3fWtsH7DfGZN+3F8AY0wAsBq4NsgapXkuXrufhhz8NRHj66RSwijvvvMB1WVJhuamV9y0cydnL\n76rI2C6PoVo+fp1ENmSa/irgcWvt+qGeG42OpbFxRDiFhailZbzrEpxqaRlPIhGFrDzDRCLqzbr4\nMk+AyYnOAamV47Zurdj8XR5D5W7b5THgKqvnO4C11l4/3BO7uvaEUE64fP8CS//8Y7Ek6RzDCJAi\nFuvyYl182/87YlNIsTGzl+G16dMrNn+Xx1A52w7pC1x5Hwu98RtjLgb2W2uvC3vbUl3a2mYDqzLX\nSHfR1naW65IkALmplResWEFvhYJtXR5DtXz8BhrZYIw5BbgJiAO9wMvAW0i/ubub9K/LZ6y1l+cb\nQ5EN9Ufz1/x9nj+4j2wI+s3dTUDt/BoUEfGAvrkrIuIZNX4REc+o8YuIeEaNX0TEM2r8IiKeUeMX\nEfGMGr+IiGdcRTaIBCaZ7Oaqq9bx5JMNwA5mzhzHsmUfKSs5MTdhsrXtFpqiEytXdACCSo/sSSZZ\n88UruW9DlAQz+JuZ41m27EOhJVPW4r6oNmr8UneWLl3P2rWX0Z+hsmbN9xk1an1ZyYm5CZMriXDu\nnSsrVHEwgkqPfGLpEp5YA7/hO0CE9jUpRo0KL5myFvdFtdGlHqk7HR0TYEAe5PjMfaVr6mgfMGJT\nR3tZ44Uhdx3KXYN+TR3ttDMjkLEL3X6t7Ytqo8YvdSce7yEdA0Xmv7uJx3eVNWZPPD5gxJ74tLLG\nC0PuOpS7Bv164nGmsTmQsQvdfq3ti2qjSz1Sd9raZnPgwN1s2NAA7GTmzHG0tZ1b1pi5CZOtbTdX\nptgABZUe2dp2C68duIKXN1xKghmcOHMCbW3nVGTsQrdfa/ui2gSazlkJSuesP5q/5u/z/MF9Oqcu\n9YiIeEaNX0TEM2r8IiKeUeMXEfGMGr+IiGfU+EVEPBP45/iNMScBDwE3W2vvMMZMAVaR/qWzDfi0\ntbY36DpERCQt0DN+Y8xY4Fbgsay7rwNus9aeCbwAzA+yBhERGSjoSz37gLmkz+z7zQIeydx+BDg7\n4BokYMlkNwsWPMicOY+zYMEDdHV1H37shRc6OPnk24jHH+Tkk29l69aOgsftSSZZvWAeT8yZxeoF\nn6WnKxlE+QUZao6VHnvTpi0Fr1mhdeU+r+PFrVWzti5U07HlQqCXeqy1fcB+Y0z23UdmXdp5FTgm\nyBokeEOlQF544U9IJK4GIuzdm+KCC27k5ZdPKmjcakphDCrpcrCx1627hn37riN7zZ5++oqy6sp9\nXmLj5fw2UR1r60I1HVsuuM7qyfuV4n7R6FgaG0eEUUuoWlrGuy6hYhKJKNlJjYlE9PD8urunDHgs\n/XNh85+c6ByQwjg50els3YaaYymyX5s79v7908lds3zbKrSu3Od1d091urauj/9qOLZcroGLxr/b\nGDPaWrsfOBZIDPXkrq494VQVonrLKonFkqRzEtP597FY1+H5NTd3smfPG481N78EUND8d8SmkGJj\n5pWwIzbV2boNNcdi5e7/3LFHj36RffsGrlm+bRVaV+7zmps7Se3BydpWw/Hv+tgKKasn72MuGv9j\nwIXA9zL/XeugBqmgoVIgH3zwPC644Ea6uqYQjb7Egw9+tOBxqymFMaiky8HGXrz4w8ybV9iaFVpX\n7vO++pUvsPKGv1bF2rpQTceWC4GmcxpjTgFuAuJAL/AycDHwXWA00AFcaq09lG8MpXPWH81f8/d5\n/uA+nTPoN3c3AYOdgswJcrsiIpKfvrkrIuIZNX4REc+o8YuIeEaNX0TEM2r8IiKeUeMXEfGMGr+I\niGdcZ/VIjUgmu1ly1TqeeXIXMTbzqZk7mbtsOU3Ria5LK0pPMskTS5dkvrEZp7XtlrLnkG/MZLKb\npUvX8+ILYzkiuZ7FE/+LvrcdxwX33AWMPPx4+tu0PbS1zSYabR52vOznp1LpALYXXhhBMtnBpElv\nZ8qxXcyNPML4zk7uSJ7Ovklncdxxrx8efzD5apHiBHF8BUGNXwqydOl6Hl17GRChnRTHrjmfcaOW\n1FyiYRCpjPnGzE7EhAt5NHE+9/3xQe5bOJKzl9+VN1mzkPH6nw9k7rsPuJpEIsKWLSnGsZoIb2Uj\nyyFz31CJokGmj/qkVlI/1filIB0dE8hOd2xnBk0dj7osqSRNHe0DUhmbOtoDGzN3zbYygwg/YdzW\nrYM+nv658PH6n5++bxy5++eNx3Kf/2b5x5ZiBHF8BUHX+KUg8XgP6RxDgBTT2ExPfJrDikrTE49n\nzYKKzCHfmLlrNp3NpIDXpk8f9PF4fFdR48Xju7Lu203u/pnO5kHHH0y+WqQ4QRxfQdAZvxSkrW02\nvQfu5pkN6Wv8rTNTNZloGEQqY74x+xMxX3xxLEfsXM+HJ3aw8m0f44IVK+g9lD9Zc7jx3pzEuSpz\njf9GJk16O1OndNNKigmd/8PW5OWZa/x7hkwUDTJ91Ce1kvoZaDpnJSids/5o/pq/z/MH9+mcutQj\nIuIZNX4REc+o8YuIeEaNX0TEM2r8IiKeUeMXEfFM6J/jN8YcCdwLRIFRwHXW2p+FXYeIiK9cnPHP\nA56z1s4GLgK+5aAGERFvufjm7g7gXZnbE4HtDmqQIoWV3lgr6YZh61+XyAsvFpy46VruMfPVq9/D\nlhuvZXKikx2xKQP2bSnHVznHSqmvLeZ1rhJP+2tMPfzA7+anUu8b7DnDNn5jzIestWsrVZS19gfG\nmHnGmD8DzcDfVWpsCU5Y6Y21km4Ytv51+RQfLThx07XcYyax8XJ+m8jsWzYO2LelHF/lHCulvraY\n17lKPM2q8dR8zynkjP9KY8xy4D+Be6y1HeUUZYy5GOiw1s41xswA7maIAqPRsTQ2jihnk1WppWW8\n6xKKkkhEyU5vTCSiZc0h32snJzoHpBtOTnTW3FoVotg59a/LVmZQyf0QpNxjprt7at59W8rxVc6x\nUupri3ndcHMKar9l15jPsI3fWvthY0wUuABYYYwB+A7wgLX2UAl1nQGsy4y92RgTM8ZErLWDZvJ0\nde0pYRPVrRazSmKxJOm8wfT5WizWVfIchpr/jtgUUmzMbAV2xKbW3FoNp5T9378u09nMxgrth6Dl\nHjPNzZ2k9jDovi3l+CrnWCn1tcW8bqg5BdkDsmvMp+CQNmPMWOBjwCJgBHAk8Dlr7ZPFFGWMWQIc\nZa1daoyJA+uste/I9/x5166ru5C2ESMiHDpUW9Pq6+ujq2sfBw9GaGxMEY0eQUNDaZ8NGGr+fX2H\n2N/VRcPBg/Q1NjI6GqWhob7+xVfK/u9fFw72srvvCFINR5S9H4KWe8w0NY2it6ebhkMH6RsxcN+W\ncnyVc6yU+tpiXjfUnILsAf01nnPbJRvLucb/AeBS4CzgAeAya+2zxphpwIPAe4qs6/8C9xhjfkn6\nF8gXiny9ONDQ0MCkSWND2M4IxkyaHPh2ak32uoxxXEuhBjtmGidNHrTplXJ8lXOslPraYl4X1v8z\nb95uusZ8TR8KOOM3xvwa+DZwv7V2f85jV1trb6xItXkolrn+aP6av8/zB/exzIVc43//EI8F2vRF\nRKTyqvPioIiIBEaNX0TEM2r8IiKeUeMXEfGMGr+IiGfU+EVEPOMinVMkNGEnJPYkk/zsi4t5fcN/\nMQnonXkGH1x2e9EJo8ONE3SKabHjK1W1tqjxS10LOyHxiaVLaF7zUxaSSWhZ81NWjhpddMLocOME\nnWJa7PhKVa0tavxS1zo6JpCdkJj+OThNHe2MHLDF9H2VHqepo73sbQy3/WLGD7oeqSxd45e6Fo/3\nkE5IBEgRj+8KdHs98Ti7BmwReuLTKj5OTzxe9jaG234x4wddj1SWzvilrrW1zQZWZa7x76Kt7axA\nt9fadgvrDhzgm5lr8wdnvp/ZbTdXfJzWtltYSSRzTX0arSVsY7jtFzN+0PVIZRUcy+yKQtrqj+av\n+fs8f3Af0qZLPSIinlHjFxHxjBq/iIhn1PhFRDyjxi8i4hk1fhERzzj5HL8x5mLgfwO9wDXW2jUu\n6hAR8VHoZ/zGmInANcDpwLnAeWHXICLiMxdn/GcDP7fW7gH2AH/voIaaV2rqZLWkKJZSx3Bzzn78\nmKO309r7EKufaiHBDP5m5ni+9rWZ3HjjptCSOvPNs6VlfFFjJJPdLLlqHc88uYsYm/nUzJ3MXbZ8\nwHoNtTa56zI38ghHd3bwfHInb504kYNvO77ix0GlUlHDPl7DSnPt3rmT1Qs+X9F55a7VVx9+MPpy\nKtU12HNdNP5pwJHGmIeBZuBaa+0vHNRR00pNnayWFMVS6hhuzgMeJ8UWfs9L3A5EaF+T4g9/uJFE\n4uq8rw/CYPM8/qEfFzXG0qXreXTtZUCEdlIcu+Z8xo1aMmC9hlqb3HUZx2p+yGZSwH2Jl7nkj1sq\nfhxUKhU17OM1rDTXNYsWVXxeuWv1PKwAPjXYc100/ggwETgfmA6sB+L5nhyNjqWxcURIpYWn2LO+\nXIlElOzsxkQiWtCYkxOdA1IUJyc6y66lFKXUMdyccx9P8q4BP3d3Txny9UEYbJ5Q3P7PnVc7M5ic\nWDfk3LPnNtjr4SdEgHEEcxwMt68K3VbYx2up/18Va9zWrRWfV+5anQDH5Xuui8b/CvAba20KeNEY\ns9sYM9lau2OwJ3d17Qm3uhBUIqcjFkuSzkGMAClisa6CxtwRm0KKjZlXwY7Y1NBzU1paxpdUx3Bz\nzn18IlvYk/Vzc3Mne/YUv2blGGyeQFHbzZ3XNDa/ab2GWpvBXk/mntcI5jgYqp5ijv+wj9dS/78q\n1u7p00ltrOy8ctfqz/Bivue6aPw/A75jjGkjfeZ/ZL6mL/mVmjpZLSmKpdQx3JyzH48ds4P3H2hn\n9VOXkmAGJ86cwNe+fh433BBeUidUZr3b2mbTe+BuntmQvsbfOjP1pnGGWpvcdWklxf0vzeD5nelr\n/CvfdkLFj4NKpaKGfbyGleb64RUrWLn/YEXnlbtW9z78wMIb8zzXSTqnMWYB8DnSv5iut9b+NN9z\nlc5ZfzR/zd/n+YP7dE4nn+O31t4J3Oli2yIivtM3d0VEPKPGLyLiGTV+ERHPqPGLiHhGjV9ExDNq\n/CIinlHjFxHxjJPP8RdjzpzHQ0tSrHfZyYNHH72dSOQg27YdU/b6ZqcCvnzMsazloyS2Ta7a/RZE\n4mMxYxabzFgtiapSP6q+8T/99PmhJSnWu+zkwfSXpr8PlL++2amAn3h6Co+STpKs1v0WROJjMWMW\nm8xYLYmqUj9q5FJPhI6OCa6LqHnpNczO7xt/+HY569vU0X541HTy4xvbqMb9ll1vJPNzmGPmJjMO\nt/0g6hW/1UjjTxGP73JdRM2Lx3tIn+mT+e/uw7fLWd+eePzwqOnkxze2UY37LbveFNATnxbqmLun\nTy9q+0HUK36r+ks9J5/8UGhJivUuO3nwmGN2AL1s21b++manArYek+J17s5c46/O/RZE4mMxYxab\nzFgtiapSP5ykcxZD6Zz1R/PX/H2eP7hP56yRSz0iIlIpavwiIp5R4xcR8Ywav4iIZ9T4RUQ8o8Yv\nIuIZZ43fGHOEMeYvxpjPuKpBRMRHLs/4/xnY6XD7IiJecvLNXWOMAd4B/NTF9utFdtpmtSZh1oLB\n0i8PpRoKWltf9oESQuuLq8iGm4DFwDxH268L2Wmb1ZqEWQsGS798mPMKWltf9oESQutL6I3fGPNp\n4DfW2o70iT95v1YMEI2OpbFxRCi1hamlZfzwTxpGIhElOwkzkYhWZNwwVFOdkxOdA9IvJyc6SVDY\n2pa6D6pp/oUYbI3KmUOtzT8ILtfAxRn/3wHTjTEfAaYA+4wxndbaXwz25K6uPaEWF4ZK5XTEYknS\neY3pfP1YrKsmMlCqLatlR2wKKTYe/isFO2JTiVHY2payD6pt/oUYbI1KnUMtzr/SQsrqyftY6I3f\nWvup/tvGmK8BW/M1fRladtpmtSZh1oLB0i/PoIFC1taXfaCE0PriNJ0zq/Hfm+85SuesP5q/5u/z\n/MF9OqfTPH5r7bUuty8i4iN9c1dExDNq/CIinlHjFxHxjBq/iIhn1PhFRDyjxi8i4hk1fhERzzj9\nHL+PksluLr98Nc8/P6au0xyrgS/JmdVIaZ7VTY0/ZL6kOVYDrbU7SvOsbrrUE7KOjglkpzmmf5Yg\naK3daepoH5Dm2dTR7rAayaXGH7J4vId0viFAinh8l8ty6prW2p2eeDxr5aEnPs1hNZJLl3pC1tY2\nm9Gj78tc46/fNMdq4EtyZjVSmmd1c5rOWQilc9YfzV/z93n+4D6dU5d6REQ8o8YvIuIZNX4REc+o\n8YuIeEaNX0TEM04+zmmMaQPeD4wA/tVa+6CLOkREfBT6Gb8xZhZworX2dGAusCzsGkREfObiUs+v\ngIsyt7uBscaYvJ83FRGRygr9Uo+1NgXszfz4OeDRzH0iFaNkTpH8nEU2GGPOAy4F5riqQeqXkjlF\n8nP15u45wNXAOdbaIb+3HI2OpbFxRDiFhailZbzrEpwKev6JRJTsZM5EIlpVa15Ntbjg+/zB7RqE\n3viNMROANuCD1tqe4Z7f1bUn+KJC5ntWSRjzj8WSpHMhI0CKWKyratZc+9/v+UNoWT15H3Nxxv9J\nYBLww8ybuingM9balxzUInVKyZwi+bl4c/dO4M6wtyt+iUabdU1fJA99c1dExDNq/CIinlHjFxHx\njBq/iIhn1PhFRDyjxi8i4hk1fhERz6jxi4h4Ro1fRMQzavwiIp5R4xcR8Ywav4iIZ9T4RUQ8o8Yv\nIuIZNX4REc+o8YuIeEaNX0TEM2r8IiKeUeMXEfGMiz+2jjHmZuBvgT7gKmvtf7uoQ0TER6Gf8Rtj\nPgAcb609HfgccGvYNYiI+MzFpZ4PAg8BWGufA5qNMeMc1CEi4iUXjf9oYHvWzzsy94mISAiq4c3d\niOsCRER84uLN3QQDz/BjwLZ8T25pGV+XvxhaWsa7LsEpzV/z953LNXBxxv8z4OMAxphTgJetta87\nqENExEuRVCoV+kaNMTcAZwKHgMXW2i2hFyEi4iknjV9ERNyphjd3RUQkRGr8IiKeUeMXEfGMk6we\nHxljjgTuBaLAKOA6a+3P3FYVHmNMBPg2cBKwH/h7a+3zbqsKhzHmJNLfVr/ZWnuHMWYKsIr0idc2\n4NPW2l6XNQYpd/6Z+64Evgk0W2v3uKwvDIMcA1OBe4CRwAHgEmvtq2HVozP+8MwDnrPWzgYuAr7l\ntpzQnQdMsNaeQTqj6SbH9YTCGDOWdB7VY1l3XwfcZq09E3gBmO+itjAMNn9jzKeBtwAvu6orTHmO\ngeuBb1trZ5H+hfCPYdakxh+eHcCkzO2JDIyt8MEJwO8ArLUvAvHMvwLq3T5gLgO/pDgLeCRz+xHg\n7JBrCtNg83/AWvtVR/W4MNgaLAQeyNzeTronhEaNPyTW2h+QbnZ/Bn4JfMltRaHbApxjjGkwxhhg\nOjDZcU2Bs9b2WWv359x9ZNalnVeBY0IuKzSDzd+3L2zmWYO91tqUMaYBWAx8L8ya1PhDYoy5GOiw\n1p5AOqH0dsclhcpau5b0Gf+vgCuBZ1FOE2gNvJVp+quAx62168Pcthp/eM4A1gFYazcDMU8udRxm\nrb3GWttqrV0MTAzzzawqs9sYMzpz+1jS+VU+8v3bo98BrLX2+rA3rMYfnr+Q/qtjGGPiwG5rrTcH\nvjFmhjHm7sztDwFPOS7JpceACzO3LwTWOqwlTLknOl6d+GTLXAHYb629zsX2FdkQkszHOe8BjgJG\nAF+11v7KbVXhyfzr5m7gncBe4GJrbd1/qiMTRHgTEAd6SX+S5WLgu8BooAO41Fp7yFmRAcoz/58D\nc4DTgI3ABmvtl50VGbA8a/AW0m/67ib9L59nrLWXh1WTGr+IiGd0qUdExDNq/CIinlHjFxHxjBq/\niIhn1PhFRDyjxi8i4hk1fhERz6jxi4h4Ro1fpEjGmC8aY/4jc9sYY57NfDNbpCao8YsUbxnwdmPM\n6aRTVhckqxjeAAAAjElEQVT4FjUstU2NX6RImXC9y4AfAputtb92XJJIUdT4RUoziXTA1ltdFyJS\nLDV+kSIZY44AVgAfAQ4YYy5xXJJIUdT4RYp3Lem/G/sX4Crg68aYmOOaRAqmWGYREc/ojF9ExDNq\n/CIinlHjFxHxjBq/iIhn1PhFRDyjxi8i4hk1fhERz6jxi4h45v8DEkYTDQ0aR+kAAAAASUVORK5C\nYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# λの予測値と一緒にプロット\n", "ax = d.sort_values(by=['x']).plot(x='x', y='lam_1')\n", "d_T.plot(kind='scatter', x='x', y='y', c='r', ax=ax, label='T')\n", "d_C.plot(kind='scatter', x='x', y='y', c='b', ax=ax, label='C')\n", "plt.legend(loc='upper left')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

多次元のglmをstatsmodelsで行う

\n", "\t

\n", "\t\t次にk=7（$log \\lambda = \\beta_1 + \\beta_2 x + ... + \\beta_7 x^6$）の計算は、$x^n$\n", " を計算するカスタム関数nPwOfをformulaの中で使用することで、実現できました。\n", "\t

\n", "

\n", " 予測値は、データのある点のみを結線したため、カクカクしているところもありますが、データの中心部分を通っています。\n", "

\n", "

\n", " k=1とk=7の結果をみるとk=7の結果は、あまりにも複雑な曲線を示しており、機械学習でいうところの「過学習」 \n", " （Overfitting）の状態になっています。過学習では、学習用データ（ここでいうサンプルデータ）に引きずられ、 \n", " 真のモデルから離れてしまいます。\n", "

\n", "" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# カスタム関数nPwOfを使って多項式近似を行う\n", "def nPwOf(n, x):\n", " return x^n\n", "formula = 'y ~ nPwOf(7, x) + nPwOf(6, x) + nPwOf(5, x) + nPwOf(4, x) + nPwOf(3, x) + nPwOf(2, x) + x'\n", "fit_1 = smf.glm(formula=formula, data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "# λの予測値をλ_1にセット\n", "d['lam_7'] = fit_1.predict()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAESCAYAAAD67L7dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOW9+PHPJJOF7BMSEoaECesDCKgoIiCKS1EUF2pt\n9WoVF9rr1mu9t6XttYvaam9al1pb+7tuVG+r1ipaF1xQVFpRkVW2hzUhZAIkmWSyrzO/PybEJGaS\nmWRmziTn+369eDGZOfM83+eZk++cnOV7LF6vFyGEEOYRY3QAQgghIksSvxBCmIwkfiGEMBlJ/EII\nYTKS+IUQwmQk8QshhMlYw92BUmo68ArwoNb6j0opK/BnYCJQA3xDa+0OdxxCCCF8wrrFr5RKAh4B\n1nR5ejlwTGs9B3gBWBDOGIQQQnQX7i3+JmAx8KMuz10M/AxAa/1EmPsXQgjRQ1gTv9baAzQrpbo+\nXQBcqJT6DVAG3KK1rg5nHEIIIb5kxMFdC7BLa302sAP4iQExCCGEaYX94G4vjgAfdTx+G/hFXwu3\ntbV7rdbYcMckhBDDjcXfC0Yk/tX49vuvBE4BdF8LV1U1RCCkyMrOTqW8vNboMAwj45fxm3n8EJk5\nyM5O9ftaWBO/UmoW8ADgAFqVUt8A/g14RCl1I1ALXBfOGIQQQnQX7oO7m4Cze3npm+HsVwghhH9y\n5a4QQpiMJH4hhDAZSfxCCGEykviFEMJkJPELIYTJGHEe/7CwevXrHDiwn1tv/Y+QtfnTn/4It7sa\nr9dLTU0N06fP4Ac/kAubhRChJYl/ECx+r4sbmHvv/XXn4/vvv4eLL74stB0IIQSS+AfF64Xf//4h\ndu3aQWtrC5deejlLllzKfffdTUaGDa13U11dxdVXX8ebb/6Dmho3jz76v4D/K+oADh0qpr6+jilT\npkVmIEIIUxnyif9v7+9jw+5jIW1z9pRRfPOciQEta7fbuf3279Pc3My3vnUZS5ZcCoDVauV3v/sj\n99zzU3bs2MbDD/+Re+/9GZs2fY7DsaTPNl988Xkuv/xbgx6HGBy3y8W6FXeSXlyE2+FgQeFDpNsy\njQ5LiEEb8onfaG63m5tvvgGrNQ63+8vq0tOmnQDAyJFZOBwFAGRmjqSurq7P9tra2vjii63853+u\nCFvMIjDrVtzJsldfxgJ4t2xiJRaWPL7S6LCEGLQhn/i/ec7EgLfOQ03rXXg8Hv7whyeIiYlh0aKz\nOl+LjY3t9bHX6+2zzc2bNzJ16gmhD1YELb24qLO8oaXjZyGGAzmdcxDKysrIyckhJiaGf/7zQzye\ndtra2gbV5u7dO5k4cVKIIhSD4XY4OP417QXcHX+5CTHUDfktfiOdeeZCvvhiK7ff/l0WLFjIvHkL\neOCBX3dbpuuZP4GcBVRZWcmYMfkhjlQMxILCh1iJpWMffwELCh80OiQhQsLS364Ho5WX10Z3gANg\n9nrkMn4Zv5nHDxGrxx9VN2IxtaNHj3DnnbfQ1uYBfPv8LRYLJ500ixtu+I7B0QkhzEASf4Tl5OTy\n7LPPmn6LRwhhHDm4K4QQJiOJXwghTEYSvxBCmEzYE79SarpSap9S6pYez5+vlPKEu38hhBDdhfXg\nrlIqCXgEWNPj+QTgR4AznP2H06OPPozWu3C5KmlsbCQvL5+0tDR++ctCo0MTQog+hfusniZgMb4k\n39VPgEeB34S5/7C57bY7AF9d/oMH93PLLaGryy+EEOEU1l09WmuP1rq563NKqcnATK31S0CIK9r7\nt3Xte7z7+4fY/dknkepSiKC5XS5eX76MdYsW8vry63BXuYwOSQxDRpzH/yBweyQ7/OCxR5lR+CvO\nq69nc4aNz35VyGlXSNljEX2kIqiIhIgmfqWUHVDAX5RSFmC0Umqt1vpsf++x2ZKwWmP9vRyQuFde\nZEp9PQAnV1dR/NJzZN9y06DaPC41NZERI+LJzu775io9Bbv8cCPj7338Wc6SbhVBs5wlw3KuhuOY\ngmXkHEQy8Vu01k6gs/SkUupgX0kfoKqqYdAdt7R3P3moqc0Tsitna2ubaGxsCao9s9cqkfH7H3+F\nPQ8vG3xb/ECFPX/YzZXZP3+IWK0ev6+F+6yeWcADgANoVUpdDnxda338jiURKcAWf80ytvzy55xY\nU8P67Gwyrr0hEt0KETSpCCoiIayJX2u9CfC7Ra+1Hh/O/o+bv+xG9IyZ/HXTRibNX8Ap0+RGJyI6\npdsyZZ++CDvTFGlTp8xGnTI75O0uXtz3/XOFECLaSMkGIYQwGUn8QghhMpL4hRDCZCTxCyGEyUji\nF0IIk5HEL4QQJmOa0znD4fDhEh555AGqq6vxeDxMnz6TW2/9D+Li4owOzVTcLhfrVtzZcdGTgwWF\nD5FuyzQ6LGEyQ2k9lMQ/QB6Ph//+7x9y550rOPHEkwB4+OHfsnLlEyxffrPB0ZmLFDYT0WAorYem\nSfxr125m+/YjzJnj4LTTpg26vQ0bPqWgoKAz6QPccsv3iImRvWeRll5c1K2wWXpxkYHRCLMaSuuh\nKbLUY4+9xw035HDvvd/kmmviefHF9YNus7i4iIkTVbfn4uPjsVpN810aNdwOR2fRJy/gdhQYGI0w\nq6G0HpoiS/397+3U1/uSdHX1yTz//N+54orBtWmxWPB42kMQnRgsKWwmosFQWg9NkfgtFm+fPw+E\nw1HASy+90O251tZWSkoOMX78hEG3LwInhc1ENBhK66EpdvVcc00iaWlbAC/Z2eu59trBH2mfPXsO\nR48e5eOP/wn4DvY+9tgjrF27pp93CiGEsSxeb0RK4g9YeXltSALcuHEXmzYVMX++Ytq00FSDdrkq\n+Z//+SUuVyVWaxyzZ8/hhhu+0+/7zH4jChm/jN/M44eI3YjF7z3NTZP4o4nZV3wZv4zfzOMH4xO/\nKXb1CCGE+JIkfiGEMBlJ/EIIYTJhP51TKTUdeAV4UGv9R6VUPvAUEAe0ANdorY+FOw4hhBA+Yd3i\nV0olAY8AXc9xvBf4k9Z6Ib4vhP8MZwxCCCG6C/cWfxOwGPhRl+du7ngeoBw4OcwxiCjlclWzYsVa\niovTcDjcFBaeg82WYXRYIsR6Vq1c+tQT+P7gHzwj16GhvP6GNfFrrT1As1Kq63ONAEqpGOBW4O5w\nxiCi14oVa3n11W8DFrZs8QLP8vjjS40OS4RYz6qVz98cx3mPPhGSto1ch4by+mtIyYaOpP8s8J7W\nem1fy9psSVitsZEJLIKys1ONDsFQ2dmpOJ026FLP0Om0mWZezDJOgCxnSbeqlSkHD4Zs/EauQ4Pt\n28h1wKhaPU8DWmt9b38LVlU1RCCcyDL7BSzHx2+3u/DVMbQAXuz2KlPMi9k+/wp7Hl42dHzKUDdu\nXMjGb+Q6NJi+I3QBl9/XIp74lVJXA81a63si3beILoWF5wDPduwjraGw8GyjQxJh0LNq5dLHHqM1\nRIVtjVyHhvL6G9aSDUqpWcADgANoBUqBUfgO7tbi+7rcqbW+zV8bUrJh+JHxy/jNPH4wvmRDuA/u\nbgKGztegEEKYgFy5K4QQJiOJXwghTEYSvxBCmIwkfiGEMBlJ/EIIYTKS+IUQwmQk8QshhMkYVbJB\niLBxuaq54463+eSTGKCCuXNTePjhiwdVObFnhckFhQ+RbssMXdBhEK7qkW6Xi9Xf/x7Pr7fhZCZT\n56by8MMXRKwy5VD8LKKNJH4x7KxYsZa33rqR4zVUVq9+jvj4tYOqnNizwuRKLCx5fGWIIg6PcFWP\nXLfiTtatho95GrBQtNpLfHzkKlMOxc8i2siuHjHsFBenQbd6kKkdzw1cenFRtxbTi4sG1V4k9JyH\nwc7BcenFRRQxMyxtB9r/UPssoo0kfjHsOBxufGWg6Pi/FoejZlBtuh2Obi26HQWDai8Ses7DYOfg\nOLfDQQHbwtJ2oP0Ptc8i2siuHjHsFBaeQ0vLk6xfHwNUMnduCoWFSwbVZs8KkwsKHwxNsGEUruqR\nCwofoq7ldkrXX4+TmUybm0Zh4fkhaTvQ/ofaZxFtwlqdMxSkOufwI+OX8Zt5/GB8dU7Z1SOEECYj\niV8IIUxGEr8QQpiMJH4hhDAZSfxCCGEykviFEMJkwn4ev1JqOvAK8KDW+o9KqTzgWXxfOmXAt7XW\nreGOQwghhE9Yt/iVUknAI8CaLk/fA/xea30WsB+4IZwxCCGE6C7cu3qagMX4tuyPWwi81vH4NeC8\nMMcgwszlqmb58lUsWvQey5e/TFVVdedr+/cXc9JJv8fhWMVJJz3CwYPFAbfrdrl4ffky1i1ayOvL\nr8Nd5QpH+AHpa4yhbnvTpi8CnrNA4+q5XPGBg1Ezt0aIpnXLCGHd1aO19gDNSqmuTyd32bVzDBgd\nzhhE+PVVBfLyy/+B0/ljwEJjo5elS++ntHR6QO1GUxXGcFW67K3tt9/+GU1N99B1zrZsuX1QcfVc\nzrnhNj51RsfcGiGa1i0jGF2rx+8lxcfZbElYrbGRiCWisrNTjQ4hZJxOG10rNTqdts7xVVfndXvN\n93Ng489ylnSrwpjlLDFs3voa40B0fW/Ptpubx9Fzzvz1FWhcPZerrs43dG6NXv+jYd0ycg6MSPy1\nSqkErXUzMAZw9rVwVVVDZKKKoOFWq8Rud+Grk+irf2+3V3WOLyOjhIaGL1/LyDgMEND4K+x5eNnQ\n8U6osOcbNm99jTFYPT//nm0nJBygqan7nPnrK9C4ei6XkVGCtwFD5jYa1n+j160I1erx+5oRiX8N\ncDnw147/3zIgBhFCfVWBXLXqUpYuvZ+qqjxstsOsWnVJwO1GUxXGcFW67K3tW2+9kGXLApuzQOPq\nudxdP/kuK+87EhVza4RoWreMENbqnEqpWcADgANoBUqBq4E/AwlAMXC91rrdXxtSnXP4kfHL+M08\nfjC+Ome4D+5uAnrbBFkUzn6FEEL4J1fuCiGEyUjiF0IIk5HEL4QQJiOJXwghTEYSvxBCmIzRV+4K\nIcSw4fV6cVY2sOdQFXtL3ewtceOub+6xlAWLBbxemGBP444rTiQhPrLVCSTxCyHEIHi9Xg4drWPD\n7mNs2H2U8uqmztdSRsQxNie1szzE8YuS4qyxuOua0SXV/P2D/Vy9aHJEY5bELwLiclVz5x1vs/OT\nGuxs48q5lSx++FHSbZlGhxYUt8vFuhV3dlyx6WBB4UODHoO/Nl2ualasWMuB/UkkutZya+a/8EwY\nz9KnngDiOl/3XU3rprDwHGy2jH7b67q81+srwLZ/fywuVzEjR04mb0wViy2vkVpSwh9d82gaeTbj\nx9d3tt8bf7EI/5wV9azfcYTPdn2Z7OOtMYypK6Kg5AuSkttY/MtfkJE58ivvzc5OxVlWzS+e3sB7\nmw5z8uQsphVE7ndJEr8IyIoVa3nzrRsBC0V4GbP6MlLi7xxyFQ3DUZXRX5tdK2LC5bzpvIznt6/i\n+ZvjOO/RJ/xW1gykvePLAx3PPQ/8GKfTwhdfeEnhdSyMZQOPQsdzfVUUDWf10eGkrrGVT3ce5ePt\nZRws8115mxAXy2lTRzF7yigOP/Bjbnzlxc4aQCtba/2uX3HWWG5aMo1fPbORp9/cxT03zmFEQmRS\nsiR+EZDi4jS6VncsYibpxW8aGdKApBcXdavKmF5cFLY2e87ZQWZi4R+kHDzY6+u+nwNv7/jyvudS\n6Pn5fPlaz+W/yn/boq3dw9Z9lXy8vYxt+ytp93ixWGDG+JHMn5HLSROziI/z7aNvKNof1Po1bnQa\nF8118NrHRTz/3l6uv3BqWMdynCR+ERCHw92xJejblilgG25HgcFRBc/tcODdsqlziywUY/DXZs85\nG8c2vEDduHG9vu5w1ATVnm95b8dztXStvlnANizABr7afm/8xWJWXq+XA2U1fLz9CJ/tPEp9UxsA\n+aNSmDc9lznTcshISfjK+wayfl08v4Ct+ypYt62MWZOzOXFiVmgH0wtJ/CIghYXn0NryJDvX+/bx\nL5jrHZIVDcNRldFfm8crYh44kERi5VouzCxm5YSvs/Sxx2ht919Zs7/2vlqJ89mOffz3M3LkZPLz\nqlmAl7SSQxx03daxj7+hz4qi4aw+OlR4vV6cFfVs3FPO+h1HOerylYRPT47n/NPymXtCLmNz+q6h\nP5D1yxobw01LpnH3yg2sXL2be2+aQ8qIuJCMyZ+wVucMBanOOfzI+GX80TJ+j8fLvlI3m/eWs3lP\nBceqGwGIs8Ywa3I286bnMq3ARmxMaC956m0O3lhfxEsfHmDOtBy+e8kJoejDmOqcQggRTWobWjh0\ntI5DR2spPlrLzqIq6hp9d4JNiI/l1CmjOHlSFidOyCIpMbLp8YI5Y9myt4JPdx7llMnZnDplVNj6\nksQvhBh2vF4vle4mio/WUXKslkNH6yg+WktVbfeLqdKT41l4kp2TJmUz1WEjzmpcMYPYmBhuXDKN\nXzz1Gc+8rZmUn0F6cnxY+pLEL4QY0to9Ho5UNnQm90NHayk5Vtd5QPa4jJR4Zk4YydicVMaOSmFs\nbirZ6YlYLP3e+jticjOTuPysCTz33l6eeWs3t319Rljik8QvhIiYdo+HppY2GppaaWv30tbuoc3j\npb3d0/lzS2s7LW0eWlo9tLS109LaTlNLO/VNrdQ3tVHf2EpDUxu1ja2465qpqW/F0+NYZY5tBNMK\nMhmbk4IjJ5WxOamkhWnrOdTOPTWPTXvK2by3gvU7jjBv+uiQ9yGJXwjRL4/HS01DC+66FmobW6it\nb6W2sZWGplYamttobGrz/d/cRnNHwm5uaae1zeP71+6hrd1DKM8libfGkJYcz/gxaeRmJuHISSV/\nVAr5o1IidiFUOMRYLNxw0VR+9tRn/OXdvUwZayMzLTGkfQzd2RFChIzH68XlbuJodSPl1Y1UVDdR\nXt2Iq6YJV20z7rqWr2xV+2ONjSEhLob4uFgS42NJTYojzhqDNdb3b8SIODxtHqyxFqyxMcR2/G+N\n8T2Oj4v1vd8aS/zxduJiSR4RR1KileTEOJITrZ0XTQ1H2RkjuPKcifz5Lc3Tq3dz5zdPDOkun4gn\nfqVUMvAMYAPigXu01u9EOg4hzMjj8XK0qoHS8noOl9fhrGzgSGU9R6saaW3zfGX52BgLGSnxjLen\nkZGaQHpyPGlJcaQmx5M6Io6kxDiSEqwkJVoZkWBlREJsv6c+RtPpnNHszBPtbNxTzvYDLj7c4mTh\nyWNC1rYRW/zLgN1a6/9WSo0G3gcic52yECbS2tZOybF6io/UdJ7dcri8/isJPiEuFvvIZHJHJpFj\nG0F2xpf/0lPiiYmig59mYrFYuH7xVH76xKe88P4+po3LZFTGiJC0bUTirwBmdDzOBMoNiEEEKVLV\nG8NRPXM4OD4vlv0Heq242e7xUFpez8GyGg6W1VJUVkNpRT3tni93z8TGWBiTlUzeqBTyslMYk53M\nmKxkbKkJYTlzpOc6c9ePT+aL++8my1lChT2v22c7kPVrMOvKQN8bzPtC8TtjS03g6q9N5vHXd/LU\nG7v44b+d3O8X8fEYva++/NkNXu9pvS3Tb+JXSl2gtX4rqGj7oLV+QSm1TCm1F8gALgpV2yJ8Blu9\nsa3dQ3VdM1W1zew+XENRaTVVtc00trRhwbfVmZgQy8HXVjG/uInk5jTGb9rPm//9Cy7/3UPDen9u\nII5X7LySS9jA7xlR10RFs4v/+p+PmTB9FIeO1NLSZUs+zhpDwehUCnLSGJvrO7PFnpWMNTZy56n3\nXGecG27jU2dH1VE2dKuMOpD1azCVVgf63mDeF6qKp6efkMPGPeVs2lPOms8Ps2h2fqAxzva3TCBb\n/N9TSj0K/AV4SmtdHFzY3SmlrgaKtdaLlVIzgSfpI0CbLQmrdfj90mdn913zI9o4nTa6Vm90Om19\njqGppY0v9lWwcfcxtuwpx1lRF9gZHSNnseOcWd2eeveBD8lKT8SencLorGTyc1Jx5KbiyPXtd46m\n87ADFcznX1PfQmNDHC/M+SZVuRfytdx3SEj2XYjkJYEDpW7G5qahHDYm5WcweayNsTmpxEYwyfem\n5zpTXZ3frXJllrOkcx6CXb/oeL+/9voz0PcG877+xhTMOvD9fzuFW3/zPi9/uJ8zT8knv4+aQV1j\n9KffxK+1vlApZQOWAo8ppQCeBl7WWrcHHPmX5gNvd7S9TSllV0pZtNa9poWqqoYBdBHdhuLBLbvd\nRdfqj3Z71VfG0NbuYdOecv75RRm7i6tpa/dtgSbGxzJpTDqZ6YnYUhMYOzqdOIvvz9ikjtPumjvO\n1f7oN//Dgo2bqE9MoTI5k4+nn0rSjNM4WtXAtn0VbNtX0a3PlBFxjMlKZkx2MnnZX+7CiObT+fr6\n/GvqWzjUcaVp0RHfLpsKdxPMuAqABKCxJgbnHjvVZRlMn/gZf3jgoq/cus/lqg/3MPrVc53JyCjB\n23D8J6iw53fOQyDrV08V9jy8bOi1vf4M9L3BvK+vMQ0kB3x70WT+sGo7v3n2c37y7Vl+D6J3jdGf\ngH47tNZVSqnngRbgFuC/gJ8rpW7SWn8SVPSwDzgdWKWUcgC1/pK+iB59VW9saW3nnQ0lrPm8hJoG\nX92TsaNSmDFhJNPHZTJhTHq3XQx9rfQ5d/+AdT/07UNtzknljh9c0bkPtaW1nWNVjTgr6yktr6e0\nop7S8jr2lFSjS6q7tZOVnkhedgr2LN9+7NyRSeRmJkXNF0JTSxtHXA04K+o7zrCpp+RYLdV1Ld2W\nS060Mn1cJqNtcdS98Rwj9UaePHIirSPPZsb4XRTec3bE79caqJ7rzF0/+S4r7zvSsY8/v1vlyoFU\nBx1MpdWBvjeY94W64ukpahSnn5DDJzuOsvqTQyyZV9BnjN5XX95wg5+2+q3OqZQ6E7geOBt4GXhc\na71LKVUArNJanxxM8B2ncz4F5ACxwF1a6w/9LS/VOaOX1+tl055yXnh/HxXuJpITrcyfMZqzTrIz\nemSy3/eFevzNre2UVdZz+JjvFEXfv3pq6lu+smxGSjy5mUlkZ4xgZHoimamJZKYlkJmWiC0lIWRJ\ntLm1HXddM9V1LVTWNFHp9p0Xf7SqkQp3I66anjfg9v0FNHZUCmM7LkQqyE1lZJSVFAiF4bL+D8ZA\n56C+qZWfPvEptQ2t/GzZbPJHpfTVh98VJ5DE/0/gT8CLWuvmHq/9WGt9f3ChB0cSf3Q6fKyO597b\ny67iKmJjLHzt1Hwunl8Q0BZ1pMZfU99CaUU9ZZX1lFU2cMTlO2fdVdOMv5UqzhpDcqKVlBFxJCfG\nMSLBSnzclxcfWSzQ3u6lzePx/d/uod3jpaW1ncZmX1mB2sZWmlt63wtqAbJtI8hKTyQ3M4kxWcm+\nv0qyU8Jegz1aDIf1f7AGMwfb9lfy8ItbyR+Vwk+vO9XvAftBJX6jSeKPLg1Nrbz80QHWbi7F64WZ\nE0Zy5bmTyM1MCrgNo8ff2tZOeXVT51WprpomXDXNVNc1U9/USl1jK/WNvhIEwUiIiyUp0UpqUhyp\nSfGkJ8djS/X9NZGZmtB5brx9dPqQ/fxDwejPPxoMdg5Wrt7FR1vLWDKvgK+fOd5fH1KPXwzeAWcN\nj72yncqaJnIyk7jq3EnMnDDS6LCCFmeNxd6xpd2Xdo+H5pZ2mls9tLf76s0AxMbGYI2xEBsbQ2yM\nBWushThrTMhv1iGEP986ZxI7Dlbx5vpiTpqYxXh7cPdIlsQvArKv1M0Dz2+hpbWdi+cVcPH8goie\nE26E2JgYkhJjSAptfSwhBm1EgpUbLprKb57bzJNv7OTny2YHda3L8P7NFSFRfKSWh/62ldY2D7cs\nnc7SM8cP+6QvRLSb6rBx3il5lFU28PJHB4J6r/z2ij6VVtTzwAtbaGpu48YlUzlFhe92cEKI4Fy+\ncAI5thG8u6GEPT1Oae6LJH7h17GqBn77/GbqGlu59gLF3BNyjQ5JCNFFQlwsN140DSzw5Bs7aWoJ\n7IQESfyiV66aJn7z3BbcdS1cee4kzjopdCVhhRChMzEvnQvmjKW8uokX1+4P6D2S+MVXuOua+c1z\nm6msaWLpmeP7LQoVzVyuapYvX8WiRe+xfPnLVFUF/ufwQLhdLl687ipWTh7La5PH8vJ1V+GucoW8\nHbfLxevLl7Fu0UJeX37dgPror/9g2g93PKJvl50xnjFZyazdXMr2g5X9Li9n9Yhu6hpb+e0LWzha\n1ciFpztYMtdhdEiDEqoKiYFat+JOMla/wc10VGhZ/QYr4xMCrhoZaDuDqUwZaP/BtB/ueETf4qwx\n3LRkGr985nOefnM3v1o+p8/lZYtfdGpsbuPBF7ZQWl7PubPyuPys8UO+XEBxcRpdKyT6fg6f9OIi\nUrv16Hsu1O2kFxcNuo/++g+m/XDHI/rnyE1l0ex8qmqb2by3os9lJfELAJpb2nn4xa0UHanljBmj\nueprk4Z80gdwONzQWaDBi8NRE9b+3A4HNd16BLejIOTtuB2OQffRX//BtB/ueERg5s0YDcAm3ff9\nrWRXj6C1rZ1HX97G3sNuTps6imWLpwyb2+2FukJifxYUPsTbLS38dv2/GAm0zT2Dc4KoGhloO4Op\nTBlo/8G0H+54RGDsHVVov+hnP7/U6jFANNUqaWv38MdV29myr4KTJmZxy9LpYb84K5rGbwQZv7nH\nD+Gdg5c+3M8b64t57YFL/W69ya4eE/N4vDzx+k627KtgWoGNmy87Qa7IFWKImzU5u99l5LfcpLxe\nL8+8vZvPdh1jYl46t399JnHD8BaXQphNQW4qmWkJfS4jid+kPtzi5KOtZThyUrnjGydG7V2chBDB\nsVgsnH1y3xdcSuI3oZJjdfx1zV6SE63cfvkMkhLlGL8Qw8lFcwv6fF0Sv8k0t7bzp1e309bu4caL\nppGZJjWHhTAbSfwm88L7+yirbOC8U/M4aVKW0eEIIQxgyN/4SqmrgR8ArcDPtNarjYjDbLYfqOSD\nzaXkZSdzxcIJRocjhDBIxLf4lVKZwM+AecAS4NJIx2BGzS3tPPO2JsZi4aYl0+QMHiFMzIgt/vOA\nd7XWDUAD8O8GxDDkuVzVrFixtuOKVDeFhedgs2X4XX7VugNUuJs496RRbLvrdoqLi3A7HCwofIh0\nW2YEI/cAaU5MAAATHUlEQVRxu1ysW3Fnx5WegcXR35i7vj46t5wFra/w+sZsnMxk6txUfv7zudx/\n/6aA5yxc48zOTg2qDZermjvveJudn9RgZxtXzq1k8cOPdpuvvuam57wstrxGbkkxe1yVjM3MpG3C\nxJCvB8Gun/4MZD0ZjFDF3Z/qykpeX/6dkI6r51zd9eoqW6nXW9XbskYk/gIgWSn1KpAB3K21ft+A\nOIa0YKpOFh2p4d3PSxiVMYKUlx+KiiqKA6nm2N+Yu72Oly/YzGH+AFgoWu1l69b7cTp/7Pf94dDb\nOCe+8lJQbaxYsZY337oRsFCElzGrLyMl/s5u89XX3PSclxRe529swws87yzlmu1fhHw9CFVV1EhX\n/YxUNdfVt9wS8nH1nKs98BhwZW/LGpH4LUAmcBkwDlgL+K39a7MlYR2GuyWC3erryem00bV2o9Np\n67XN9nYPv3x2I14vfO9bJ3Pote5VFLOcJYOOZSCynCVBx9HfmHu+7mJGt5+rq/P6fH849DZOCO7z\n7zmuImaS5Xy7z7F3HVtv74d/YAFSCM960N9nFWhfA1lPBiPQ36vBSjl4MOTj6jlXk2C8v2WNSPxH\ngY+11l7ggFKqVimVpbXutY5oVVVDZKOLgFDU6bDbXfjqIFoAL3Z7Va9trv60mAOlbubPyMVuS2ST\nPQ8vGzreBRX2/IjXTcnOTqViAHH0N+aer2fyBQ1dfs7IKKGhof85C6XexgkE1W/PcRWw7Svz1dfc\n9PZ+Op6pIzzrQV/xBLP+D2Q9GYxAf68Gq3bcOLwbQjuunnO1F/zegd2IxP8O8LRSqhDfln+yv6Qv\n/Auk6uSx6kZeXXeQ1KQ4vnXOJCB6qigOJI7+xtz1dfvoCs5oKeL1jdfjZCbT5qbx819cyn33Ra5S\nJ4RmvgsLz6G15Ul2rvft418w1/uVdvqam57zsgAvLx6eyZ5K3z7+lRMmhXw9CFVV1Eivr5Gq5nrh\nY4+xsrktpOPqOVfPvPryzff7WdaQ6pxKqeXATfi+mO7VWr/hb1mpzjkwXq+XB1/Ywo6iKr5z8TRO\nj6IbpZu9OqOM39zjh8jMQXZ2qt/qnIacx6+1fhx43Ii+zeKTHUfZUVTF9PGZzJmWY3Q4QogoIlfu\nDkO1DS08995e4uNiuHaRGhZ30hJChI4k/mHo+ff2UdfYytIF48nKGGF0OEKIKCOJf5jZcdDF+h1H\ncOSmct6peUaHI4SIQpL4h5Hm1naeeXs3MRYLyy6YQmyMfLxCiK+SzDCMvPrPg5RXN7HotHwcuZG/\nKEsIMTRI4h8mio/U8s5nJWRnJHLpGeOMDkcIEcUk8Q8D7R4PK9/ajcfr5drzp5AQN/xKXAghQifq\n77m3aNF7EaukOFSt+fwwxUdqmXtCLieM81/hr2vlwdzcciyWNsrKRg96frtWBSwdPYa3uARnWVbU\nfm7hqPgYTJvBVmaMdIVKMfxFfeKPGTOGj7eezg9/+E7YKykORRXVjaxad4CUEXFcee7EPpftWnnQ\nd9H0c8Blg65C2LUq4De35PEmvkqSkaqAGaxwVHwMps1gKzNGukKlGP6iPvGPnuTErpzUVaRRVlnP\n6JHJRocUNbxeL//37h5aWj1ce74iNSm+z+WLi9PoWnkQUjsf+14bmPTiLyt++io/ftnHYNoNl67x\nWjp+jmSbPSsz9td/OOIV5hb1+/g/+r+FVBzKYkSWhZ89+Rl/fXcPdY2tRocVFTbtKWfb/kqmOmzM\nDaAWj8PhxrelT8f/tZ2PHY6aAcfhdjg6W/VVfvyyj8G0Gy5d4/UCbkdBRNusHTcuqP7DEa8wt6jf\n4p+Qt5ZRrTVcd/7JvPFpGWs2Hmb9jiNctmA8C0+2m/Zc9fqmVv7vnT1YYy1cs2hyQGUZulYeHD26\nAmilrOyVQVch7FoVcMFoL/U82bGPPzIVMIMVjoqPwbQZbGXGaKmoKoYPQ6pzBqNrdc7WNg9rNpbw\n2r+KaGppx56VzJXnTmT6uJFGhhi0UFTme/rNXazbVsbXzxzPknkFoQksQsxenVHGb+7xg/HVOYfU\n5nKcNYbFcxzc/925nHminbKKeh58YSu/e3ErR1zD74Yt/uwscrFuWxljR6VwwZyxRocjhBhion5X\nT2/Sk+NZtngK58waw3Nr9rJ1fyXbD7o495Q8LplfQFJinNEhhk1zSzsrV/vKMlx/4VSssUPqu1sI\nEQWGdNYYm5PKD//tZG65bDq21ATe2VDCj/7fJ6zdXEq7x2N0eGHx8kcHqHA3ccGcsVKWQQgxIEM6\n8QNYLBZOnTKKXy2fw+Vnjae13cOzb2vufnoDO4tcRocXUvtL3az5vISczCQumV9gdDhCiCFqyCf+\n4+KssVw0t4D7v3M6Z8wcTWl5Pb99fgu/f2kbR4fBDdtb2zw8vXo3XuD6xVOIl7IMQogBGpL7+PuS\nkZLADRdO7dz/v3lvBdv2V/K12fksmVtAUuLQHPIb64twVtRz9qwxTM6PrhIIQoihxbAtfqVUolJq\nn1Lq2nC0X5Cbxo+unsW/X3oCGSkJvPXpIX7yv+v5cEspHk90n8LaU8mxOt5YX0xmWgLfOGuC0eEI\nIYY4I3f1/BSoDGcHFouF06bm8Kvlc1h65niaWz38+S3N3Ss3sLu4Kpxdh0y7x8PTb+6i3ePl2vMV\nIxKG5l8sQojoYUgWUUopYArwRiT6i4+L5eJ5BZwxYzQvfbifj7cfofC5zZwyOZsrzpnIqCi+L+27\nGw5TdKSWuSfkMHNCVrfXulbbjNZKmENBb9Uv270xAc2tWT4DqRA6vBi1+fgAcCuwLJKd2lITuGnJ\nNM49JY/n1uxl455ytu6vYNHssVw01xF1W9NHqxpYte4AqUlxXHXe5K+83rXaZrRWwhwKeqt++SqX\nBjS3ZvkMpELo8BLxTKeU+jbwsda62LfhT59FZmy2JKzW0J7Bkp2dyuwZdj7aXMrKN3by5ifFfLzj\nCNcunsq5s8cSE9N/3ZtQxNAXj8fLgy9upbXNw/evmsW4sV/dunI6bXSthOl02vptN1pEU5xZzpJu\n1S+znCU4CWxuB/oZRNP4A9HbHA1mDENt/OFg5BwYsYl7ETBOKXUxkAc0KaVKtNbv97ZwVRhPxZyW\nn869N57G258e4s1Pi3nkb1t49cP9XHXepLCeORNInY4PtpSyfX8lJ0/KQtl7X95ud+Gr1+irr2+3\nVw2JGijRVqulwp6Hlw2ddymosOdjJ7C5HchnEG3jD0RvczTQMQzF8YdahGr1+H0t4olfa33l8cdK\nqZ8DB/0l/UhIiIvlkjPGccbM0fz9w/18suMov/7LJk6dMopvLpxAlgH7/101Tby4dh8jEqxcs0j5\nrbzZtdpmtFbCHAp6q345nxgCmVuzfAZSIXR4MbQ6Z5fE/4y/ZbpW54yE/aVunntvLwecNVhjY7hg\nTj4Xnu4gMT5035F9fdt7vV4e+fs2tu6vZNniKZx5oj1k/UYLs2/xyfjNPX4wvjqnoUcztdZ3G9l/\nbyaMSecn3z6FT3ce5e8f7Of1j4tZt62Mb5w1gbnTc4kJoO79YHy26xhb91cyZWwGC2aODmtfQghz\nGjYlG0IpxmJh7gm53Lf8dC6ZX0BDUxtPvrGLX/75c/Yerg5bv7UNLfzl3T3EW2NYtnhKQDdXEUKI\nYEni70NCfCyXLRjPfctPZ860HIqO1HL//23iT69up9LdFPL+nntvL3WNrXz9zPGMsiWFvH0hhIBh\nWKsnHEamJ/LdS07g3Fl5/HXNHj7bdYwteyu4YM5YFs9xkBA/+NNNt+6r4JMdRxk3Oo3zTs0PQdRC\nCNE72eIPwsS8dO667lRuvGgqIxKt/ONfRfzk8U9Yv/0InkEcJG9sbuOZtzWxMRauv3BKRK4jEEKY\nlyT+IMVYLMyfMZr7v3M6S+Y5qG1o5fHXd3LfsxvZ73QPqM0XP9hPVW0zS+YVkJedEuKIhRCiO0n8\nA5QYb+XrZ07gvuVzmD1lFAecNfzqmY08/toOXDWB7//Xh6r4YHMpY7KTuWiuI4wRCyGEj+zjH6Ss\njBHcfNl0zi2p5rk1e1m/4ygb95Rz4RwH588ZS0IfN0xpaW3n6dW7sVjg+sVy/1whRGRI4g+RyfkZ\n/HTZqfxrWxkvfXSAV/55kI+2OfnGwgnMmZrTeWqmy1XNbbe9zp49I8g/uRmvLZFFs/MZb08zeATD\nj1kqZ0YjqeYZ3STxh1CMxcKCE+2cOmUUb6wv5p0Nh/jff+zk/Y2lXHXeJMaNTuus5pieU82YhR9h\naW1n6ZnjjQ59WDJL5cxoJNU8o5sk/jAYkWDlGwsncOZJdl5cu4+Nupx7//w586bnUlKWhiXGy4mL\ntmCJgfItMX3uDhIDV1ycRtfKmb6fRSSkFxd1q+aZXlxkYDSiJ9mpHEajMkZw69IZ/PCqk8kflcLH\n249gnwezL/uEtOwairc5yE2vMTrMYcvhcOOrJQngxeGQuY4Ut8PRZebB7SgwMBrRk2zxR8AUh42f\nL5vNum1OXvpgPxSU097sZVzqNgp/PTyrOUYDs1TOjEZSzTO6GVqdMxCRrs4Zbg1NbWzaX0n+yCQc\nuea8GYXZqzPK+M09fjB5dU4zSkq0snThRNOv+EII48g+fiGEMBlJ/EIIYTKS+IUQwmQk8QshhMkY\ncnBXKVUInAHEAr/WWq8yIg4hhDCjiG/xK6UWAtO01vOAxcDDkY5BCCHMzIhdPR8CV3Q8rgaSlFJy\n5xEhhIiQiO/q0Vp7gcaOH28C3ux4ToiQkcqcQvhn2AVcSqlLgeuBRUbFIIYvqcwphH9GHdw9H/gx\ncL7Wus9LWG22JKzW4Ve9MjvbnOUajgv3+J1OG10rczqdtqia82iKxQhmHz8YOwcRT/xKqTSgEDhX\na93vTWqrqhrCH1SEmb1WSSTGb7e78NWFtABe7PaqqJlz+fzNPX6IWK0ev68ZscX/LWAk8LeOg7pe\n4Fqt9WEDYhHDlFTmFMI/Iw7uPg48Hul+hbnYbBmyT18IP+TKXSGEMBlJ/EIIYTKS+IUQwmQk8Qsh\nhMlI4hdCCJORxC+EECYjiV8IIUxGEr8QQpiMJH4hhDAZSfxCCGEykviFEMJkJPELIYTJSOIXQgiT\nkcQvhBAmI4lfCCFMRhK/EEKYjCR+IYQwGUn8QghhMpL4hRDCZIy42TpKqQeB0wEPcIfW+nMj4hBC\nCDOK+Ba/UupMYKLWeh5wE/BIpGMQQggzM2JXz7nAKwBa691AhlIqxYA4hBDClIxI/LlAeZefKzqe\nE0IIEQHRcHDXYnQAQghhJkYc3HXSfQvfDpT5Wzg7O3VYfjFkZ6caHYKhZPwyfrMzcg6M2OJ/B/gG\ngFJqFlCqta43IA4hhDAli9frjXinSqn7gLOAduBWrfUXEQ9CCCFMypDEL4QQwjjRcHBXCCFEBEni\nF0IIk5HEL4QQJmNIrR4zUkolA88ANiAeuEdr/Y6xUUWOUsoC/AmYDjQD/6613mNsVJGhlJqO72r1\nB7XWf1RK5QHP4tvwKgO+rbVuNTLGcOo5/o7nvgf8FsjQWjcYGV8k9LIO5ANPAXFAC3CN1vpYpOKR\nLf7IWQbs1lqfA1wB/M7YcCLuUiBNaz0fX42mBwyOJyKUUkn46lGt6fL0PcDvtdZnAfuBG4yILRJ6\nG79S6tvAKKDUqLgiyc86cC/wJ631QnxfCP8ZyZgk8UdOBTCy43Em3ctWmMEk4DMArfUBwNHxV8Bw\n1wQspvtFiguB1zoevwacF+GYIqm38b+stb7LoHiM0Nsc3Ay83PG4HF9OiBhJ/BGitX4BX7LbC3wA\n/JexEUXcF8D5SqkYpZQCxgFZBscUdlprj9a6ucfTyV127RwDRkc4rIjpbfxmu2DTzxw0aq29SqkY\n4Fbgr5GMSRJ/hCilrgaKtdaT8FUo/YPBIUWU1votfFv8HwLfA3YhdZpA5sC0OpL+s8B7Wuu1kexb\nEn/kzAfeBtBabwPsJtnV0Ulr/TOt9QKt9a1AZiQPZkWZWqVUQsfjMfjqV5mR2a8efRrQWut7I92x\nJP7I2YfvrmMopRxArdbaNCu+UmqmUurJjscXABsNDslIa4DLOx5fDrxlYCyR1HNDx1QbPl117AFo\n1lrfY0T/UrIhQjpO53wKyAFigbu01h8aG1XkdPx18yRwAtAIXK21HvZndXQUInwAcACt+M5kuRr4\nM5AAFAPXa63bDQsyjPyM/11gETAH2ACs11r/yLAgw8zPHIzCd9C3Ft9fPju11rdFKiZJ/EIIYTKy\nq0cIIUxGEr8QQpiMJH4hhDAZSfxCCGEykviFEMJkJPELIYTJSOIXQgiTkcQvhBAmI4lfiCAppb6v\nlPrfjsdKKbWr48psIYYESfxCBO9hYLJSah6+KqvLzVZqWAxtkviFCFJHcb0bgb8B27TW/zQ4JCGC\nIolfiIEZia/A1lijAxEiWJL4hQiSUioReAy4GGhRSl1jcEhCBEUSvxDBuxvffWP3AXcAv1BK2Q2O\nSYiASVlmIYQwGdniF0IIk5HEL4QQJiOJXwghTEYSvxBCmIwkfiGEMBlJ/EIIYTKS+IUQwmQk8Qsh\nhMn8f9j4/Vqo5jP7AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# λの予測値と一緒にプロット\n", "ax = d.sort_values(by=['x']).plot(x='x', y='lam_7')\n", "d_T.plot(kind='scatter', x='x', y='y', c='r', ax=ax, label='T')\n", "d_C.plot(kind='scatter', x='x', y='y', c='b', ax=ax, label='C')\n", "plt.legend(loc='upper left')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

さまざまな逸脱度

\n", "\t

\n", "\t\t逸脱度は、「あてはまりの悪さ」を表現する指標と久保本では説明されています。\n", "\t

\n", "\t

\n", "\t\t最大対数尤度を$log L^*$とすると、逸脱度(D)は以下の様に定義されます。\n", "$$\n", "\tD = -2 log L^*\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\t4章で登場する逸脱度を以下の表に示します。\n", "\t\t最大の逸脱度は、一定モデルで、最小の逸脱度はフルモデルです。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "

名前	定義
逸脱度(D)	-2 log L*
最小の逸脱度	フルモデルをあてはめたときのD
残差逸脱度	D-最小のD
最大の逸脱度	NullモデルをあてはめたときのD
Null逸脱度	最大のD - 最小のD

" ], "text/plain": [ "<__main__.Table2Html instance at 0x7ff6581e4050>" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 様々な逸脱度\n", "hdr = ['名前', '定義']\n", "tbl = [ ['逸脱度(D)', '-2 log L*'], \n", " ['最小の逸脱度', 'フルモデルをあてはめたときのD'],\n", " ['残差逸脱度', 'D-最小のD'],\n", " ['最大の逸脱度', 'NullモデルをあてはめたときのD'],\n", " ['Null逸脱度', '最大のD - 最小のD']]\n", "Table2Html(tbl, hdr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t各モデルのglmのfit結果を求めて、その結果を表にまとめます。\n", "\t

\t\n", "\t

\n", "\t\tここで出てくるモデル選択の基準となるAIC(Akaike's information criterion)は、以下の様に定義されています。\n", "$$\n", "\tAIC = -2 { (最大対数尤度) - (最尤推定したパラメータの数)} = D + 2k\n", "$$\t\t\n", "\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 逸脱度（Deviance）は、-2 log L* と定義されており、glmの結果の-2*llfから計算できる\n", "# 各モデルに対するfitを計算してDevianceの違いを整理してみる\n", "fit_f = smf.glm('y ~ f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "fit_x = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "fit_x_f = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\tデータの個数と同じ数のパラメータを使ってあてはめたモデルがフルモデルです。\n", "\t\tこれは、実際にそのようなモデルを作って計算するのでは無く、フルモデルはすべての観測値が予測値と一致するモデルであり、\n", "\t\tその結果、残差は０であり、対数尤度は、以下のようになります。\n", "\t

\t\n", "" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "-192.889752524496" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# フルモデルの対数尤度の計算\n", "fullLogL = sageobj(sum([r.dpois(y, y, log=True) for y in d.y])); fullLogL" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t結果をみやすくするために、小数点以下3桁を表示する関数prf3を定義します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 結果を小数点3桁で表示する\n", "def prf3(f):\n", " return '%.3f'%f" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "

モデル	k	log L*	deviance -2 log L*	residual deviance	AIC
一定	1	-231.863	463.727	77.947	479.727
f	2	-237.627	475.255	89.475	479.255
x	2	-235.386	470.773	84.993	474.773
x+f	3	-235.294	470.587	84.808	476.587
フル	100	-192.890	385.780	0	585.780

" ], "text/plain": [ "<__main__.Table2Html instance at 0x7ff65b94b488>" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# フルモデルは、k=Nで結果がすべて観測値と一致します。従って残差は0となります。 \n", "# log L*: fit.llf, residual deviance: fit.deviance, AIC: fit.aicで取得できる\n", "hdr = ['モデル', 'k', \"log L*\", 'deviance -2 log L*', 'residual deviance', 'AIC']\n", "tbl = [['一定', '1', prf3(fit_1.llf), prf3(-2*fit_1.llf), prf3(fit_1.deviance), prf3(fit_1.aic)],\n", " ['f', '2', prf3(fit_f.llf), prf3(-2*fit_f.llf), prf3(fit_f.deviance), prf3(fit_f.aic)],\n", " ['x', '2', prf3(fit_x.llf), prf3(-2*fit_x.llf), prf3(fit_x.deviance), prf3(fit_x.aic)],\n", " ['x+f', '3', prf3(fit_x_f.llf), prf3(-2*fit_x_f.llf), prf3(fit_x_f.deviance), prf3(fit_x_f.aic)],\n", " ['フル', '100', prf3(fullLogL), prf3(-2*fullLogL), '0', prf3(-2*fullLogL+2*100)]]\n", "Table2Html(tbl, hdr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

平均対数尤度を使ったモデルの予測の良さ

\n", "\t

\n", "\t\tいよいよ図４．９の計算に入ります。久保本ではモデル選択の指標としてAICが使えるのかを一定モデルを使って\n", "\t\t数値例で示しています。（これが久保本のすごいところ！）\n", "\t

\n", "\t

\n", "\t\t平均の対数尤度は、真のモデルを使って生成されたサンプルから求めた対数尤度の平均値であり、実際のモデルづくりでは\n", "\t\tこのようなことは不可能ですが、平均対数尤度と推定されたモデルの関係とこの差（バイアス）bの分布からAICの意味するところを\n", "\t\t表現しています。bの定義は以下の通りです。\n", "$$\n", "\tb = log L^* - E(Log L)\n", "$$\t\t\n", "\t

\t\n", "\n", "\n", "\t

観測データが一つの場合

\n", "\t

\n", "\t\t図４．９の(A)観測データが一つの場合について、真のモデルλ=8のポアソン分布から50個のサンプルを生成します。\n", "\t

\n", "\t

\n", "\t\t準備としてサンプル生成関数genSample, 対数尤度を計算するlogLを定義します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# サンプル生成λ=8ポアソン分布のサンプルをN=50個生成する\n", "# sample = [ float(y.n()) for y in r.rpois(N, lambda_true)]\n", "# sageobjを使って変換する方が速いのでgenSample関数を使用する\n", "def genSample(lambda_target, sample_size):\n", " return sageobj(r('rpois(%d, %f)'%(sample_size, lambda_target)))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Rの乱数の種をセット\n", "r('set.seed(101)')\n", "# 真のλ=8から50個のポアソン分布のサンプルを生成する\n", "N = 50\n", "lambda_true = 8\n", "sample = genSample(lambda_true, N)\n", "sampleDf = pd.DataFrame({'y' : np.array(sample)})" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 対数尤度の計算\n", "def logL(sample, lambda_estimated):\n", " return sageobj(sum([r.dpois(y, lambda_estimated, log=True) for y in sample]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\tglmを使って一定モデルの結果を取得します。ここでは対数尤度が-116.7、 $\\beta_1$の値は2.06と求まりました。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# k=1でのglm回帰を実行\n", "fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-116.764129142 2.06686275947\n" ] } ], "source": [ "# fit.summary2()\n", "# 最大対数尤度\n", "mxLogL = fit.llf\n", "# 切片\n", "beta_1 = fit.params[0]\n", "print mxLogL, beta_1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\tこの結果を図化し、logL_pltにセットします。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# β= [1.95, 2.20]の範囲の対数尤度を求める\n", "logL_lst = [ (x,logL(sample, exp(x)) )for x in np.arange(1.95, 2.20, 0.01)]\n", "# log L*の曲線プロット\n", "logL_plt = list_plot(logL_lst, figsize=4, plotjoined =True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\tつぎに、真のモデルから50個のサンプルを200セット生成し、平均の対数尤度を計算します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 真のλ=8から生成されたサンプル200セットにβの推定値を使って対数尤度を計算\n", "# logL200 = [ logL(genSample(lambda_true, N), exp(beta_1)) for i in range(50)]\n", "# とても遅いので、\n", "logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t求まった結果を一つにプロットします。logL*を赤○で、平均対数尤度を緑○で、各対数尤度を小さな青○で示し、\n", "\t\t各βの対数尤度曲線と一緒に表示したのが、以下の図です。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAECCAYAAAARlssoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt8U/X9P/DXOUmaJm0T0tLSi2BLkVu5tNVCVUSFfYFt\ntfBDwH21Oh1U9vXr0Ll+5ffbZA+v4BwqOi/o6NCturmv+gWZDipf3ZwoolCY3G/FWuroLbRNkzZN\n8vn9EXpKSYtNctq0yev5ePSR9pyQvN+cJK+cz7lJQggBIiIiAHKoCyAiosGDoUBERAqGAhERKRgK\nRESkYCgQEZGCoUBERAqGAhERKRgKRESkYCgQEZGCoUBERIqwD4Vt27YhOTkZN99880Xv53K58Mtf\n/hKZmZmIi4vDd77zHVRWVirzr7vuOkRFRcFoNMJgMMBgMCAnJ6e/y+83VVVVWLhwIYYPH46UlBTc\ncccdaG5uDnVZAfGnl2effRbjx4/HsGHDMHPmTOzZs0eZN1SXsT/9t7a2oqioCLIs4+jRowNcaeDU\n6jESlvH69esxfvx4mEwm5Obm4p133vHrucI6FH7961/j3nvvxdixY7/1vmvWrMEf/vAHbN68GfX1\n9bj66qsxf/58Zb4kSSgtLYXdbofD4YDD4UBFRUV/lt+vbrjhBsTHx+Prr7/G7t27ceDAAZSUlIS6\nrID0tZctW7bgoYceQllZGc6cOYOCggIUFBTA4XAAGLrLuK/9f/PNN7j88suh0+kgSVIIKg2cWj2G\n+zJ+++238fOf/xyvvPIKrFYr7r77bixZsgSnTp3q83OFdSgYDAbs2rULmZmZ33rfLVu2oLi4GJMm\nTYJer8eDDz6Iuro6fPbZZ8p9wuXcgU1NTcjLy8OaNWtgMBiQmpqKH/7wh/joo49CXZrf/Onl5Zdf\nxh133IErrrgCer0e//Vf/wVJkrBlyxblPkNtGfvTf11dHX7961/jwQcfHFJ9qt3jUOod8K9/h8OB\nNWvWID8/HxqNBj/60Y8QFxeHnTt39vn5wjoU7r77bsTFxfX5/ud/s5AkCWazGXv37lWm/elPf0JW\nVhZMJhPmzJmDkydPqlrvQDGbzdiwYQMSExOVaVVVVUhLSwthVYHxp5fdu3cjNzdX+VuSJGRnZ+Pz\nzz9Xpg21ZexP/1OmTMENN9wwkOWpQu0ew3kZ33LLLVi+fLny99mzZ9HS0uLXezusQ8EfBQUFeOml\nl7B//344nU688MILqK6uRmNjIwBg4sSJmDx5Mnbs2IFTp05h+PDhmDdvHlwuV4grD94XX3yB5557\nDg888ECoSwnaxXppaGiAxWLpNi0+Ph719fUAwmMZh9Oy7E0wPWZlZUXUMi4uLsaVV16Ja665pu9P\nICLA7bffLv793//9ovdpa2sT9957r0hNTRXJycniF7/4hfje974nnnjiiR7v39LSInQ6nfjggw/6\no+QB8/HHH4v4+Hjx3HPPhbqUoH1bL1FRUeLdd9/tNq2oqEjcfvvtPd5/qC3jvi7LU6dOCUmSxJEj\nRwaoMvWo3WO4LuOOjg5x8803i6ysLFFbW+vXc2iDCKywotfr8fTTT+Ppp59Wpk2ZMqXX1a7Y2FjE\nx8ejpqZmoEpU3ZYtW3Drrbfi+eefxy233BLqcoLSl14SExPR0NDQbVpDQwMmT57c4/2H0jIOp2XZ\nm/7oMRyXcVtbGwoLC9HW1oZ//OMfPmvH3yqY1Boq+rKmsGfPnm7fFqqrq4VGoxGVlZWiublZ3HXX\nXeKbb75R5tfV1QlZlsUnn3zSb3X3px07doj4+Hixffv2UJcStL72Mn/+fHHPPfcof7vdbjFixAjx\n9ttvD+ll7O+yPHXqlJBleUitKajRY6Qs4/nz54sbbrhBOJ3OgJ4rokNh/PjxYseOHUIIIV555RWR\nnJwsjh8/LpqamsT8+fPFokWLlPvm5uaKRYsWicbGRtHY2CgWL14scnNzB6wHNblcLjFx4kTx29/+\nNtSlBO3bejl/GW/dulVYLBaxc+dOYbfbxUMPPSQuvfRS0dbWJoQYmsvYn/47VVZWDqnhIzV7DPdl\nXFZWJsaMGSMcDkfAzxfWoRAdHS0MBoPQarVCq9Uqf3eSZVls27ZN+bukpEQkJCSIYcOGiVtvvVU0\nNzcr877++mtx4403iuHDh4u4uDixcOFCUVNTM6D9qOUf//iHkGVZGAwG5f+k87aqqirU5fnlYr18\n9dVXPst4/fr1YtSoUcJgMIiZM2eKAwcOKPOG4jL2p/9HH31UREdHi+joaCHLsnK/xx57LMRdXJya\nPYbrMi4vLxdCCDF79myh0+mEwWDodr8777yzz88nCTHEdtolIqJ+w11SiYhIwVAgIiIFQ4GIiBQM\nBSIiUjAUiIhIwVAgIiIFQ4GIiBQMBSIiUjAUiIhIwVAgIiJFxJ06WwiBlpaWUJdBRDQg4uLi/Lom\nd8SFQktLC8xmc6jLICIaEE1NTTCZTH2+f8SdEI9rCkQUSfxdU4i4UCAiot5xQzMRESkYCkREpGAo\nEBGRgqFAREQKhgIRESkYCkREpGAoEBGRgqFAREQKhgIRESkYCkREpBjSoSCEQHNzM3imDiIidQzp\ns6R2nvHU37MA1tX1fEI8WZYQHx+DxsZWeDzhHzSR1i8QeT2z3/B3sZ4TE+P8fzy1CgsHsixBkiTI\nct/PKDiURVq/gLo9f/CBBvPnG3DTTQYcOjQ430qRtowjrV9A/Z6H9JoCUaicPi3h9tsNaGvzvhEP\nH5ZRUdEKeXBmA1Gf8SVMFIBTp2QlEADgm29kNDWFsCAilTAUiAIwaZIbqake5e8rrnDDYglhQUQq\n4fARUQDMZuDdd+149VUdjEZg2TJnqEsiUgVDgShAaWkCP/85w4DCC4ePiIhIwVAgIiIFQ4GIiBQM\nBSIiUgQdCi6XCyUlJdBoNCgvL/eZf+bMGcydOxeyLMPp7L5Rzmq14qabbkJycjLS0tJQXFyM9vb2\nYEsiIqIABRUKdrsdM2bMgNVq7XH+/v37MW3aNCQlJUGSfA/BXrZsGRwOBw4dOoTdu3fj0KFDWLly\nZTAlERFREIIKBZvNhqVLl6K0tLTHM5XW1tbijTfewLJly3qct3nzZqxZswYWiwXJyclYtWoVNm7c\nCLfbHUxZREQUoKBCISkpCcXFxb3OnzVrFvLz83uct3fvXmi1WmRlZSnTcnNz0dLSgsOHDwdTFhER\nBShkB681NDTAbDZ3mxYfHw8AqK+v79fnluWezyio0cjdbsNdpPULRF7P7Df8qd1zSI9oDtXFceLj\nY3rcxtHJZDIMYDWhF2n9ApHXM/sNf2r17FcolJWVKcNFkiTBbrcH/MSJiYloamqCEEL5gG5oaADg\nHZbqT42Nrb2uKZhMBjQ3O+B2e3r4l+El0voFIq9n9hv+LtazxRLj9+P5FQpFRUUoKiry+0l6kpOT\nAyEE9u3bh+zsbADArl27YLFYMG7cOFWeozcej7joVZncbg9crsh4QQGR1y+gXs8nTkjQ64FLLhnc\nV/mKtGUcaf0C6vU8IANvncNE5w8XJSQkYNGiRXjggQfQ0NCA6upqPPLIIyguLobMK5XQEHDXXdG4\n8spYXH55DJ59NirU5RCpIqhP37KyMhgMBhiNRkiShMLCQhiNRixfvhwAcOedd8JgMGDevHkAgGHD\nhsFoNOK1114DAKxfvx4mkwkZGRnIzs5Gfn4+Hn300SBbIup/e/bIePNNHQBACAmPPRaF1tYQF0Wk\nAkmEamuvCpqbm2E2m9HU1ASTydTnf1dX19LjdK1WhsUSA6u1NSJWPSOtX0C9nisqZMyd2zVeK8sC\nJ0/aYDSqUaV6Im0ZR1q/wMV7TkyM8/vxOE5DFICcHA9uucV72hZJEnjwwfZBFwhEgeBFdogC9PTT\n7fjpT52IjgaSkobsCjdRNwwFoiCMGsUwoPDC4SMiIlIwFIiISMFQIApCVZWEM2d6P2UK0VDDUCAK\n0E9/qscVV8RiypQYrF+vC3U5RKpgKBAFoKJCxmuveY9iFkLCgw/qEcSpwIgGDYYCUQAuPORTCN9p\nREMRQ4EoALm5Hixe3AHAe/DaL37hRIz/J6QkGnR4nAJRgJ5/vg0//Wk7oqMH/1lSifqKoUAUhLQ0\nAS3fRRRGOHxEFKA1a6KQnh6L0aNj8d//zWSg8MBQIArA/v0ynn5aDyEktLdLuPfeaLS3h7oqouAx\nFIgCcOZM9787OiR0dISmFiI1MRSIAhAV5XsUM0OBwgFDgSgAY8Z4YDR27XE0cqQHZnMICyJSCUOB\nKAApKQKvv+7ArFkufP/7Hfjzn+3gpcUpHHCXCaIAXXWVG1dd5Qh1GUSq4ncbIiJSMBSIAuRyAW++\nqcW2bZpQl0KkGg4fEQXA5QKys2NQW+v9XpWX58K773IoiYY+rikQBeCVV7RKIADA559rUVsbwoKI\nVMJQIApAXZ3vW8fBFQUKA0GHgsvlQklJCTQaDcrLy33mnzlzBnPnzoUsy3A6nd3mpaenQ6/Xw2g0\nwmAwwGg0YsGCBcGWRNTvLr/cfcEUgeHDQ1IKkaqCCgW73Y4ZM2bAarX2OH///v2YNm0akpKSIEm+\nR4BKkoTt27fDbrfD4XDAbrdj06ZNwZRENCAaGy98PUtoawtJKUSqCioUbDYbli5ditLSUogeLjtV\nW1uLN954A8uWLev1MXr6d0SDXUaGx2eaXh+CQohUFlQoJCUlobi4uNf5s2bNQn5+/kUfY926dRgz\nZgxMJhMWL16Murq6YEoiGhB5eR585zsu5e+7725HbGwICyJSSUh3Sc3NzcW0adNQVlYGq9WK2267\nDUuWLMGHH37Yr88ryxJk2Xc4S6ORu92Gu0jrF1CvZ5cLGDUKyMnx/j18uAytdvD9P0baMo60fgH1\new5pKLz11lvK70ajEc8//zwmTpyIyspKZGRk9NvzxsfH9LiNo5PJZOi35x6MIq1fQJ2eS0vP/0t3\n7mdwirRlHGn9Aur17FcolJWVKcNFkiTBbrerUkSn9PR0AEBNTU2/hkJjY2uvawomkwHNzQ643b5j\nxuEm0voF1Ou5okJGScn5b0KB996zD7rtCpG2jCOtX+DiPVssMX4/nl+hUFRUhKKiIr+fpCdVVVV4\n/PHH8cwzz0Cn837DOnjwICRJwujRo1V5jt54PAIeT+8buN1uD1yuyHhBAZHXLxB8z8eOyaioOH+K\nBKvVM2h3S420ZRxp/QLq9Twgw0edexidv6dRUlIS3nnnHWi1Wjz++OM4e/Ys7rvvPhQWFiIlJWUg\nyiIK2OnTvmuavBwnhYOgtkyUlZUpB51JkoTCwkIYjUYsX74cAHDnnXfCYDBg3rx5AIBhw4bBaDTi\ntddeQ3R0NLZt24YjR44gLS0NkyZNwpgxY/Dqq68G3xVRP8vMvPAbmYDFEpJSiFQliSF8oEBzczPM\nZjOamppgMpn6/O/q6lp6nK7VyrBYYmC1tkbEqmek9Quo1/MHH2jwgx8Yu007cqRl0AVDpC3jSOsX\nuHjPiYlxfj9e5Oy3RaSir77yHT5yuXq4I9EQw1AgCsCXX/q+dVpbQ1AIkcp4PQWiADQ0SLgMR3Ed\n/oYO6PBXfBdWawzO7VVNNGQxFIj8JNlasGpPEaZjC2R4N8k5oUP9hmXAb9YAMlfAaejiq5fIT3E/\nXoorz7yjBAIARKEDqf/9ImLWPBLCyoiCx1Ag8oPm4AHoy7f2Oj+69GVItp73biMaCjh8RGFHCKCt\nDbDZJNhs3tvWVgmtrYDDIUOrBRobtXA4BJxOoL1dQns74HQCTmfX7+3tEpxOwOPxPqYQgFwZDw02\nQUBSfjyQld9hAzw3ytCMiIZOB+h0gEYD6HQCWq33b++tUH7X64GYGAGjUcBoRK+3BoOAwQBc5LRd\nREFjKNCg5HIBTU0SrFbvBW3OnpXQ2CjBavX9vbW168O/Mwjc7m/75NQjKkogKgrQ67235/+u1wNR\nUUL5UJck74/GIwPQKJGggfu8ePAOJ7VpAafLe9Edlwvo6ABcLvncbec0SZnX1ibBbvcG0reRJIG4\nOGDYMAGzWVxw6zs9IUFCRkZnP8EvFwp/DAUaMC4XUF8v4cwZCbW1EmprZeV3762Mujrvh31TU88f\nkDExAhaL92fYMO/tyJEexMQAsbFCufX++P5uNktITo5Ba2trQCdM05xohuWqQki9HPPpsVjQ8OZh\nwOD/BZtdLsBuB+x2qYfbrt+bm72BefaspNyePi2f+xs4e1aCx3Ph/18MYmIEEhM7fzxISjr/765p\nI0Z410goMjEUSBXt7d7zAZ0+LeP0aQnV1d7bmpquD/6GBglCdH1YSZJAQoJAUpL3Z/RoD6ZP98Bi\nAeLjvR/68fHdQyDYs5BqtRL0eu+HbyDcmZeh/f/ciOi33+xxvv0/70Ggn6haLWAyASZTZ+AEdrIB\nIQCbzRsONpsG7e0GnDzZhn/9C6ir8wZvba2Ezz/3/t7QIPmsWcXHe5CWJpCW5r1NTe36PS3Ng+Rk\n73AYhR8uVuoTmw2orJRRVdX9Q//0aRnV1d5v+ecbPtyDkSMFUlK8H/Sd30CTkjznbgWGD/cOzww1\nLetewO5Dccg7VIYodAAAmhGH5uUroF9xX4ir8w5zxcUBcXECWq03ZK1Wd6+nffB4vEN0tbUS6uok\n/Otf3jCvrvbefvKJjOpqGS0tXcEhywLJyd6wGDnSg4wMD9LTPUhPF8jI8C5vbvsYmhgKpLBavR/8\np07JqKyUz/0uobJSRl1d14e+Xi+QliZwySUejB/vxuzZ3t87p6WkhPnwQ3Q0di59HgtL1mAGPkYH\ndPgAs1DxM4FBdjmFPpFlYPhwb0hfTHMzcPq0jJqarjXC06dlfP21hJ07dfjmm67XiNEokJ7eGRbe\noOj8SUkR0Gj6uysKFEMhwrS1ASdOyDh2TMaxYxpUVwNHjkTj5EnvmHSnhISuN/PMmR3KG3rUKO/4\nc6R/C9RogFqMwNu4UZnW0RHeu6J6h7Y8mDABANw+8+124Kuvun+ZqKyUsWWLBtXVXds59HqBMWM8\nGDfOg7FjvT/jxnlfXxySCj0ugjBlswFHj8o4etQbAEePanD0qIyvvup6cyYmCowfD0yY4MG8eS7l\ngz893QM/TjobkV5/3fetU1UFJCaGoJhBwmj0vpYmTPAdpnI6ga+/lnDqlIzjx2XltfnBB1HKlxGd\nzhsW5wfF2LEejB7t4Z5TA4ihMMS1t3s//A8ckHHggAZHjnjfbDU1XavyI0d6cNllHsyZ41LecGPH\nupGY2HnKXWfEnGZYLa2tvqtKVmsIChkioqKAzEyBzEw3Zs/uWssQAqirk3D0qKy8do8elbFjhw71\n9d7XsFYrMHasB5MnezBpkhuTJnlvzeZQdRPeGApDyJkzkvLhf+CAjEOHvGsBLpf3AyojwzvGv3hx\nh/Lhn5npQWxsiAsPQ7fd1oH/+3+7v31yckJUzBAmSTi395kbM2Z0H5JqaJBw7JiMgwe9X3r279dg\n82Yt2tq8r/dRo7qHxOTJHowaFYouwgtDYRDyeICTJyXs3avBl19qlDdF5zenmBiBCRM8mDbNjTvu\n6EBWlhsTJvDDfyDdeKMLq1Z50NHhXSapqW4kJIS4qDCTkCCQkOBGfn5XWLhcwPHjMvbv94bEl1/K\n+O1vo2C1eoPCYhHIzQUmT9YhO9uNyy93Y8SIIXsdsZBgKISYEEBVlTcA9u7VYN8+Gfv2aZTd/0aN\n8mDiRDd++MMOZGV5kJXlxqWXCp6IM8SOHJHR0dE1hFRTI8PlAjeU9jOtFhg/3oPx4z1YtMh7VSMh\ngJoaCfv3yzh4UIuDB6Pwpz9psW6dd0NEWpoHublu5Oa6cfnlHkyZ4obReLFniWx8CQ8gIYBvvpGU\nD/+KCg327dMo33IuucSD7Gw37rnHialT3Zg61Y1hw0JcNPXo97/XATh/u4KEAweAqVNDVVHkkiSc\nO6jOje9/X8BiiUJjowNVVQJ79miwe7cGe/bIeOIJPRwOCRqNd03bGxJu5OZ6t7nxi5YXQ6EftbUB\n//ynjC++0ODzzzX44gsNzpzxvvKSkjzIyfGguNiJnBw3pkzxIDGRq7lDRXKyC0DUBdNCUwv56goK\nF264wbtG4XIBhw7J2LNHgz17NPjsMw3+8AcdhJBgsQhMn+5Cfr53uGryZM+QPLBSDQwFFZ05I2HX\nrq4A+Oc/ZTidEoxGgZwcN37wgw7k5nrXBlJSGABDmdns+7WyrS0EhVCfabXA5MnevZh++MNzR6I3\nAxUVGuzc6Q2JX/3KuzZhNApccYUbV17pDYncXHd4H5B5HoZCgNxu4MABGZ9/3hUCVVXeD4qRIz3I\ny3Pjxhs7kJfnxsSJPCgn3Iwd64b33ESdQ0gCKSkhLIgCYjIB117rxrXXejdmO53Avn0ydu7UYudO\nDV58MQq/+pUEnU4gO9uD/HyXEhThumMHP6r6qKPDOxT0ySdafPqp91tFS4v3xTJligff+54LeXlu\n5OW5kZzMtYBw9+mnWly4TaG6Ghg9OlQVkRqiooC8PA/y8pz4yU+8X/4OHZLx2WfetYk33tDhN7/R\nQ6sVyM11Y+ZM78/ll7vDZriJodALp9O7Wvnppxp88okGu3ZpYLd7Vyvz8ty4+24nrrzSjexsN6Kj\nQ10tDTSbLdQV0EDQaHDuOAgPli7tgBBAZaWEjz7S4qOPNNiwIQpr10qIiRG46io3Zs504dpr3Rg3\nzjNkTwXDUDhPRYWMHTuA//3faOzaJcPhkBAbKzB9uhv33efEVVe5MHVq5G6Aoi4Oh+873uUKQSE0\noCQJGD1aYPToDtx+ewfcbuDLL2V89JEWf/+7Bo88oseqVRKSkjzn1iJcmDnTjdTUoTN6EPROWC6X\nCyUlJdBoNCgvL+82TwiBhx9+GKNHj4bZbMbVV1+Njz/+WJlvtVpx0003ITk5GWlpaSguLkZ7e3uw\nJQXEZgPmzInGU095L3u4cmU7ystbcfSoDX/8owMrVjhxxRUMBPKSpKHzJqf+o9EA2dkerFjhxFtv\nOXD0qA1//rMdS5Z04OhRGffcE43s7FjMnGnEww9H4dNPNejoCHXVFxfUmoLdbsesWbOQlZXV4/yn\nnnoKGzduxF//+ldkZmZi9erVWLBgAU6dOoXY2FgsW7YMHR0dOHToENrb27Fo0SKsXLkS69atC6as\ngMTGAseO2TFqVAyam9t5LiC6qIwM31DgMCIZjcB117lx3XVuAE40NEj4+GMNtm/X4k9/0uG55/Qw\nmQSuv96F73zHhVmz3INuV/Sg1hRsNhuWLl2K0tJSiB4uT6jT6bB27VqMHz8eOp0OJSUlaGxsxP79\n+1FbW4vNmzdjzZo1sFgsSE5OxqpVq7Bx40a43b6n5R0Iw4aB53mnPomLE+h+ZTSB+PhQVUODVUKC\nwPz5LvzmN23Yv78V27a14s47nfjqKxk/+YkBkybFYN48I9aujcLevTI8g+C7aFBrCklJSSguLu51\n/ooVK7r9XVVVBUmSkJqair1790Kr1XZby8jNzUVLSwsOHz7c69oH0WDw979rcOHeR8eO8aR41DtZ\nBnJyPMjJceL++52orZXwwQcavP++Fi++GIUnntAjKcmD2bPdmDvXheuuc4XkdBwDtqHZ6XSiuLgY\nt956K0aNGoUdO3bAfMG5b+PPfdWqr6/v11pkWYIs+24o1GjkbrfhLtL6BdTrOSVF8gkAnU4edMej\nRNoyHkr9pqYCRUUeFBU50dHhxGefySgv16C8XIs//lEHg0Fg9mw3CgpcmDu391OFq93zgLyEbTYb\n5s+fD51OhxdffFGZ3tOQ00CIj4+BdJH9xUymCDl08ZxI6xcIvudXX+1pakxQj9mfIm0ZD8V+b7jB\n+wMAR48C//M/Et5+W4sf/1gLrRaYNQtYuBCYP7/nU6qo1bNfoVBWVqYMF0mSBLvd/q3/pr6+HnPm\nzEFmZibKysqg13uvYpuYmIimpiYIIZQP6IaGBgDeYan+1NjY2uuagslkQHOzA273IBjc62eR1i+g\nXs9//rMWL72kv2Ba66A7fXakLeNw6TcxEbjzTu/P6dMS3ntPg7/8RYv//E8Z//EfwLRpHnz/+y4U\nFLiRmSn12rPF4v8XFb9CoaioCEVFRX2+f3t7OwoKCpCXl4eXXnqp27ycnBwIIbBv3z5kZ2cDAHbt\n2gWLxYJx48b5U5bfPB4Bj6f3tRS32xNRex9FWr9A8D1v2yahoqL7tJoaz6C9GlikLeNw6nfECOCO\nO9y44w4nGhuB8nIt3n1Xh8cei8Ivfylh0iQ3nnzSeyS2Gj3368Db2rVrodfrfQIBABISErBo0SI8\n8MADaGhoQHV1NR555BEUFxdD5jlsaZCrr/dd0xzs+5/T0BcfD/zgBy784Q8OHDpkw4YNDkycKHD2\nrHrPEdSnb1lZGQwGA4xGIyRJQmFhIYxGI5YvXw4A2LhxI3bu3Kncp/N29erVAID169fDZDIhIyMD\n2dnZyM/Px6OPPhp8V0T9rKe9pvswmkqkmthYoLDQhfXr27FkiXqPG9SG5m8bTjp+/PhF/73JZMLr\nr78eTAlEIdHTThLhetZMiiwcpyEKwMmTvm+dw4dDUAiRyhgKRAG45BLfNYWRI0NQCJHKGApEAbjl\nFt9Too4ZE4JCiFTGUCAKQFOT7zSeOpvCAUOBKAAjRlw4RUCv7+meREMLQ4EoALm5nddo9pJl7n1E\n4WGQnb6LaGjYtEkLxNQCo3YAbh08p65HfX3P56QhGkoYCkR+anO14a2Oe4D7ygDNucOY2+PwTMU9\nWPPd+0NbHFGQOHxE5Kd7P7wLx+I2dgUCAOhbUFr5KJ7d81ToCiNSAUOByA8nzh7D/xx7q9f5z1c8\nA4fLMYAVEamLoUDkh+1flUOg9zPsWtut+Pxfnw1gRUTqYigQ+cHThwtDeUR4nLKZIhNDgcgPMy+5\n7qLzY3VxuGJE3sAUQ9QPGApEfsgaPgn/duncXucvnXwnYqPiBrAiInUxFIj89NK//Q7mmkJAnHeh\nHbcO30/4Mf7f9FWhK4xIBTxOgchPsVFxMP31LTTZTwDpfwM8OuDY97Dk+Rj0cOlvoiGFoUAUgPHj\n3fj6/XGocRaIAAANxElEQVRAQ9f1xMeMaQlhRUTq4PARUQB0Ot9VAl6Ok8IBQ4EoAE6n77Q+7K1K\nNOgxFIgCEB3tOy0qauDrIFIbQ4EoADU1vm8dmy0EhRCpjKFAFICODt+xIol7HlEYYCgQBUCj8Z3W\n2jrwdRCpjaFAFID4eN81hfj4EBRCpDKGAlEANBrfUOCGZgoHQYeCy+VCSUkJNBoNysvLu80TQuDh\nhx/G6NGjYTabcfXVV+Pjjz9W5qenp0Ov18NoNMJgMMBoNGLBggXBlkTU744f9x0/OnEiBIUQqSyo\nULDb7ZgxYwasVmuP85966ils3LgR7733Hurr6zFnzhwsWLAAtnO7aUiShO3bt8Nut8PhcMBut2PT\npk3BlEQ0IC67zPf02JmZISiESGVBhYLNZsPSpUtRWloK0cOROzqdDmvXrsX48eOh0+lQUlKCxsZG\n7N+/X7lPT/+OaLDjUBGFq6BCISkpCcXFxb3OX7FiBW688Ubl76qqKkiShLS0NGXaunXrMGbMGJhM\nJixevBh1dXXBlEQ0IFwu32kOXoWTwsCAnRDP6XSiuLgYt956K0aOHAkAyM3NxbRp01BWVgar1Yrb\nbrsNS5YswYcfftivtciyBLmH01lqNHK323AXaf0C6vX89NMdOHNGB0A693gCeXmD7/8x0pZxpPUL\nqN+zJFQav5FlGVu3bsWcOXN85tlsNsyfPx8ejwfvvfceDAZDj49x+PBhTJw4ESdOnEBGRsa3Pmdz\nczPMZjOamppgMpn6XKsQAhKPNCIi8uHXmkJZWZkyXCRJEux9OC1k5wbmzMxMlJWVQa/X93rf9PR0\nAEBNTU2fQiFQjY2tva4pmEwGNDc74HaH/3V2I61fQN2ea2slvPuuFno9sGBBB4xGlYpUUaQt40jr\nF7h4zxZLjN+P51coFBUVoaioqM/3b29vR0FBAfLy8vDSSy91m1dVVYXHH38czzzzDHQ6HQDg4MGD\nkCQJo0eP9qcsv3k8Ah5P7ytIbrcHLldkvKCAyOsXCL7n5mZg+vRYNDV5v1w8+6wGH388eM+dHWnL\nONL6BdTruV8H3tauXQu9Xu8TCIB3I/U777yDn/3sZ7Db7aipqcF9992HwsJCpKSk9GdZREF77z2t\nEggAcPSoBrW1ISyISCVBhUJZWZly0JkkSSgsLITRaMTy5csBABs3bsTOnTuV+3Terl69GtHR0di2\nbRuOHDmCtLQ0TJo0CWPGjMGrr76qSmNE/Smmh7Vy7qZK4SCovY++bTjp+PHjF/33WVlZ2LZtWzAl\nEIVEbKwbgEDn3keAQGxsCAsiUknk7LdFpKJXXolCVyAAgIQvvwxVNUTqYSgQBeCbb3z3XuPBaxQO\nGApEAUhN9d17bTDukkrkL4YCUQCKi53wblPoJDB1aqiqIVIPQ4EoALW1Mi7cpnD2bKiqIVIPQ4Eo\nAJWVvm8dhgKFA4YCUQCuvbZzl1QvSRI4d55HoiGNoUAUgEmTPNCcd/E1s1lAO2DnHCbqPwwFogA8\n9JAObnfXNoWzZ2WcOhW6eojUwlAgCsChQ77XaD5zJgSFEKmMoUAUgCuv9L302qWXhqAQIpUxFIgC\n0NTk+9bhEc0UDhgKRAFoafE9zUVP120mGmoYCkQBGDbM9zQX0dEhKIRIZQwFogAMH+4bCrzsN4UD\nhgJRAMaP973sodkcgkKIVMZQIAqA3edyzAJudygqIVIXQ4EoAK2tvmNFHD6icMBQIArAwYMXvnUk\n1NaGpBQiVTEUiALgdvtuaNbpQlAIkcoYCkQBsNt9x4psthAUQqQyhgJRAK644sK9jwQyMkJSCpGq\nGApEAcjMvHD4SILTGZJSiFTFUCAKQGamB9HRXcGQmurhcQoUFoIOBZfLhZKSEmg0GpSXl3eb53Q6\nsWLFCqSmpsJsNmPatGnYunWrMt9qteKmm25CcnIy0tLSUFxcjPb29mBLIup3kgR0dHT9bbNJ3CWV\nwkJQoWC32zFjxgxYrdYe569cuRJffPEFdu/eDavViptvvhkLFy5E7bl995YtWwaHw4FDhw5h9+7d\nOHToEFauXBlMSUQDYssWTbeL7DQ3Szh9OoQFEakkqFCw2WxYunQpSktLIYTvLnqzZ89GaWkpUlJS\nIMsyli5dira2Npw4cQK1tbXYvHkz1qxZA4vFguTkZKxatQobN26Em4eG0iA3apTvhuaEhJCUQqSq\noEIhKSkJxcXFvc4vKCjAhAkTAADNzc1YvXo1xo0bh9zcXFRUVECr1SIrK0u5f25uLlpaWnD48OFg\nyiLqd07nhWNFEjjySeFgQC41PnfuXLz//vuYOnUqNm/eDL1ej8bGRpgv2DIXHx8PAKivr+/XemRZ\ngiz7DgBrNHK323AXaf0C6vUcGysjJ8f3sbUD8o7qu0hbxpHWL6B+zwPyEt62bRtsNhteeOEFXHPN\nNdi3bx8A9DjkNBDi42MgXWSroMlkGMBqQi/S+gWC73nxYu9PdzFBPWZ/irRlHGn9Aur17FcolJWV\nKcNFkiTB7nuqyF7Fxsbi/vvvx+9+9zu8/vrrmDJlCpqamiCEUD6gGxoaAHiHpfpTY2Nrr2sKJpMB\nzc0OuN2+p0YON5HWL6Bez3/7mwaPPNL9qjpvvdWKYcOCrVBdkbaMI61f4OI9Wyz+f1HxKxSKiopQ\nVFTU5/vn5ubi4YcfRkFBgTJNlmXodDrk5OTA4/Fg3759yM7OBgDs2rULFosF48aN86csv3k8Ah5P\n72spbrcHLldkvKCAyOsXCL5nu11CRUX3aR0dnkF7Sc5IW8aR1i+gXs/9OvCWn5+PVatW4eTJk3C5\nXHj55ZdRWVmJefPmISEhAYsXL8YDDzyAhoYGVFdX45FHHkFxcTFkOXLGA2loMpsv/FIhYDSGpBQi\nVQX16VtWVgaDwQCj0QhJklBYWAij0Yjly5cDAJ588klcf/31mD59OuLj47FhwwZs2rQJl112GQBg\n/fr1MJlMyMjIQHZ2NvLz8/Hoo48G3xVRP6up8T11th+jqUSDliRCtbVXBc3NzTCbzWhqaoLJZOrz\nv6ura+lxulYrw2KJgdXaGhGrnpHWL6Bez888o8Vjj3XfsFdR0YK0tGArVFekLeNI6xe4eM+JiXF+\nPx7HaYgCcPSo71unuTkEhRCpjKFAFIDJk32/hfKIZgoHDAWiAPzoRy4kJHQFQ3a2C/28JzXRgGAo\nEAXA7QaSk7s2x40cOWQ3zRF1w1AgCsDevRocOKBR/t6yRYemphAWRKQShgJRAJKSPJDlrrWDuDge\np0DhgaFAFIDMTIEnnmhHUpIHl17qwe9+54BOF+qqiII3yM7pSDR03HZbB777XRf0egE/DpMhGtS4\npkAUoJ//XI+srFiMHRuLV17hagKFB4YCUQD27ZOxYUMUAMDjkbBypR4OR4iLIlIBQ4EoAG++2X3k\nVQgJdXUhKoZIRQwFogD0dCLftraBr4NIbQwFogDcc4+z2y6psbECY8eGsCAilTAUiAIQHw/s2tWK\nefM6sHixE//8py3UJRGpgrukEgVo1CiB3/+eY0YUXrimQERECoYCEREpGApERKRgKBARkYKhQERE\nCoYCEREpGApERKRgKBARkYKhQERECoYCEREpgg4Fl8uFkpISaDQalJeXd5vndDqxYsUKpKamwmw2\nY9q0adi6dasyPz09HXq9HkajEQaDAUajEQsWLAi2JKIBUV0tYc2aKDz9dBRsPPURhYmgzn1kt9sx\na9YsZGVl9Th/5cqV+OKLL7B7926MGDECzz77LBYuXIhTp04hKSkJkiRh+/btuOaaa4Ipg2jANTUB\nBQVG1NR4v1dt367Fu+/aQ1wVUfCCWlOw2WxYunQpSktLIYTwmT979myUlpYiJSUFsixj6dKlaGtr\nw4kTJ5T79PTviAa7/fs1SiAAwOefa2C1hrAgIpUEFQpJSUkoLi7udX5BQQEmTJgAAGhubsbq1asx\nbtw45ObmKvdZt24dxowZA5PJhMWLF6OOl6+iISA93YPo6K4vNCkpHpjNISyISCUDcursuXPn4v33\n38fUqVOxefNm6PV6AEBubi6mTZuGsrIyWK1W3HbbbViyZAk+/PDDfq1HliXIsuQzXaORu92Gu0jr\nF1Cv50svBTZvbsNrr+kQHQ0sX96BqKjB9/8Yacs40voF1O9ZEiqN38iyjK1bt2LOnDk9zrfZbHjh\nhRfw5JNPYt++fUhOTva5z+HDhzFx4kScOHECGRkZ3/qczc3NMJvNaGpqgslk6nOtQghIkm8oEBFF\nOr+ipaysDAaDQdlTyB+xsbG4//77YbFY8Prrr/d4n/T0dABATU2NX4/tLwYCEVHP/AqFoqIiOBwO\nOBwO2O3fvqdFbm4u/vKXv3R/QlmGTqdDVVUV7rrrLnR0dCjzDh48CEmSMHr06D7VExcXh6amJsTF\nxfnTBhER9aJfB97y8/OxatUqnDx5Ei6XCy+//DIqKysxb948JCUl4Z133sHPfvYz2O121NTU4L77\n7kNhYSFSUlL69PiSJMFkMvGbPxGRSoIKhc7hJKPRCEmSUFhYCKPRiOXLlwMAnnzySVx//fWYPn06\n4uPjsWHDBmzatAmXXXYZoqOjsW3bNhw5cgRpaWmYNGkSxowZg1dffVWVxoiIyH+qbWgmIqKhL3L2\n2yIiom/FUCAiIgVDgYiIFAwFIiJSMBSIiEjBUCAiIgVDgYiIFAwFIiJSMBSIiEjBUCAiIgVDgYiI\nFP8fMnmPjabeAtEAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 4 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sample_plt = list_plot([ (beta_1, y) for y in logL200])\n", "logL_pt = point((beta_1, mxLogL), size=40, color=\"red\")\n", "logL_mean_pt = point((beta_1, mean(logL200)), size=40, color=\"green\")\n", "(logL_plt + sample_plt + logL_pt + logL_mean_pt).show()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 実際の解析では不可能ですが、今回は真のモデルが分かっているので、評価用のデータを生成し、E(logL)を計算できます" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

繰り返しは簡単

\n", "\t

\n", "\t\t次に上記の処理を繰り返した図４．９(B)を試してみます。\n", "\t\tSageでは、一度試した処理をコピー＆ペーストしてループに入れたり、関数にすることで結果を確認しながら\n", "\t\t処理を進めることで効率良く目的のグラフを得ることができます。\n", "\t

\n", "\t

\n", "\t\t最初に空のGraphics=gを生成し、上記の処理と同じことをして、プロット結果をgに加えていきます。\n", "\t\tたったこれだけで、図４．９の(B)が再現できます。\n", "\t

\n", "\t

\n", "\t\tRのseedをセットすることで、久保本の平均対数尤度に近い曲線が表示されました。（2016/07/17）\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAECCAYAAAARlssoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzsnWd4VEXbgO+zJZtNTyiGIghIR0wQECkKSu8iWBBeQAmo\nKDYU8RUVG36KYIMXUEAEETuIIKAoFqpUKQm91xRIL1vO92PJJqfsZneT3RBy7uviCnvmlJlT5pl5\n5imCKIoiGhoaGhoagK68K6ChoaGhcfWgCQUNDQ0NDSeaUNDQ0NDQcKIJBQ0NDQ0NJ5pQ0NDQ0NBw\nogkFDQ0NDQ0nmlDQ0NDQ0HCiCQUNDQ0NDSeaUNDQ0NDQcKIJBQ0NDQ0NJ5pQ8II1a9YQGxvL0KFD\n3e5ntVp5+eWXadCgAeHh4XTt2pVjx445yzt37kxQUBAhISGYzWbMZjPx8fH+rr5fOHnyJIMGDaJq\n1arUqFGDUaNGkZGRUd7V8hpv2vHhhx/SpEkToqKiuP3229mxY4ezrKI+W2/an52dzbBhw9DpdBw8\neDDANfWesmpbRX223qIJBQ959913eeqpp2jUqFGJ+06dOpVFixaxfPlyUlJS6NChAwMGDHCWC4LA\nvHnzyMnJITc3l9zcXHbu3OnP6vuNfv36ERMTw6lTp9i+fTv79u1jwoQJ5V0tr/G0HStWrGDKlCks\nXryYCxcu0LdvX/r27Utubi5QcZ+tp+0/d+4ct9xyC0ajEUEQyqGm3lNWbauoz9ZbNKHgIWazma1b\nt9KgQYMS912xYgUJCQm0aNECk8nEq6++SnJyMlu2bHHucy3EIUxPT6dNmzZMnToVs9lMzZo1GTFi\nBH/++Wd5V80rvGnH3LlzGTVqFK1bt8ZkMvHcc88hCAIrVqxw7lPRnq037U9OTubdd9/l1VdfrRDt\nLOu2VYQ2lxZNKHjI448/Tnh4uMf7Fx9pCIJAZGQku3btcm5bunQpzZs3JyIigu7du3P06NEyrW8g\niIyM5NNPP6VatWrObSdPnqRWrVrlWCvv8aYd27dvp1WrVs7fgiAQFxfHP//849xW0Z6tN+1v2bIl\n/fr1C2T1SkVZt62iPVtf0ISCH+jbty9z5sxh7969FBQUMGvWLE6fPk1aWhoAzZo146abbmLDhg0c\nP36cqlWr0rNnT6xWaznXvHRs27aNjz/+mJdeeqm8q1Iq3LUjNTWV6OhoybaYmBhSUlKAa+PZXivP\nUY3StK158+YV/tl6hKjhFSNHjhQfeOABt/vk5eWJTz31lFizZk0xNjZW/O9//yv27t1bfOedd1T3\nz8zMFI1Go/jbb7/5o8oB4e+//xZjYmLEjz/+uLyrUipKakdQUJC4cuVKybZhw4aJI0eOVN2/oj1b\nT5/j8ePHRUEQxAMHDgSoZqWnrNtW0Z6tpxjKWyhdi5hMJmbMmMGMGTOc21q2bOlSrRIWFkZMTAxn\nz54NVBXLlBUrVjB8+HBmzpzJgw8+WN7V8RlP2lGtWjVSU1Ml21JTU7nppptU969Iz/ZaeY5q+KNt\nFenZeoOmPvIDO3fu5Pfff3f+PnPmDPv376d9+/ZkZmYybtw4zp8/7yxPSUkhOTmZ+vXrl0d1S8XG\njRsZOXIk3333XYXuSDxtR+vWrdm+fbvzt91uZ8eOHbRr165CP1tfnmNFsT4qi7ZV5GfrNeU9Valo\nuFIfNWnSRNywYYMoiqL42WefibGxseLhw4fF9PR0ccCAAeLgwYOd+7Zq1UocPHiwmJaWJqalpYlD\nhgwRW7VqFbA2lBVWq1Vs1qyZ+Mknn5R3VUpFSe0o/mxXr14tRkdHi5s3bxZzcnLEKVOmiHXr1hXz\n8vJEUayYz9ab9hdy7NixCqE+Ksu2VcRn6wuaUPCQ4OBg0Ww2iwaDQTQYDM7fheh0OnHNmjXO3xMm\nTBCrVKkiRkVFicOHDxczMjKcZadOnRLvuecesWrVqmJ4eLg4aNAg8ezZswFtT1nw119/iTqdTjSb\nzc77Ufj35MmT5V09j3HXjhMnTiie7ezZs8U6deqIZrNZvP3228V9+/Y5yyris/Wm/W+88YYYHBws\nBgcHizqdzrnfm2++Wc6tUKcs21YRn60vCKJYCQxvNTQ0NDQ8QltT0NDQ0NBwogkFDQ0NDQ0nmlDQ\n0NDQ0HCiCQUNDQ0NDSeaUNDQ0NDQcKIJBQ0NDQ0NJ5pQ0NDQ0NBwogkFDQ0NDQ0nmlDQ0NDQ0HCi\nCQUNDQ0NDSeaUNDQ0NDQcFLp8imIokhmZmZ5V0NDQ0OjXAgPD3cb9rzSCYXMzEwiIyPLuxoaGhoa\n5UJ6ejoREREuyytdlFRtpqChoVGZKWmmUOmEgoaGhoaGa7SFZg0NDQ0NJ5pQ0NDQ0NBwogkFDQ0N\nDQ0nmlDQ0NDQ0HCiCQUNDQ0NDSeaUNDQ0NDQcKIJBQ0NDQ0NJ5pQ0NDQ0NBwogkFDQ0NDQ0nFVoo\niKJIRkYGmlO2hoaGRtlQoQPiFQa3KynAk5zkZP/GPtLpBGJiQklLy8Zuv/YFltbea5fK1Fa49ttb\nrVp4iftU6JnC1YpOJyAIAjqd66BT1xJae69dKlNbofK1Vw1NKGhoBJiFC420bx9C794hJCZqn6DG\n1UWFVh9paFQ0du/W8fzzJkTRMRIdOdLMli3Z5VwrDY0itGGKBhkZMHy4mWbNQhk5MpisrPKu0bXL\nyZM6p0AAOHVKQM1OwmKB9ev1bN2qfaIagUV74zT4v/8zsWaNgZQUHatWGXn//aDyrtI1y2232YiN\ntTt/DxxoRZ7vxGKBIUPM3HtvCH37hjJ5sinAtdSozGjqIw3OnRNkv7Wxgr+oWlVkzZocli0zEB0t\nMmSIVbHPP//o2bix6NOcMyeISZPyCQkJZE01Kiul/vqtVisTJkxAr9ezdu1aRfmFCxfo0aMHOp2O\ngoICSdmlS5e47777iI2NpVatWiQkJJCfn1/aKvlEdrZD35uaWvmsDh54wIJe79BhGI0i995rKeca\nXdsIAuh0jn9qWRHDwqT6JLNZJEibvGkEiFIJhZycHDp27MilS5dUy/fu3Uvbtm2pXr26ak7Q0aNH\nk5ubS2JiItu3bycxMZGJEyeWpko+cf68QOfOoXTrFkqbNqFs3qwPeB3Kk27dbCxfnsPEiXmsWJHD\nHXfYyrtK1yypqQI9e4YweXIwjz9uZvz4YMU+LVvaeeaZfHQ6kZAQkY8+ysOgzek1AkSphEJWVhYP\nP/ww8+bNU/UqvnjxIl999RWjR49WLVu+fDlTp04lOjqa2NhYJk+ezIIFC7DZAtspzZtn5MQJx63I\nyhJ4993KNSw7ckRgzBgz//d/wYwZY+bkyco3WwoUGzfqOXu26LP7/nuD6kLzCy8UcOpUFkePZtG/\nv1LFpKHhL0olFKpXr05CQoLL8jvvvJN27dqplu3atQuDwUDz5s2d21q1akVmZiZJSUmlqZbXyEdh\nlW1U9t57JmdHdfKkjg8+qFxCMZBcf71d8rt2bVFVhQRgNDpUTBoagaTcXrnU1FQiIyMl22JiYgBI\nSUkJaF0SEgpo1swxO6la1c5LL5XPukZ5YZf2U6ojV42yIS7Ozttv51Gvnp1WrWwsWJBb3lXS0JBQ\nrmPi8gpkp9NJ3dirV4f16/NITRWIihIxmaA08lKv10n+Xu0884yV06cNWCwCQUEi48dbMRg8r3tF\na29pkbc3IwMOHNBRs6ZIrVolv9NjxtgYM6a4MFDetz//1LNkiZHgYJFx4yw0bGhX7BMIKvuzrYx4\nJRQWL17sVBcJgkBOTo7PF65WrRrp6emIouhchE5NTQUcail/EhMTqrrwXa1a2V4nIsJctif0Ex06\nwJYthb8EwDfbx4rS3rKisL3R0VC3btmee8AAxz8H5a/PrKzPtjLi1ds2bNgwhg0bViYXjo+PRxRF\ndu/eTVxcHABbt24lOjqaxo0bl8k1XJGWlu3XgFd6vY6ICDMZGbnYbOUzwvOGKVOCOPXnKWpwnrPU\noFH32kycWFDygVeoaO0tLcXbO2+enkWLitZgmja18fHHeZL9d+3S8eyzRZ1M3bo25s+X7lOcHTt0\nPPectFNasSLbpZ/CsWM6xowJxm53vNNVqtj5+uuyUUt5+2xtNhg0KISsLEddBEHk00/zuOEG79+L\n3bt1nDmjIz7eRo0agdEqFLZ36FBISnLUf8aMPG666dp4r6OjQ0vcJyBDkEI1UXF1UZUqVRg8eDAv\nvfQSCxcuJDc3l9dff52EhAR0fl5ds9vFgITFtdnsWK1X98uk37+Pl1Y9SdVDW53bUvLaI/b/EFvD\nRl6dK9DtPXdO4NAhHc2a2alaNfCqSJvNTnq6wM6dRdsiIkTFPUhM1En2OXhQ59wnN9ehKoqIcHg7\nA0RH29m5U8Qxa4OgIJGgIDtWF0ZIBw8KbN9efJCjIyvLTrDS2tXJjh06li0zUquWnYcftpRoXOHp\ns83Kgr/+Kl4XgePHRWrXlh67fr2er782UqOGnaefLiAsTHqeefOMTJrkaEBEhMjq1dnceGPgnnFS\nEleemcCvv+po2rTyWICVqvddvHgxZrOZkJAQBEGgf//+hISEMHbsWADGjBmD2WymZ8+eAERFRRES\nEsIXX3wBwOzZs4mIiKBevXrExcXRrl073njjjVI2ScNTdKdPETWoj0QgAFRN3EjUwN4IFy6UU81K\nZssWPe3bhzJ4cAgdO4Zw4ED56IBHjrQ4O7zgYJFnn1UaKXTubKNGjaJOsW1bG506hdC7t5n27UMZ\nPjyEAQNCGDfO0Ql+9pmRQoEAUFAgkJzsug5t2tioU6fo/H36WNwKhMREHQMHhjB7dhCTJwfz/PNl\nF0YjLAzuuafI+bFFCxu33CI1Md+zR8fQoWa+/dbIRx+ZeOIJZWUXLTI6/5+RIbB8uVGxTyAQBJHW\nrSuX306pZgolqZPmzp3L3LlzXZZHRESwZMmS0lShbMjNxbRiGfojh7BXq07+3UMQq1Qp71r5HfPc\n/6FLS1Mt0yVfxLxgLjkvTA5wrTxj1iwj2dmOjjMtTccnnxiZNi3wVmM//mjg9GmHQMrLE5g/P4h2\n7aSqoWrVRNauzeHnnw2kp8Obb5oo3ukX8s03Bj74APR6ZZna5PnoUYG33zYRHi7y8st5TJ9uIixM\n5LXX3N+Hv//Wk5dXdI1ffzUAZXfvZs7Mo29fKzk50KuXFbNMPb9tmx6rtej6W7YonUVjY0X27y/6\nfd11gZ0J3n9/AXXqCAwYYHXO4CoL5b+CVc4Y//qDiIQRks4xbMpksl59g7yHx5ZjzfyPafVK9+U/\nr/SrUJg1y8j69QaaN7fzwgv5V6y+PEOuXw8tWVXqF44dk/bWx4+rz1iuu05k5EgLr7wShJpAKEQU\n4emnC5g7N4icHMd+jRrZkI9RLlyATp1CsVgc+zhG1o7/Dxtm5o8/XBuBBAd71sFu3arjxAkjjz3m\n0e5OdDro08e1uiU+3oZeL2KzOeqrNhJ/5508xo41c+yYQJ8+VoYODWzolYQEC6NGXd2qX39Ree2u\ncKhPIoffrxgtC/n5hE96jqB1ylhO1xSulNSFWPz3IS5ZYuDVV4NZv97AzJlBV0bPnjNpUj716zs+\n2ptusjF+vPrCeGKijvHjg5k8OYjLl0tdbQU9e1oJCirqZPv1c39P27e3AeqdcmFMJLMZli7NoWFD\nG/HxVr77Trlo/M03RqdAuHK083+JiXryXK9jS2YJoPRTAVi7Vk+/fiF8+qljEf2nn8ou9EtcnJ2F\nC3Pp18/CmDEFzJyprGydOiI//5xDUlI2772XrznxBZBKPVMwL5yPkOM6wYl51kcU3NU9gDUKLJb2\nHdF//aXr8g63++3af/8t7WT++su7TqdOHZFNm7LJzARX6bkvXBDo1i2EggJHJ7h2rYEtW3w3o1aj\nVSs7K1fm8NtvBho2tNO3r7pQyM6GTZv0XHedyNNPF7BggRGDAVJSino7g8EhGM6fFxg4MMRpTdSh\nQwiHD0vfU/d+C6LbhWPHekPRQrZahzttmkmS92HJkiB69iy7QUL37ja6d69caplAsGePjvXrDTRu\n7Pv9rdTy17Blk9tyYwnlFZ2cseMQjeoLeGJwMLkJj/jt2hkZ0tFqtg/JxwTBtUAA+PZbg1MgABw7\npvfqOunp8NBDwbRrF8qkSSZcheS6+WaHBY0rgZCVBV27hjB0aAjdu4diNsPBg9n88EMuglA0a4iI\nENHpYPZso1MgAGRk6Dh3TnrOHj1sDB5sQacT0emUMw/3MwUoPrNQ8yGVW+eVtZ/pjBlBtG4dSr9+\nZo4evfpibb3xRhA9e4YwZ075LHD7ws6dOnr3DuH1100MG+Z73Su1UChJiS2a3JhwXAPYbmpJxicL\nyTJGS7ZnmKqS/tkX2Br5z19EbvfdrJl3+lubDaZPD+I//wl2+fJnZso7GxG9FxOSV14x8dNPRo4e\n1TFvXhDz5vn2kS1dauTIkaILFwZcPHFCkIzGU1N1WCyoZr5TG83PmpXH+fNZHDyYhdlc1Gu3aGFT\nmHgWp0EDqQorPFzZ48tnIsUTA5XExo164uNDadIklAULlPfst9/0TJ1q4uRJHVu2GBg3Tt1R7OJF\ngX//1bkVcP7i99+N7NihZ/LkYNaurRhRk1euNJCfX/Q+LVumCQWvye/T3315X/fl1wIFvfuy/5cD\nPBHzORN5myeqfsGR9UlY7uzm1+s+8kiB05QzIkLkhRe8s36ZMSOIt982sXq1kcmTg1m4UPkBOPT3\nRQgCXgmFwsi5rn57yr590uMKl2puucVOaGhRh9y2rRWj0bFdikh4uOvz5+YKFJ/wuXJys1oddVmy\nRGryqrY4LrcY8jTBj80G991n5swZHWlpOiZONJGYKD3/yZPS32ozhWXL9LRsGUrXrqHEx4fiwkgu\nIBw8WDG6ybp1pcK9uJmyN1SM1vqJvCH3Y3UxGrZHRJL7xNMBrlH5UK9ZMBN3D+TebeN4YVc/ajfw\nf5TUH380Ok05MzIEFi/27pry0duvvyp7+44dbU49uCCITJmSjwttmSp3312kDtLrRZfqoZKQj7IL\nI6zs3q2TqLP++cfRBqXPhUBmpvK8//ufkSZNQunUKUSijtu61UCubG06Lw8GDTLTpUsoK1aUfBPa\ntpXOJho18kw/nZaGZLQKAn/+KX02clVUcfPUQp57zuxUoaWm6njhhfKbtastxF+NPPighUceKaBe\nPTs9e1p46y3fzIwrtVAgNJTL368kv99AxGIrc5ZbbyP9h5+w3diwHCsXWEwmx+JtoDJ87d8vffXk\no+mSkJuC7tunFAo6HYwfX0CHDhb69rVIOnlPaN3aRlSUowerU0ekUSPfeocuXaQdbPXqIsuXG5g+\nXTpiF0WBjAx1m3z5c/nnHx2vvGIiLU3H5cvyeycq1j9WrTKwebPBeZ2i+ojceadyAXn+fKnp7Nq1\nnknTqCiQW1cVWokVsmuXtL5q6jJZkka/WI55ijezy/JEp4PXXstny5ZsPv88jypVfFsIqtTWRwBi\n9epkzPsc4cIFcvcew3xDNcQGDcq7Wtc8nTtbmTevqKdzdJyeExzsWAguRE0vfuaMQJ8+IU69/caN\nBhITPV9pfustE5cvO449dkzH7NlG/vtfz2NCFbJzp57iHeyFCwIJCWp6dBGjsbjZaqF1kKhYUN+0\nSXpOOXKPZmXHJjj//vabEbnzmjzWpadZcuWL2KBcutu3T1qutojdr5+Fr78ufD8cFlvlRXmEUClP\nKvdM4QqZmdBrVD3qPNCNuMEtK4wOsSLTpo2NoCDHCFIQRDp08M7c0bEGITqPnzxZ2WvNnGlULOSq\nqWFcsWOH9D3YulV9yDhiRDD16oXRvn0IVwL9SlC+T646c4GCAjh1SifZx24XFIutjlmLemellpyn\nd28rdesWCl7pcWruKk2b+rbQHBKCzBpKxGiUXk8tQrGcQtXilSN8Xs8pCy5erFz9QeVqrQvmzAli\n2zbHB3/mjI4pU8ouFkxFwWp12McHKhOqw3/A8fqJosDAgd65JD/4oJXffsvhnXdy2bAhm65dlRVX\nejm7znKmRlqadGf5ginA008H8fPPjpAbhw/r6dxZ2Q7PzWBFLBawWpWdvVyvfe6cgCvh4jiHdNt3\n3+mLdawl3wS5AHTlqS0nNxeJOS0IijUDh6OhVJ0mZ/t26fW//rq8TENFqlevIIsKOHJ7bNqk58wZ\n3818NaEAikU5+e9rnRMnBG67LZSWLcPo2DGUs2f9bzdePE8xeK6eKE6LFnZGjrS6jJ7pCMQmLXNn\nqlkSamqOlSulndWFC8p750ninUJCQ+HSJeU55MLMMZtwhaAQIp984j68hpzkZOm+ubmeHRsWBtWq\nFV3cYBBp2FAqsOPipBFce/ZUTlXCwqT3zNf1nNIjsHlzxVhUOHtWoHPnUAYMCOG220L57Tff6q0J\nBWDECItzehwcLLoMmXCtMm2ayTmKPHJEx/vv+3+1uXhoCF/JzobDh5WqlUIci9fSUas3C5byBdKu\nXZWdl5rjmJzWrb1b4JaqThzIF2MzM91dV1Soj7xdLPV1xpiVBcnJRRe3WgUOHZKP+g2SUBtffaWc\nBcivX1JEFn9SmlF3IPn8c6MkOOP06b59x5V+oRkcliV//pnN3r16brjBTu3alWthSW7pIf/tDyIj\nRcmMzGDw7p7v3evw3szLE4iIsLN+fY7iuXXoIF2w1evFK9YxniEP/SSPGQSOrGvFbejVnMzWrPFc\n9SGKjiQ5sq3I0pkTFOSuoxKwWKQWS+78DNRUagaDtP2eqt0ci8pF9xxEhfps1y6pkFCbJco93stz\ntB4TUzH6A7lviTcBJoujzRSuEBXlsGuvbAIBoHdvizPcgiCIbiNclhVSvTN4o9oAx+JuYSedkaHj\noYeUduyO0NpF57XZXM8q1JCruArXnYpjs5X8vnjq+OXQv8PmzcqFaXlqi5L03PL8C0WLzErU1GLe\nrL0U59IlkM/ODh2StsdkKvmeya2nyitHNfjutBho6tWTPuPYWN/6sorRWj8jivDWW0F06hTCo48G\nq9pNX8v8/HORlY4oCqxa5f8JZGlj6Zw7J3115SoK4Iq3cNGF9HrRbfKZkpBb0QDOxfJC1NpV3GvZ\nPQ6T1KT9yh754gXpOeQjaTlydcvDD1twZa2kWhPR/W9XqAXik5u3njpVssR54omi6arBIDJhQvmp\ndCuKNaL8u5UHnfSUitFaPzN7tpH33zdx4ICe774z8thj13bMIzmldSTzBfmI3VsdtryDDglR9lqO\nj8L3mYK8ToV5h4sjXxBVUx/t3auMwaSOQNaeQ9xz9mPF/q1fvgchq8ie9sIFd89IpH596RZHCGzP\nh/++rimohaPYsUPaOZ08WXJndfBgUV2tVuXMJ5CU53qGN8jVm77WWxMKcCUWTBHr11cMa4Oy4vBh\n6WugZnpZ1pTW9FU+IlVTSZw/rzzOGysnuQWPMsCeci1EbUS9YYN8+OyqcxZp+NQD5FrlaxACuTsP\nEDbxWeeWktRHcuH077/ePVNfQzs4nOykNyE0VHqyBnXyGMw3vMUkXmAq9TiqOM+PPxa/BwJvvVV+\nZuK+egYHmptvln5UWuyjUiBfWC3MCFVZkI8opLFr/IPcTNNzvbsD+ahdbeScmancVpo4Nmodvjy4\nm9r5PVeVCVgOn+QC18nPQDgZmJZ9h3DxIgDx8e4aInD6tHTL7bcHxgHFsagsfTZhYUW/9fv38euJ\nJnzDvUzibabyIodpSOhLEyU3Sj5oOH++/L7J1q0rhp9Cofd9IWozW0/QhAI4c/0WEigHrsrMY49J\nHZj69/fOo1kZVE25j9LfxH3ymZJQiwtVPF+DK2rV8rxTOUdVgpDrzwUuE4VgsWBIciQudjivuUZu\nZdW5s3cvtSemtmqoWbw4z5WbS+T9gwhPOyUtRyRk7v8wz53l3CZPGRoXV34dc1nnkvAXR49Ku/PU\nVE0o+Iz8oftqeVFRKR6LH9TjCJU1jkidRTf6r79Kt7it9uHKFzjBO49tuW1/RISyY9LrS75X8pAR\n7qhBCr1YpdheHUf8DPGK992vv7q7X6LCSS8pyVv1kW8fgSOXtPSetGnjuOnBP3yL/vw55UFXMM/9\nn3OqJfdG9zRKqz+oKM6scmsjb2ffhWhCAejVSzpKLYzzX1kwGKQdgDfhpX1FvrjtD/WAmqDw5gOX\nq4LUpuPyEbnagMLzBT8Rg0GgCUks5kFOcj2HacB7PE0YWVjr1ccafwsAERHuOkmlCevq1YFZJ1Oa\npBap9gxbN7s9Vn/qJLqzZwDlKHfRovLLgNa8ecVQHcjfT1+iBIAmFAClPbqnLv3XCvIFKkd4CP8i\nH8X7I2b9rbfKg8aJxMZ6frxcqKittTz6qPTLa9FC7d55PvPK69aTeP7lQZZwPadpwFGe4X1s6Ml9\neKxT6hTX06shV3Xt3BmYT11NxeYcZHjgTSUGq2dhy8gov65Kr68Y/YHc81rNMMITNKGAUheXnl4x\nXoKyQr7QLvfk9QdyNY4/9LZ//SUPLy2UyrRRTfUkDxwnX3gGWLfO01GuQNq6vRhQXsiAjeDl3zt/\n//67+5G//BmqJbLxBw7Pa6kgdiTsgfx+A90eW3BbB8SqVVXLPFHT+YuKko7z8GHpM1aqTz1DEwpA\ngwbSYWphYpXKgtyOfNMm/zuvyZ3I/LGOoxbYr6zz/crvldqAwhv1UWyB0jyzEOM/W9AfOujyOsWR\n399AGU84HD+lgnj3bkc3Y+l4OwVd7lI9TjQayZn4X5fnDVTyJzX27q0Y3aQy2ZJvVIzW+plhw6Rf\n7X33BWCofBUh77QCEfsoP18qeP2hPoqMVAp3eXwYb1BzTPPkXtWv72mPLHAR9ZGysw5XdO7durl/\nR+Ve+XInO3+hNjotbiqZ/tkSFoY9RhZFK8k7iSN9ybdY2nd0ed5AvJOuqCjpOD33nHePJhSAXr2s\nzJ2by733Wnj55TwmTapcUVLlBMIELyvL/69eTIxyW2lGnGphLjyZooeEeDoNEonBfRhX+/XXA3Ds\nmHuVhnzm1bx5YHo2Nec1iUAym5kY+iE1OUtbttCERFqxA8sdXdyetzzNxMtzluINDz4o7bfq1dOc\n10rFiRNkUxhAAAAgAElEQVQ6TpxwZHgqaxWDhv9RCw196JCyM/Ym85octYVmT0w31dYZXJEd4nol\nvKBDJ2z1bwRKUhUo03eWn/oIp/qokNq17WQSwT+05QBNPDqvMmFS4KgoQsFikd5ntZmtJ2hCAVi6\n1MCbb5rYssXAwoVBvPpq5cu8Fmi6GX/nS+5nG7ewgr7cw7dezdMjIqSjUTUHsYsXlR12Sor3dXWH\nJx+e6PHUS+DCm++Th/L9y46MJWva+87ft9/uTn0kSPJXAxw8GJjF0uPHldv27pVe+8QJ7+uiNksL\nFL52roFm61ZpRU+f9m2hTsunAKxcKb0Nq1bpeffdcqpMJSDkrddYa5km2daXleSPHkDGJ595lBFG\nmetAuY9aX1zWummbTaQbv9CURJKpxnL6K/YJCfF8hpLX/i5uNfzDo9aZ3MEfFBDEMgZSd9J/6NOg\nunO/4rmnPSEjw6vdfUZNZRcSIhXYvjiDyVOjBpKK4swqH1P56oCoCQWUAfCKZ47SKFuMG/8m9P1p\nqmWmn5YTvHA+eQ8llHgeuSpHLc5RWJiyM/QmyY4c+WhVd/gwBxhCIw45t10mEt3St8m//0Hntptv\ntnls1mixwJnIpjyaOluy/a8OUqmycaP788nVRQZDYKJ9qgk/eZhvX+pRnqEmKkqYC3nGPl9Ny0vd\n+1mtViZMmIBer2ft2rWSMlEUee2116hfvz6RkZF06NCBv//+21l+6dIl7rvvPmJjY6lVqxYJCQnk\n++qGVwoslgoyFPATchNcfyYqD/58vtty8+cLfDqv3a78clu2VLajWjWfTg8ofa/CX3hGIhAAokgn\n/MnHMK7/zbnNmxg0BoN6roQtW6S/3QdtFBWTrejowPRsqanKbWlp0m7GF6HgbTrRsiQqqmKYH5WU\nY8NTSiUUcnJy6NixI5ccvu0Kpk+fzoIFC1i1ahUpKSl0796dgQMHknXFXm706NHk5uaSmJjI9u3b\nSUxMZOLEiaWpkk94kgnqWkZuQeNrdEVP0B9zbYcPoD/uvrwQ+VRZLTBdjx7y3kcZE8gb5PdF5+K9\nF0SRkI+K9P9y50jXONKFWq3K97F2benvpk3dh7mQ5zUoi5zYnqBm8isfafvSwZenWag84sHVSlmF\npylVa7Oysnj44YeZN2+e6mKa0Whk2rRpNGnSBKPRyIQJE0hLS2Pv3r1cvHiR5cuXM3XqVKKjo4mN\njWXy5MksWLAAW4Dtz0oTOfNaQN6h5uT4TyjYr3MfZ8JeXR422jPU9L716jneST1WCvMGB+rVMm74\n09mTudL/hyJbDUbAmHxWdf/t26W/t21z/+nKO4hAdWxqSXbkMwNf8hOUpwqnoiTZadBA+nKXi/VR\n9erVSUhwrf8dP34899xzj/P3yZMnEQSBmjVrsmvXLgwGA82bN3eWt2rViszMTJKSkkpTLa+pVk36\nxgVqVFUZybt/mPvyB9yXFyJ/RmrJ1fOX/sgG2mPFSD4mvmYIpqQ9nle2NOh0TkklT1wPMJQvCEFe\nIBLapw86lL2QPJNaSfG55LM/UYTqXKA7a2jPBgT8M/S+4kYhoWpVWZKdBt5/X+WpPqooC83Hj5ec\n28MTAjZGLigoICEhgeHDh1OnTh02bNhApCNQipOYK6YLKWVtNyhDpxPQ6Yqe9KhRNr7/vuitq1fP\njsHgu7zU63WSv1c78fHKbd6035v22vv1I//ZCZh+W6coszRtSsHj4z26dvPmMlWOTpAcZ/r2Kx46\n8zHEA8QTBAzhCLbJz5H97nRsLVqUeA1X98XZziaubewtHW7HYHS8U/JLhZHJJ8yiP12RD6yDqUFH\ncsiUmaWazTrJjLZDBzspKa57yvr1i+2fk8PKVv9HN3EtxisC5yw1+JAnWUt3Z7sKES5cID4+TNLE\nJk0g6PRJ7Dfc4PKaALGxyvt2882y8wuCYh/5M1eWe/dO+krhsy3+aI3GwFy7tDRtqpwJ+1JvQfTc\niNotOp2O1atX0717d0VZVlYWAwYMwG63s2rVKsxmM19++SVPPvkkF69kkgKw2WwYjUZ+//137rjj\njhKvmZGRQWRkJOnp6UTIvXXcIIoiQkUR/xoaGhoBxKuZwuLFi53qIkEQyPHAx79wgblBgwYsXrwY\n0xUTjmrVqpGeni7poFOvmC5Ur17d5fnKgrS0bMlMYdiwYM6dKxp1GQwia9b4GGIQx2gjIsJMRkYu\nNtvVb7lw110hSL1QRdat87z9vrS3tNd0d3zQsu8J/WiG2+PT5y0qcdTr6hqF7c08dprEwW/QVtzq\n3OMEdYh5+0msbdo6t3XvbsZmKxqxPcosEviUHvxEMjUl1/yDDtzNMtKQmklNmpRN165Fv0eODOLU\nKdcri599ls3114Nh53bCJzzlcr8DNOIBvnTeO9OXiwn5dA5dWUs/fqJFEyt3LRnN+KHJfJjUg9xB\ng8kb96TL8509C8OHS92Pb7zRwpw5RQ4iffoEk5dXfJajfPbKe29n3Tr/Z7spfLZDh0KhFluns/PL\nL1d/pp1evcw0LfiXGzjOJWLYSDvWrJPapUZHl+wa7pVQGDZsGMOGeabzBcjPz6dv3760adOGOXPm\nSMri4+MRRZHdu3cTFxcHwNatW4mOjqZx48beVMtr7HZRYsK4fr2g0MFaraXvzG02e5mcx9/s3Knc\n5ku9vWlvaa+5e7dUZxocLDqPD9q/X/0CxRBPncJau06p6miNimaM7SNydx+lCUkkU40twq2cj8+B\nYvvl59vYs6dIKGzGRAI7+Z1I5KbkqRxlF1bklp2bN9vp3Lno95YtOkUiHVrOgfB0SBrI6dM1qFED\nTF8tdXsvGrOTVN7Gaq0CQNi8T2DnTlZSCxMFOPRvo5mR1At27iTowkWyxj7h8ny//668XGKijpkz\ni+7Hli2CQt8tf/Zl9U76SlJSUR30+sBe2xf0Bw8wZ+tobhZ3O7ed5zp0P0yjoN8Ar87lV0XZtGnT\nMJlMCoEAUKVKFQYPHsxLL71Eamoqp0+f5vXXXychIQFdgP3KK1tSHTlyTVogbn9pF/PlnUpx9xZb\nY/fxdESDAVv9BqW6fiHHj+s5SGN+ZACbaI9dVN48efC/r7iPdCKwK8ZkIibyKFAJcyH3kVB152n9\nP+g2ER5vwqyzj2C1WxE88GAyGoqsoIRLjoB8JlmeaP2VhWkhy71rtprJr1xBXVGijhYSqAizviJc\nvkTkPf0kAgEglgtEjBmJcfNGr85Xqs9/8eLFmM1mQkJCEASB/v37ExISwtixYwFYsGABmzdvdu5T\n+Pett94CYPbs2URERFCvXj3i4uJo164db7zxRmmq5BMVxWPRX5RH+z1JeO8Nxc0483v1xebG9DW/\nd78STWM9xZOETBlBSdDrCUhoCyO6kN16ESP1czAaL8Etc+Hu4dAvAW5cTX5MHfJNVtDnQ43tUH0P\nICLG/is5Z16BDRovh+7Pwl2ToPZmsF1JyiuIrDyzhDc3T6GgQye3dTseCUef7cT8vZ84NhiK1Dqp\nZjggi+QtmmSJGmTI8ziA0hO8opmAl8bhMRAEL1mM/sJ51TLBZsP84XSvzleqx1OSOunw4cNuj4+I\niGDJkiWlqYJGGSAPgVCewcfKhKAgMuYvQug7hChRGoo6r2Fzst5+r+RzXL7EsEajSItbyeWwPKqn\nRWDaNhKsr4LB87CZK44sJ/Xeh0BfbMRebz2bbgshLygcwovpgG75lJbEULDpNWj5FYRescJLvZG3\nD1tpffpDbq/dmVOZJ8kbNQSqJRYd2+ltSJcKukW75/LsyAOENrgRw5HD7L4O/qkFYQXQ5yCEF8C0\n9iCaM3nhz2epElyFkXY758PgqZ7wfVNoUQt24Pj9xHGoV0LEAZsN0Fng+k2gL4AzbQCpHruiBJgr\nxGS6uqc2QeuVlnzS8t8cIz8PjWsqmMz2Fw7HpspKaKgoGfF6Ych11WJtcysdIvYzOH0Bt/MneQTz\nA3czdnZvGlZ1P9oV0y8xfupNfDW0WOCkOhkQ9yGNJ//MU2/LPMn0+RB+DnKjIb/o5l3KS+PxdWOk\nAuEKF6rkAMqF9UzS4LaZ0o1VHIOrJ9aNZfeIA/xn1QNSgVBIpHS0mCHmsjdtH7kLP2HEJ33ZVq/I\nLyIsHzqehNUNi/afsX0aA6tH0rmXcobwZ134+iHYvDIMd3mKfsv8BJ55B8KuCLuCUAr2jMZiexmj\n3rEwXp4Jc3zBl6iuVx1eWFtWMJntLyqvQAClo5Oas9VVSfBliDkEQep67ryIakzlRXqxmrtZxueM\nILSKe4EA8OW8BL6qqx5W9LXrD5H0+2LH+a150OMZmBALT9WD56vCvYM5etnRiX9zYCm51rKzWjmX\nfY4Fez9lX6rnDnjBR49x787xEoEAkGWSCgSA/al7mXKXQSEQnNcPh1d7hbi81vy9nzD/wrNFAgEg\nKJuCWz7g2T/GF22rsQMGjIJHboZRnaDNTHIsspfQmAPNv4I2s6DuH4o1lUBytXs0F9zVrVTlcrSZ\ngoYiIKA/w1w4id3l0KdXOQSZNWDXCOBWjw49k3kaBk+Bpt87RuEWM+y5n7S8l4kJruLczy6K0PhH\nqPMX2EyQOAi9vqGbMztYkP07uJEdPx1Zxm08yot/PQ+3FQvwp7dCs+/p+8Of/HzPbxxNP+JRe7zh\nl+OrPd637mW4mH3eKyHyZchBcDOS/yrqJP+nst1is/DeP2olV45LWsLTtzzHprMbIOEJ0BVTydT9\nm74/fMIPA34i0hTF0qQv4JlJYC5S/dkvteBg2gIaxfjXMtFJ24+g1l64cDMcHApcvbOFvAeGYf7f\nx+ivpGotjmgwkPPEM16dT5spaAScmTs/hEdaQZvZUH8d3LwYRnTjqd/GlZiQJjknmX4/9IAWXxWp\nZYy50GoBdy/rQ5bFEWzxVOZJzt59MzwwADpMg9vfhLG3MO7PYeTb3OvFD0a4t9g5aXBcY3eyuqln\nSm4KH2x/j+vcZFHzleJCryRe2RjEb4YTXp0/0+beVyRXVJcY2y/8Q3LuRdUyABGRr5O+ZMIfT0oF\nwhX2pvzLG5unsO7EWp787TGJQACwRO9l8Ir+ZBb4LzGEKIp8sOPKomz8Amg9F/qMI+/RG9hw5i+/\nXbe0iBGRpH+/gq20kWw/RW0y5i/G2tazwVYhmlDQCCjbzm9lyqaXVMuWJC3iy6TFbo+f++8sTmed\nUi1LTNvPl4mLEEWRYSvvQ6yi1Lv/mfwjUzaqX7+Q6wrcLyTH6MPdlgN8f+hbhjS+H30ZfmI6Qc+j\ncY+jF9yPWmunw4JlcH+LEQhe6F10go4Io/u2mfTq57PYSzZ93XZhK1a7a13Mtwe/Ysb2aYioDwzO\nZ5/jqyT/GabM3zuXHw//oCwITmf4qvtJzVWJC36VYKt/I7eyhVvYxjAW0Z013MAxCnr29vpcmlDQ\nCCif7Zvnvnzvp27Llx3+zm35D4e+44/Tv5OYts/lPksSF5ORL49QWsQD0Xe5vUb3xne7LQfIsWZT\ng0je/Fu9E1UZLDtwM1G6t/EDNK96Ew/fNEa1XBDh9XVw/H24O6o3Wa++Sdc6yrAzCrIcEQR61+vH\nfU2Gut21xw3qnUxc9XhCDO69ZYP07oVttiWLrec3u93nj9O/uy33FVEUmfvv/1yWZ1kyWZK0yC/X\nLkt2cAtfMIxf6I7dR5WXJhTQoqIGksOXDrotP1hCefYV9ZC78m3nt7rdJ8eazf5U10Jj9LA5xOWq\np2jrb29Ku3buO06AFlVbErrseyb+mstPX8BdRx0WP7GZMH4z7JkF/PQhpBZb4zh6J3dEDkV3saXi\nfLdG9uGd2x2qjdc7vE3ktjcgq1iY8YvN6P3Ng+j/eoHbxC3sfWMpBAdzZ51utKp+i+uKbngOZpwk\nJjiGF299mWdaP09sSA0okC0o2w1EBEXwesepqqcJD4pgRPOHXF7mtpodaF7lJtf1wPUspDiC4J8u\nKz3/MsfS3efy2Hlhu9vyawVtoRmoXdsuSYRSnmF6rwb8GSuwJJ14FbP78jhzQ37NdR1FNz6sMWaD\nawuZQswG14aV4cFRfPf4HroPn0h63DIuh+dSPTUK066HmPP9S+iudEzuRr5jWj6K4bMdAPQ55Pin\nIOVR2Pa4w1rHaoK8KN59PItOL5rJr/0r1P0DbEGQNJCnZtxI8JWvVRAE9Bufh1XPQ9UkxyJ6aiNW\nAiuvnLp69Uznvkv6fkvzFx/H1mCVYzoBkFMF1r0J28cCIisH/UKDKIeAWnfv3zRvEgUtl0CbK+sv\ni1ex7vNzxIbWcNnmybdNITnrMt8eXlx0HaCWtRPze3xOesFlPtjxnkv1UL8GA0nNTeH3U67t7rvX\n7emyrDSYDMEYdO67w7CgUmRoqkBoQgE4elQqBQKc46fcMRhErFYBBBuIer/OnO5t/AC/nlzrsnxI\n4/vdHv/kLjO/ughZpLfDY3tCCOrUj9c3v4xdVNfR3BBRj5bV4txeJ9wUyZH1i2G9Q0KexZGhT68r\nmqn899ZX6f/VIMdCdzHG3jyO+5s8iD3KkyxyAmQVLUjn5kJ+ngCHezr+XWHPnkzuKqbVuv56SEsz\nwkX10Xfx6BYxwVWIXLWcNI5CzW1QEA5Hu4K1yMSqUCAAVAupBnlhsPVxnMGZcqtQN8L9SN6gM/D4\n9TP5dvyrDqsvfQEc78x118VTZXweVcxVmHTrZN7a8pri2DrhdZncbgonMo7z15k/VNce6kc2YFDD\nIW7r4Ctmg5nudXtxLve0y30G3jjIL9cuK5zfcSnR1EeVnCxLFtZOr8CzNeEVAzx9PfntXi9T+/ri\n9Gsw0KWeu0lMUx69+XG3x/dauZ9pa5Q6eaMN5i2Htr/s4YbIevyn2SjV4wUE/tvuFZ9Cp8udeTvW\n7gQzjsPad2DnKPj7eZi5h9c7OFQs+ffcV8IZlaOPKHWtFZelxjjUretecIfK1Pv5+UBaQ9j7ABzs\nKxEIZUl6OnCpPmx+CjY8D2facvFi0b1+6pYJGL5aAYe7O2Yrl+rBn5NYPfh3aoTVpF3N9szvsRgu\nyyT/sS58138FIcaSZ4G+8nzbF12ui3S5/i66XN9Vtexqoaw0HNpMoRKTY8lh8PJ+cEcxXWnkaejy\nKvetWMs3/Zd7pOf1Br1Oz8JeX1Jr8CcOP4WYo5BdDXaO5Mf5jxNpctErFuPZTTB4PyyMg9MRUP8S\njNwFsVlgucnRAb19+3us+ro2F2+YBaHJjgOTmzKtz4sMuNG7qJHuEHKrIW58rtgWEXDMJmwNG5Ez\neiwhnyoDQmYRitqYzFWwuCoyrVp8vI0ff3Rdr8uXIaRY/xkcHBinRDVv+IgIqQALPtWLrMS+xbaI\nVDUXzcB61usNHwyGun9CSAokN4fkJtR6178NaFalOdPv/OjKryuCLDcadjzMwoXPX/U5WPIbfQlt\n/gdVDkBONdg9nCzLSMKM3qm9NKGAMvZPZeHz/fPZcVF98WzzuY0sTfrC7eKhrxj1RscocsNEp8oK\nRKKC3S8ig8M707xkEXXT4eU/1MvBYV7ZPOVFLn4/yRESwhoMKU24+1n3UT69pU4du9swCNlvvcuL\nC2/iMcsHNOIQVvT8RF9eYQqgzIdgt0NQkDIUhDx725Ej7jsoeXyhWrVEUgNgUalmAVtDtgxhMHjQ\nuYo6ON6l2IbAxB9qWKhG+/wXSEqFjOvBGkSw4ep283/uj6fhnmKWfaEp0PVF7l72DT8M+ImwoJLN\nqAvR1EdUvFC+ZcXXB5aWUP6l/yshejfnzX30CcQQdRWCPTqavFGjnb/T0gTHIuz5OEhxhNOWq2FK\nS2pqyR3cF5FjacwBYkglnEzuZhn/crPqvhkZIMtSCyijj37/vesEO2q0bRuYhbKqKiEy4uKk187w\nn/9Z2ZEb41C3+UnNVpZsOPMXC12Yeu9O3snHO9/36nyaUKDyWhulurHi8aS8PLA1bkL6F99wBGkm\n+700J/3bH7HXKMpkZjIpO2yzu2huPuBJSJC8PB0gcIkY8tyGk3Po/tWicgbJDJ0MBvdrCvKZQqAG\nPo4ZjrRu8nwlvtSlbdtKOJX3kC8SP3dbvqQEh1A5mlAAoqKkL/FVrjosMxpGu48j06iE8vLC0qET\nDTnInaxjGItozwZuYg/Wm6Sj73r15L2PqFiALS2edHAhIZ5ac4lUraoeWlqullEbkRc/T0yMdMsv\nvwRm5JOXp9y2f7+0QSaT99ZtBw5U0pGbB5zPPldiuStLPDU0oQBkytTMlSXpzqgWo92WjyyhvDwR\n0fE7d/IFw9hEe9V91EI0l5AOwC/UqOHpCyUgCJCcrPwsz8tyqJSUUE0urByzFf9z7hzIow4fPy79\nrdd7OuoSCcLxwDxJZlRZqRtxg9vyOuF1nb41nqAJBSAvr3K+cH3q92PszeNUy55qNYE761zdJnjF\nUZvdKXXz3j1nuVoxNtY3HUzdup4fFxSkPtq+KIs1FxPjTtAIpKVJt7RoERj1i9pMTD5TKskPRnf+\nHLN4lHQiySeYQ9zIU7xf+RyIPGR4s5GlKpejCYVCoo/AjT87Yr1XIl7vMBXm/eWwsz96F+x4GD7Z\nxIvtXi7vqpWaxx7Lp7h+W6cTVRdxXSHvg1JSfPtc8vM9HxmbTOozVaNsXblnT/edvHwxN7bsA7aq\noqZOy86Wtt/dTFx3/hxRvbvyKHOIwDGFv5EjzOBZwh99uPJM472g1XWtmdD6BdWy22t34ZE4974/\nciq9SerpzFMw/Cmo/2uRa/6Fm9hy7j1urdGufCsXKE51gFMdi20osrWvKKj1FS+9FEzx2YHd7ugs\nfc0sp2a2HBQklphvOtdjP0ABi8Ux65G3R241dfiwex27fE1hy5bAjP9SVGwTLl+W3p/MTNf3K2Ta\n/6E/rR4FN3jZ9+Q9OALLHV1Uy/1FRUgf+nzbF5n2RFdoPdsR+iSnGuz6D1+u7u/MeOcplVooZBZk\ncPfyPtDguLTguj3ct2Igq+5ZR7MqzculbhqlR9lBCWRnl226UZOp5PSS3qgnrVaH4cOlS9LtnTpJ\nfx8+7J0q7PTpwPRsctNZgLAwDw05RBHTd1+7P/83SwMuFCqKdWJMemfSvike4VfEqPd+cFcBZKD/\nWJr0BScyjquW5Vhz+GjHjMBWqJyoUkU6569R49rQ3dasqbQ+chVGQg155xUaqpyOeGKSGhvrucoj\nOFiZCQ/g8GHp7+rV3Z/TIBvuhYcHRu1SvbpyW7Vq0udw880u3q+8PHTZ7jsxIS3wOQ0qisZKrqbz\nlUotFFYf/7mE8lUBqkn5cvmy9DVISfH/0Mho9P+XduSIvB2CqnrDFX37Sk18Jk1STgk8WfvMyvK8\nrWFhkKXSL8ot5LKy3HUAosIfQz5a9xdqC83NFZNtF3Uxm7HVucHt+W2NmvhSrVJRUYSC52tX7qnU\nQsFaQraoksqvFeQdW0nmjmWB2mi4rImMVH7N0dGeH//OO/nExNjR6URuvNFKQoJvN+bkSc+F7KVL\n6qEi5H4J11/vzqJJUAgWuQOZvwhTCbPTpYt0MWbPHtda69xRrs2gRb2e3P+oBzr0J5XN6KlSC4X2\nNTu6Le9Qq5Pbco2rGzU9enKy58d36xZCWpoOu13g8GEDI0Yoe2tPHB2bNvW8V7HbpYHsCpFvO3jQ\n/acrX9wO1EzBMUMpfi1R4Rzqzpkv95Fx5PcbqNhuRU/mjI+x129QNhX1gorizKrTlc0zrtRCYUTz\nh1xG5dQJOsbFPRngGlUeAvGhxcQo1xTURrKuOHNG+nn88YfSisMT1cItt3jup+AqDIfcd8H9wrFI\n7drSLbVrB0YoOKykij9cgc8/l84M3N4zvZ6MTxfSg59ZzIOspgfTeJa7auwl//4H/VDjkqkoQsFu\nL5uKVmrro9jQGnzZ51t6LxjpCBldSH447/d6h0617yi3ugUSnU5qXy5fpPQHpdXTys025Xb8ABkZ\n8o9EKJVHs9o6iCcRdr//Xv2GNuAwjzGLTvxFHsH8wN2Il+9Fr1dKBvkCbv36NjdrPwIXLsB1xbJ1\nerIgXhy93je1yYYNym1//mkAitZjSrTGEgRSb+nK8O1FSYY+fT0HtfwTGmVPpRYKAK1j28IHR6Dh\nz1D1gCPvbeIg7n+6vGsWOORetGaz/0eVckHkLUaj1BRUTSWhthir5i3sCrngUROWLVrY2bVL73af\n665TNrQ7a/iBuwmhSM/Tib8pGDyT66y/kUJdyf7yBfKSnNHkZpQ2m3LW9DYvkE4kS7kPkEods1mU\n3D9PzTLr1VNuCwqSXrtHDwsrVrjP07Fnj/SCX31lpH//8hEKcvXX1YuIt177alRq9VEhRr0RDgxw\nxPjfPQIKKkcu1kLkHWUgkrGUNuWnPKiamrlo69byIbzI9dd7fg252kBtdtOpk/QaagvAcuctMzl8\nyQMSgVBI0ImjfGRXeqDWrCn9XaOGe4kqf6YnTiil1UTe4S3+y2EaEjp5krQesqisngoFNcupKlWk\n7fckDpPc92PjxvIbv7ZqVVFi62vWR2WGXl9RRgKBoax0k+4orfpIns2rWjW1EyrVR+npnl9D7smq\nJnj27JHudOGC8pOSmwoO4WtiuKTYr5DbM1ZRkzOSbXKnsIwM95+uvBOPinLdsekQCZkzk+B5RRni\n5F7HnlqkxceD3OS0Z0/pwQcOeN/t2Gzl9426N/+99tCEApU3IF55UlqTVPno+/x55avcrJlS3eCN\nN7Ncp642gyquOgLIyVHuc9tt0hPdyBG319Vjpx7HJNvkI+d27dyrUuQC7Z57pPvrVfTzIbNnOqW1\nXGh7KsRPnwa5MP7jD+ko/957pY3xZBYSMBXOFZ3mm7zIz/TkbSaSs/dYCQddW5RaKFitViZMmIBe\nr2ft2rWSsoKCAsaPH0/NmjWJjIykbdu2rF692ll+6dIl7rvvPmJjY6lVqxYJCQnkl0dsY42AU9rF\nbPnCqdoMIChILnhEr2Yo8n3VBJnZLB2Bq1mqrFsn7fXOUUO5kwz5PnKT1PBwOyASxSVCVeJUyQP/\nPUnBB/cAACAASURBVPJIAWZdoU5JZChfKI7RnziO7uIFQKne81R9pJZVTR4bqonM/6xWrZLVM57H\njyoFFguhLzvUaL1YTU/WMJF3+CerOaYfvg1ABUqHL3kq1CiVUMjJyaFjx45ckgdqucLEiRPZtm0b\n27dv59KlSwwdOpRBgwZx8Uoc4NGjR5Obm0tiYiLbt28nMTGRiRMnlqZKGhWEu++2UFzN0KiRd3pb\neYet5s2ZkqLcVtYZyNq0kZ5QzWFO3qEt5X5y3GRgOxDbkaO4t8dP/eAbdhHHJWLIIIK1dKMDfzvL\n5TOLtWsN5NoLdVACq+ijOKcoCIhXPOd8Ve/d3MJGf5Yxn1Es5kFG8wnPjZP2D3KrME88cY1G/8/m\nQz6cTtCmjYrtQVgIf3wsujOnVY66elBzevSFUgmFrKwsHn74YebNm4eo8hbdddddzJs3jxo1aqDT\n6Xj44YfJy8vjyJEjXLx4keXLlzN16lSio6OJjY1l8uTJLFiwAFtlcyGshNx4o53iagZlljTvUOvE\n1GIOeWNzLj+nmppRPoJXCwgnnxVdIobxfIhdZWHQHh3NiQnTkDuA1Sg2cTB/8B4Tdo/iZv4FHGsC\n3fiV37iTbqxVveahQ9JPPQOlHs3S8Q7EKIfLt3xW5MknKaRfJqJ3N5ZzN6P4jAdZwieMoe8z8egP\nJDn3a9TIJmlf1arK5yQXrj16+DkfhN1O8ML5LosFi4XgRZ/5tw6lRP6++zobL5VQqF69OgkJCS7L\n+/btS9OmTQHIyMjgrbfeonHjxrRq1YqdO3diMBhoXiwwSqtWrcjMzCQpKcnVKf1CIOzyr2bk007P\n00f6zpw5UvOWX3/17iHIn5max65cl28wqHfanqImeE6dkn5CamGh5WGsAeYxmrtYxwr6cplIznMd\ns3iUE9+sZ5+uJXIHsML1DCE5mdB3p6rWLwgLM3DYUsvXFORWUUFIpxJiUBDZE//r/O3LjCrs+acJ\n3bNNsT0m9wwRI4c6T/r11waKty8xUdkNyb1zIyP9awEkpF9Gf959WkvDgcD2S94iTwIltyDzlIAs\nNPfo0YOoqChWr17N8uXLMZlMpKWlESlTfMZc+XpSvIlaVgaUlAT9WkdubWSz+X+qLnf48lZd8fLL\nxdeeRGbMUDogOCxmik4sit7Fxpc7xCmjrkL9+tJt0dHKhjRpot6hracL/VlBNJepwXnGMROhQT0O\nH1ZWstDE1PTjDwhuYnU3Zz9x7FSYpLZuLR2d16YoZ8E/tCb9qx+wtr3VuU3e9pLum+7CeUwrlrss\nNxw5TNDvvwJw/Lj0ZGoC6NIl6T5Ll/rYw3mIGBKKWMKIwR5Txa91KC2JiSUbPXhCQMbIa9asISsr\ni1mzZtGpUyd2794NoKpyCgQ6nYBOV9TxNW8uKF5Mg8F3eanX6yR/r3buusvOhQtFL9T119u9ar8v\n7b3tNlGSd9ho9O6eZ2friY8vctbJztYrhPvx4wbi46Uj7txcHeHhnl3joYesbN1a1DsOH27DYNBJ\n2jt5spXz5w2kpOgwGETeeCNf0Y6xY20kJXniWCSg0+kYMsTOpk3Ft4tUrarDYABDqLnQ7tMlbbEQ\nFqaTzKZ+/DFIdi+a8CD7yCKMU1zPujtyJJ1B8+YOIVq4KNyokfvnYzh3BuGmm9zWy5iagt2go317\nadRXnU55bnkTY2PFUn2TJWIwU/Dk05hOHnf8lq+GA5bRY/xbh1LSvLnS4MKX+gqiFz3z4sWLneoi\nQRDIKSaKdDodq1evpnv37m7P0aRJE8aMGUPLli3p06cPeXl5CFcUvRcvXiQ2NpZ9+/Y51U7uyMjI\nIDIykvT0dCK8sDUURdF5TQ0NDQ2NIryaKQwbNoxhw4Z5vH+rVq147bXX6Nu3r3ObTqfDaDQSHx+P\n3W5n9+7dxMXFAbB161aio6Np3LixN9XymrS0bMlMoWfPENnCmsi6dT7OvXCMICMizGRk5KqEF7j6\n+PprPXPmmHCMZEWefjqPvn09r7cv7R040ExmZtEoRqcT+eUXz+959+4hEjWXySSyapX0+FdeCeLv\nv6V6kO++y/Y40c6oUcGSsNddu1qYNKlA0t5Fi3TMm1dk9tGggZW5c6Vm1e+9F8SqVZ6lRFy+PJtF\ni4x8+61UXbJ0aTbVqgE2OxEjHkB/7qzq8avpwYu8xdq1ORIz0pdeCmLTJld1UL7vgwYFk56up0kT\nWLIE/vMfOx984MYuVBQJH/UghlPqqTRFo5HLS7+HqChmzjTy/fdF7dPrRdaulV7/qadMkhDb48bl\nM2iQnxebAcPlS4TXq03u0FGYk3ZzjBswDh9E6IiBV31kPPk3ofZco6NVEl7I8Kv6qF27dkyePJlm\nzZpRp04d5s+fz7Fjx+jZsydVqlRhyJAhvPTSSyxcuJDc3Fxef/11EhIS0Pk5KardLmK3F02Q9u8X\nFV6LVmvpO3ObzV4m5/E3jz0WUmxdQWDUqGDOnfM+jZ837T11SuRYMZ8gs9m7ex4UZL0S+sAhyIYO\nLVAcX7OmlZ07i3eEIkajvcQAdoWsXq2TBNBLTtbx3HNF17DZ7EyZYpKE4961S8+sWdJ6zJunl6jK\nXCMSHGzn4EGRnTulJVZrUb3Tn38JU/+BhGZckOzzD615gOe5jEB6ul3iqNexo4VZs1wLJvm927lT\n6v194IBQ4vPJvn84Ef95AEHlBudMmow1LAKsdv79F0n7BEF5/REj8hk6VMeFCzruuMNK9+4FHj+3\nUnHF+qpD0nz27rRgIYhesRYWDsvDZXKgq4T9+0WF/44v/U+pet/FixdjNpsJCQlBEAT69+9PSEgI\nY8eOBeC9996jS5cu3HrrrcTExPDpp5+ybNkyGjZsCMDs2bOJiIigXr16xMXF0a5dO954443SVMkn\n5I46V/mAoMyRr6cEwiJYfk1vrV127ChuwSLw88/KDm/HDqXHlTcB8eRjk+BgZafgiVOX1ep5Z5Kf\nD507yx+AKFkHsTVrTvLfO3hC9zFfM4RFDGMQ33Ibm7hMNEajqPDcvvVWmyItpvt6eG+SWtC1B+nf\nLOf0jbc7t52JbkbGh/8j5+nnnNvkCX/UFNiiWJRe8tIloRyynwlYcMxmig8gr2ZatpQ+JF+tj0o1\nUyhJnWQ2m5k+fTrTp09XLY+IiGDJkiWlqUKZoNdLX9LKJhTKg+RkuQOTd8fLO3dlmGxlXCJwhKqQ\ne/u6wmQSJR2Y2gL1uHEFTJ5cZLUSH6/sPXv0sPLFF565BBsMat67gmKU/NMf0XxsH8fHjFOcw2IR\nyMuTmt/++KOB5GTPx4Bms528vKI6e5rAxdKhE6aNnTh57jLWXCvh9asif7TykBVq39v995uds/d/\n/9Xz6qtB/N//uba68ifx8Vf/bB9QWK25MVJzy9W7lB5A5DFtytrrVUNJaWMfycN7R0crH5rDQU6K\nWg5hTxFFZZ3lHZyaR7Pc69kVJpNj5rFpk3KsJrcqkcd+Ks5NN9kU/hgHDrgWSmq5KOQzBbW2u2L7\ndh39/1OTnsPr8O23yrbI7enV/IRSU6XX89aPpSypV69izBTkZry+ogkFtIB4cjVJIPw2ShuZtnNn\n6dB54EClwvmJJ6TROXU67wLiybOgVami7NwXLpT2qGod+vr1ns0S8vMd/htq9uXyjnv4cIvMoavo\nfl64oCa85HW304gkOvEnL5jeQ7iUJttfurenTn92Owwfbmb3bj2HDul54olgjh6V1ufff6X3Qy0C\nq3z2oJaiNFBUlHwKZfXdakIB9ZFSZeKGG6QdRuPG/p8qyXXx3qrsrFbpq6sm2OPj7bzwQj5BQSJR\nUXaWLMn16jqFOu1Czp1Tfi7yYG9qs0zPvagduvO6dZUnUeskd+/O5rXX8njkkXyK+0BcvKgrUb0m\nAAdoyp/cwWtZzxHV806EKzHJAD7/PEeiMpo0ybPFmJwcSEkpuk82m6BIHSpPDapmVyKfCTZvXn7T\n9717PYwGWM6Uxlu/OJpQADp2lI4yq1SpGCODsuLee6Xtv+8+/5t5yD96bw3OBg0qGl7q9SL9+qnX\n+ZlnCjh9OouDB7O5807vVtCrV5d2RA0aKDsm+YxF7uEM3gk8iwWuu075/qkFOwsJgUcesfDsswXU\nqVN03V69LIoOYv9+6Q0WZY50hmNHCZ36mvN3y5YiO3Zk8dxzDmFw222edcqOAZa0/vKF2nvukQZD\nVAsDLhe2arOfwCCq5Pq+OpHHPvJ1sKsJBZSZlbyN2FnR+esvvdvf/kA+WpQnzSmJe+6x8t13Obz8\nch4rVuTQpYt6hz97tpHmzUO55ZZQ/vnHu9f9s89ynaGxq1Wz88EHytHyiBEWZ6woQRB57DHl6p58\nJuYOoxHCVBL/uRIsR44IZGcLrFqVw6uv5vHee3l8+qmynp6EVQ7+4VvnKvexYwLduoXy7rsO6SIP\nqOcOZXIi6e82bezMnZtHly5Whg4tYP58pf+DfMalFpI7MAisXVsxgqMNH15AcWHbtKlvgztNKAC/\n/irtBPftqxjTxbJCruPdts3/7a9RQ9pJeau3TU0VmDEjiOnTTXzwgUk1Ac7mzTpeftlEcrKOU6d0\n9O/vnWJ66VIjubmOTyQ5Wcfq1crOYe1ag9M2XBQFVq5UDs9Gj7a4SZ9Z1O6bb7ZiNMLp08o8EHJ1\nmyjCsGHB3HZbGHFxoXz5pYHHHrMwfLhFdYSoNvt9iE95jZfIv2J6KeTkoLvsCHP9wgsmibXS2297\nZt9oMsGUKfkIguN6Q4ZYVBfaBw60smBBLtOn56sGDJRHze3Tp/wiJ6tFcb0aWbLESHE1YlKSb8JM\nEwooF7EKX+jKgtwc1Btbfl+JiSmdUHjzzSD+/ttAVpbAmjUGPvpI2Wl98kkQxT8Sm00gLU2xm0vk\nfg5qfg9y09qLF5VD+vBw2LAhm0WLcnjhBbXZRgFPPJHPt986RszKqKEC8hiRv/yiZ+1ao7P8zTdN\nbs1609KUn/oCHv7/9u4/Loo6/wP4a2Z2F1hgcaGMJFErfxRaoJWYnnb90EzNPDP7gb8SwtKvlVrW\nPUArvevCsqtLs04zk/yRXuqVidSZdqUGSupl2A+z/K2phMAu7O7MfP8Ymd2ZWXCX/c2+n49Hj7vF\nBWaG2XnPfD7vz/uN2ZiDe7AeACAkmuSib+qeyIcPe36pyM+3Y+/eOnz9dS0WLNDuL88D+fmx6Ngx\nEddeG4+dO7XHdeNGCzIyeJjNAh580IZp00KTjgoA/fsHY9Wc79TrP1qaRUlBAdIjuCt35Y9bM3UO\nvHo8NxDUqZvepoqqc+7VF2cASEnRLgLzpky6eqxbXYobAEaMcCAx0Tl8NHas+2bGCQnAoEG82zv2\n6moGtbWMnO3kbrhJvd0//6wNHM1VxezXT73tzuP1GW4DADSMGi1PXqiHgLxNBEhNFZtM5fzoIx3W\nrZMC2rlzLKZN006Y7Nihw48/sqiqYvHll7oQzikAH38cGZkoEycqz72sLBo+arHTpy9eyrc1C0Wx\nWnUevLfHPCfHJqfgxcaKGD1aezHu31/7Q73J0Bg71n6hYqyIrCwed9yh/ZB17Spgy5Y6vPaaFR99\nZMFDDzXf4f7kSe3Fbf16PZYuNaBPH+mRddAgZZlrnU7ZZAcABg92KLKDzGYBZnPTvzcnx+7Ss1p5\nXPRwwN7jetQ9Wyh/rUMHZRAxmfz3oVBnQp07pz0m8+cb5JuTw4dZrFwZugvzZZdFxgVBHbi9Sb92\nRUEBobkohpMrrlCe9O4yaPxNnfbaVM+BpgwaxKO01ILXX7diy5Y63HST9vvbttVmY3jaaxgAnn8+\nBkeOsHA4GHzzDadpDNSoQwcRDzzgcLsNjdas0WHIECO2blU/qjg/yUeOsLBYgF69BMyZ04CUFAGX\nXy7ggw+0qbQdOohYscKKXr0cGDDAgY0bmy8maDQC995rR1qagPSUOiSgcTWcgD76cvz+7xKISc7F\nCeonMXV6ri/UZWXczYGo/07qbLVgMRhE5OY2H+jDxfbtyoOmniv0VGRMqwcYx4lwOKJryMhV164C\njhxxnkDXXBP4oDBunB2nTjHYulWHjAxe1TTHM927C+jeveltdTfk4c0wiHrRWUkJh8cf9/z7G+3Z\nw2LKlFh5VbBeL15Y0a3sscCyzvmt/Hw7hg51ICam6YnOW2/lceutnnW0376dwwsvND4mudbrYLEv\n6Q9AvJuZegX/fT7Uw5PuFq/p9cp9dld3KhhsNgZlZRw6dgz/eQX1MWspelIA1TpSZ1vt3Rv47COG\nAWbOtGHTJgtefrkhICtW1a0y7Xbt/Elz1Avi3A1zeOL771lFmQhniQ8GOp0IlhVhMIh45RXnpOzj\nj8ciKysB3bvHY8kS34dOfv216W13t/5iyBD/jE+748lqdvWcSTDOyaacORMZFwh1Bh7VPvJBS6sJ\nthbqu7D4+MgYQ70YbVdXbWG55qgnfHv1atlx6dOHlyej1QwG4MSJWhw9WouHHpI2rrzcOYYuCAwK\nCprPLGp04gSDqir3/zZgAK8o05GVxaNDBwF9+zrwxhvaDKHnnrNhwACHPG8zY4b/hlAOHFBe4N0l\ndqgnukNVdeCSS4QmF0aGG4tFtUDRi3pVrmj4CNLwSUWF84CGss5KKHTqpOxtECkFwC5Gm/+uzfdv\nzocfWnDPPUYcPcoiM5N3e/H0RHq6iI0bLVizRofkZBGffqrD9u06MIyIZ59twIEDLBISRLRvLx13\ndaYNzzc/7yWKwNSpsVi9Wg+dTkRRUQNycpQX8XbtRGzebMHHH+vQtq2IkSMdzT4hJyQAa9ZYL7Rz\njEdysuh1PwNRdP8Unph48eDavr2gSAAJxpCmq7//3YpvvwX69+fdrjAPR7fe6sBXXzkv6RkZLVvb\nQU8KbnhaJri1uOkmvtnXgWCxADk5cbjiigQMGmTEiRP+f0R318fAm5TUlBTgv/+14NChWqxbZ22y\nFMfGjRzGjInFyy83/cjZrZuAwkIbJk+2Y+1aK0pL6/Dll3X46isOAwbE44Yb4vH229LtsLqxPcBo\nqqS62rmTw+rV0vc6HAyefTbGbf+DTz7RYflyA5Yv12uK1PlTSQmHzp0T0L59Al57TXtM1E/m7v4m\n6pz7qqrgDuH06CFg1ChHxAQEQJqHGjbMDr1eRPfuPF5/vWU3MRQUoG0g0tLHrkg1ZYoNN9/sQHKy\ngAEDHHjkkcBnWyxaZEBpqQ42m5TZM3u2m+I+zaitBR57LBZ9+xrx9NMxbicrT53yLpdf7fhxBgMH\nGnHFFQkYMybO7feuWaPDhAlx2LxZj6KiGDzwwMVzXnU6IDNTwLFjLEpKpIu5KDKYPVu6mPfooUxJ\nZVkRl1yi/TkWi5Tzry7fwfPaFN/t2znMmhWLgwdZ7NihwyOPqErA+onDAeTmxqG6moHNJi2q279f\nuX3qTDN3cxrqdSzelkGJRgYDsGRJPY4dq8WWLRZ06NCyY0ZBAZB637oIVfpbqLz5pgHbt+tw7hyL\nbdt0WLw48AO46klbbydx//a3GKxdq8ePP3J4910DFi7U3pHedpvy4mowiG4b5TRl1qwY7NnDwWaT\nVk27S0mdN6+xt7VEm3KqJYrAzz8zKC9XjmU1DhMNGCDgkUfsiIkRER8v4h//qNcMe9XXA/fcY8TE\niXGYOzdWTitmGBGFhQ2aMXj1RPPFViiLotRbesYMKVi7C7ru1NRos4vKypS/6847eTz2mA1ms4hr\nrnE/LKcsgyIGJU3a1ebNHIqKDG4bNbV20bfHbpxStrrVPLq2dt98ozwN9uwJfKbH/ffbkZAgffB1\nOhETJnj3dHLokHKbf/lF+zc7cICF6wXbZmO8KuGhrgHlrhyDwaC8WF0sk43ngdGjY5GdnYB587RB\npvH7585twJEjtTh0qBajRmkH88vLOcXf6ehRFiUldSgvr8Ojj2qPpXqiecSI5o/3kiV6vPRSDL75\nRgpy773n2Y2CtDhQ3XhI+77nnmvA99/XYts2C669VnvBV16MmaDU43JVVBSLl1+OwdChRk2F2dYu\nuva2CeqyxMFoMhNO1CUQtCUR/K97dwFbt9bhrbes+OwzC4YM8W4Wc+hQ50WNYUS333/8uG/BXZ3i\n565u0s03Ky9o6nLbahs2cNi61VmzSInxeGW3unZUYqKI668XkJ7u/txt107Exx9bMHFiA557rh48\nD3TokIA+feLx7bfay4B6yEfd6rEp0mS0cr9aUopenQ4cjHpc7litDLZsia4CmZR9BCA/34Zdu5yH\nYsSIyEhB85fcXDtiYqQ2itnZPO6/Pzj7n54uIj29Zb/LdfhDFJkLr/0bzEwmZRtMd4vILr1U+bV2\n7Zr/mZWV/rnAZGQImDWrHq+8EoP4eBGvv64dYnJlsQCTJ8ehooIDw4jyvNnBgwyeeCIWn32mnDC5\n5RYe77/vfH3DDZ4d28RE4KabHCgrkz5Pl10maBrKA8D69TqsXKnH5ZcLmDVLWyk1Pl69eM2jXx8Q\nV10VWTeJNTVS9lhL119RUAAwfDgPjrOipIRD9+5CUCZaw82YMXaMGRPqrfCcejHTvn3aO9n27X0b\nh87Ls2HWLOlqxDAiJk3Snhdjx9qxerUehw+ziI0VMWNG8wsK7rzTgddec300da5qTk8X5LmA8nIW\nixYZYDQCzzzTgLQ07YVpyhQ7pkzx7Fz95BOdXOVVnUjhbj5n+HAHWNaKI0d0APQYOdLhcUrqqlVW\nvPWWARaLdHzUNZnKy1k88kgsGvf75EkWq1YpV2bfdZdDXs+g14u4/fbQ3KgZDCK6dAld2W5v1NVJ\nGX1ffaVDWpqAlSutXpePASgoyAQBYBjp8T3aayFFgn79HIpm7u6GvC67DHC96Op0old3nJMm2ZGW\nJmL/fhb9+/O4+WZ3v0PE1q112L+fQ3q6oOkTodahgwiTScD581IQu/56HgkJ0vDP889LAeXECQb3\n3WeU6w1VVLD46isv0qbcUA+RcpwInpd+fm6u+6Wvw4Y5oNMJALxLPEhIAKZPb3o5rVSjxxmIduzQ\nPuI884wNnToJOHSIxcCBDmRlhWZBpc3GoKREh8mTw/9GcckSg7xO4dgxqZfIBx94VgbFFQUFSHnV\nubmNKXp6/P47gz//OXT128nFPfaYHfHxwN69LPr25TFypPZOUhqHdl58HA5potmbwDBsmAPDhjX/\nnoQEoHdvz+4md+zg5IAASP1/J0yww2QS5eGpH35gFQXofvyRQ00NvMqcUrvrLgeGDrXj44/1MBql\nkhosK801NLXtVqtU2rq56qstoZ4vaGoeJRhtYT2h7hIYrtQr+N1V5PUEBQVAsQoQ0FYbbO1qa4En\nn4zF7t0cevfmMX9+vVzbP5yNG9f83duNN/K48kpBrqMzYoS2d3GwpacLijF9UQQWL5aykHbu5LB+\nvRXXXCOgTRsRv/8uvad7d96ngABIVUffeaceZ882ID7+4k9MR48yuOceI1JSWFRUSLWIXPtA+yIz\nUxmEOnfW/lxBAN55R4+ff2YxeLADf/hDcIdwevRw4ORJFiNH2jF8eHgEp4vx1wgHZR8BSE1VnpQt\nrUMeqYqKYrBhgx5Hj7L417/0blehRqLERGDTpjoUFdXjzTetWLgwMCkspaUcpk+PwaJFercriV1d\nf72AefMa0Lkzj86deQiC825u+3Yd7Hap5PeHH1oweLAdo0bZsHq190MATUlJ8WwIbeFCg2Iyf/ly\n/90/DhrEY+rUBlxyiYAePXi8+ab27/LCCzH4859jsXixAffdF4evvw7ujdpPP3E4cYLFzp06t61e\nw5F6fUhLC+LRkwK0teNra0O0ISGizvHXllmIXGYzMH584MaDt23jMGZMnHznf+oUi9mzm59sHjvW\njrFj7Th0iEH//vFyw6GuXXno9dJdclGRQV7tfOWVYrNj9JGooMCGgoKm9+nzz51BgOcZfPEF5/EQ\nnT80rlXatYvDqlV6TVezcKTuU9HSQp+t59PvA3UedUvyqiOZ+rEz2jrP+eLLLzlFNs+XX3p+R9up\nk4hly6z44x8dGD7cjhUrpCeCsjJODgiAFCCs/ntY8Ii61EawiySqF7RlZITupIyUxBN16myXLi07\nZvSkACn1cNcuFlu26NCtm9T1ihBPXHedoHrt3d2su0Y56mYpOp13HeP8YfNmHVwn6cvKODz4YPB+\nf1FRPYxGEYcOsRgyxIG77gruuL5UTp5Bjx6821av4WjMGDsqK1mUlurQubOAv/ylZdcxCgqQPnQp\nKSKSkkSYzaImfa+1U5dqbmnWQjQaNsyBoqJ6lJRIH8Rnn/X9hqJXLwHjx9vw7rsG6HQiXnqpwW89\nP+x2z3oTqM+BxknvYElMBF55JXQ3Z++/b8XhwyKuvlqImH4rHAe89FIDXnrJt+NGw0cAFi3So7jY\ngNOnWWzdqsO0adEVFU6eVJ4GJ06E/2nR0AAUFMRgyBAjXnzRENIhr/Hj7Vi1yoo5c/zXQa6oqAHf\nfVeL77+v1fRGaAlRBKZNi0H79gm49tp4t2sDXGVnK594unXz73j+li0c8vNjMXt2TFjO4bVpI+La\nayMnIPgTPSmg8VHZKdiZDqHWowevCAyZmeGfgldUZMDbb0uf2PJyDsnJIvLzI+Mx31NN9WZuiZIS\nHYqLpeN15gyDqVNjsWqVBcnJott1CE88YcOuXRzsdumz4c9ju28fi5ycOLkv+q+/Mnj33RAVNyIa\nPt8SOhwOzJgxAxzHobS0tMn3VVRUQK/X47333pO/VlVVhdGjRyM1NRVpaWnIy8tDgyd9B/1MvXy/\ncaVntLj7bgeck4oihg4N/6CgriEkVUSNHOfPw6MWm/5SU6N8fewYgz59EpCZmYDSUu1NUJs2wKhR\nDvTvLwUDfyZf7N7NyQEBkOYrSPjw6ZNksVjQr18/VDXVGPYCURQxadIkJKpW4OTm5sJqtaKyshK7\nd+9GZWUlZs6c6csmtYh68cxll0VX+k1JieukIoPS0hA1xPWCuh6Nu/LL4UgUgSeeiMHVVyei0+iU\n5QAAC8JJREFUS5cEfPJJcB7W77zTga5dnces8aJstTJyfSdXK1bo8NRTsfjiC+lcaOwZ7Q+9evGK\nSsQ33hgZtYWihU9Boba2FhMnTsSSJUsgNpO3tXDhQpjNZmRmZspfO336NDZs2IAXX3wRZrMZqamp\nKCwsxNKlS8FfbAWQnzX2xm3Utm2E5KD5ibrcsq+F5IJBvZbihx8i40lh2zYOK1ZIwzhWq1ShNBhM\nJmDTJgtWrrRg2jTlI4q7j6767t1dee2Wuu46AcuXWzFihB2PPmrDggU0dBROfPpLt23bFnl5ec2+\n5+TJk5g7dy7eeOMNReDYs2cPdDodMjIy5K/17NkTNTU1OHDggC+b5bUzZ5TDRTU10TV89PTTDRg5\n0o6OHQXcf78dTzwR/gul1KtyQ12+wlPqvgANDS3Pg3/vPT2ysuLRt6/RbQMgtYQEqRvd1Kk2+e48\nNlbEc89px7HUfbr9vU7gttt4vPVWPZ5/vgEJCX790cRHAX92ffLJJ5GXl4fOnTsrvn727FkkqVoy\nJV8oqn5GXdnJz1iWAcs6L/xJSQyyspz/fvnlDHS6lsdLjmMV/xvukpKAf/7TNRB4t92h2N/nn7ej\nvl4qRdCpE4+nnnL49Dfzhi/7O3CgiDFjHPj228YJXBv0eu9/zi+/sFixIgZt20rn8fz5sVi71grW\ngx9lMgEbN9bj2DEWSUki2rQRof6bjx0roE2bBpw5wwLQY8wYHqIYGeezLyLtsxsIAQ0Kn376KcrK\nyrBs2TK3/97ckFMgJSfHg3HpQPHqq+p3sADiff49JlMEVJXzo2Dur9kMbNzY+IoD4KdcUC+0dH9d\nci0AxFz4zztmM7B7t+tXvD9n1b3J1caNc/7/Nm3oXI4WXgWF4uJiebiIYRhYLE3XeLfZbJgyZQre\neOMNGNwk+1566aWorq6GKIryBfrs2bMApGGpQDp3rk7xpHDypDS2+9tvUqOUF16oR69eLX9c5jgW\nJlMczp+3gufDf3zeV7S/wWexAA8+aJSHOrt1c2DBAv+nM4XDvgZTa99fs/niNw5eBYWcnBzk5OR4\n9N6dO3fi4MGDGDdunPxEUF1djYqKCqxbtw6LFy+GIAjYu3evPAFdVlYGs9mMrl27erNZXhMEEYLg\nfEpJSgLi43ns2MGgWzceqak8HA7fn2J4XoDD0fpOrKbQ/gbP//7H4osvnDc2+/dzePllwaPVygAw\nb54Bq1frkZYm4LXX6tGxY/PnO/1to0fAho/69OmDw4cPK7527733YvTo0cjJyUFKSgpGjRqFgoIC\nLFu2DFarFXPmzEFeXh5YTwZG/ejNN/VYt076NJWXSyua33+fMiJI+GqsrNrIbgd43rMSFps3c5g3\nTxqyOnyYxf/9Xyw++ijIFfdI2PLp6ltcXIy4uDgYjUYwDIO7774bRqMR+fn50Ov1aNeuneK/2NhY\nmM1mpKSkAAAWLVoEk8mETp06ITMzE9nZ2Zg7d65fdswbpaXK2LhrFy2mIeGtd28et93mXGQ4fbrN\n4wysI0fYZl+T6ObTk4I3w0kAsGXLFsVrk8mEFStW+LIJftGxo4CyMufraCudTSIPxwHFxVbs2cMi\nPh5eNWi/4w4Hioqcnd3uvbd1lQchvqHaR5BWe37wgR6Nq3oHDAj/Mg8kclksUt8Fs1nEjTf6ktCA\nFiVEdOggYvPmOpSW6tCunYhhw+h8J04UFCDNI7jWjt+/n4aPSGDU1QHDhhnx7bfSOTZtWgOeeSb4\niwU7dWp9BQSJf9BgIoC0NOXdVrt2NHxEAmPbNp0cEABgwQJDxHT2ItGBggKAhx+2IyfHhtRUAbfc\n4qDOayRgkpKUESAxUQQTplVVTp1i8PXXdImINjR8BKnz2vz5DQCiMxiIIrB4sR67d3O46SYeDz9M\nwwqB0rcvj0mTbFi8WA+TScTCheGZ+vzddyyGDzfiyisZ3HknsHcvG9I+ySR4KCgQLFqkx+zZUj7j\nhx/qwbJSNzESGC+80IDZsxuC3nfZG0uX6lFd7XyEWbtWh4wMmpCOBvRseME337B49VUDSkrC+JMa\nIFI/BSd1Jzrif+EcEABoKpf6q80oCX/06YfUfvNPf4qD3S7dGf31r/XIzY2eO2XXO0Ig+E3aSfiZ\nOrUBO3dyEAQpek2cGD2fh2hHTwoANm7UyQEBADZsiK5YmZXFN/uaRB+zWWrKs3FjHYDoazwVzSgo\nQFrR7CotLbo+ADNn2pCRIQWCzEwe06aFf5MdEhyR0ryI+E903RI3Qac6ChwXXUEhNVXE559bYLUC\ncdFbRp4QAnpSAAD89JPyMBw+HJ2HhQICISQ6r34qAwc6wLLOp4NBgyj1jhASnWj4CEC/fjzWrbNi\n61YO3boJGDGCggIhJDpRULigTx8effpQ1g0hJLrR8BEhhBAZBQVCCCEyCgqEEEJkFBQIIYTIKCgQ\nQgiRUVC44PhxBv/+tw7ffUeHhBASvSglFcAPP7AYMsSI6moGHCfi7bfrqZk5ISQq0W0xgJUrnQ1F\neJ7BkiX6EG8RIYSEBgUFAGazsgBemzbRVRCPEEIaUVAAkJdnw6BBDuh0IjIyeMyZE529mgkhhOYU\nIFUHXb7cClEEGGo6RgiJYvSk4IICAiEk2lFQIIQQIqOgQAghREZBgRBCiMznoOBwODBjxgxwHIfS\n0lLFv02YMAF6vR5GoxFGoxFxcXFITk6W/72qqgqjR49Gamoq0tLSkJeXh4YGyvwhhJBQ8SkoWCwW\n9OvXD1VVVU2+p7CwEBaLBRaLBVarFefOnZP/LTc3F1arFZWVldi9ezcqKysxc+ZMXzaJEEKID3wK\nCrW1tZg4cSKWLFkCUfRuwdfp06exYcMGvPjiizCbzUhNTUVhYSGWLl0KnqcOaIQQEgo+BYW2bdsi\nLy+v2ff85z//Qc+ePWEymZCdnY2KigoAwJ49e6DT6ZCRkSG/t2fPnqipqcGBAwd82SxCCCEtFNCJ\n5quuugpdunTBpk2bcPz4cfTr1w933HEHqqqqcPbsWSQlJSne3zjfcObMmUBuFiGEkCYEdEVzQUGB\n4nVRURFWrlyJ9evXIzY21ushJ39hWQYsG7iVahzHKv63taP9bb2iaV+B6Ntfd7wKCsXFxfJwEcMw\nsFgsXv0ylmXRvn17HD9+HL1790Z1dTVEUQRzYSnx2bNnAUjDUoGUnBwv/85AMpniAv47wgntb+sV\nTfsKRN/+uvIqKOTk5CAnJ8fj90+fPh3jx49Hjx49AAB2ux0HDx7EVVddhaysLAiCgL179yIzMxMA\nUFZWBrPZjK5du3qzWV47d64u4E8KJlMczp+3gueFgP2ecEH723pF074CrX9/zeb4i74noMNHhw4d\nwuTJk7Fq1SqYTCYUFhbCYDBg+PDhiIuLw6hRo1BQUIBly5bBarVizpw5yMvLA8sG9tFNEEQIQuCH\nrnhegMPR+k6sptD+tl7RtK9A9O2vK5+uvsXFxYiLi4PRaATDMLj77rthNBqRn58PAHjnnXfQuXNn\n9OrVC6mpqdi3bx8+//xzxMVJj2aLFi2CyWRCp06dkJmZiezsbMydO9f3vSKEENIijBiq2V4/OH/+\nPJKSklBdXQ2TyeTx9/32W00AtwrQ6ViYzfGoqqqLirsN2t/WK5r2FWj9+3vppYkXfU/0TrETQgjR\noKBACCFERkGBEEKIjIICIYQQGQUFQgghMgoKhBBCZBQUCCGEyCgoEEIIkVFQIIQQIqOgQAghREZB\ngcjs9lBvASEk1CgoEJw6xeD2241IS0vE4MFGnD0b+F4ThJDwREGBoKjIgH37OADA7t0cXn3VEOIt\nIoSECgUFgupqptnXhJDoQUGBIDfXDqNRqqCekCDi4YdtId4iQkioBLTzGokM2dk8vviiDpWVLLp3\nF5CWFrEtNgghPqKgQAAA6eki0tP5UG8GISTEorLzGiGEEPciOiiIooiamhokJiaCYWhylBBCfBXR\nQYEQQoh/UfYRIYQQGQUFQgghMgoKhBBCZBQUCCGEyCgoEEIIkVFQIIQQIqOgQAghRPb/VSj2v172\noMkAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 75 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 上記の計算を25回繰り返しプロット図4.9(B)\n", "g = Graphics()\n", "for i in range(25):\n", " sample = genSample(lambda_true, N)\n", " sampleDf = pd.DataFrame({'y' : np.array(sample)})\n", " # k=1でのglm回帰を実行\n", " fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit() # 最大対数尤度\n", " mxLogL = fit.llf\n", " # 切片\n", " beta_1 = fit.params[0]\n", " logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)]\n", " # プロット\n", " g += list_plot([ (beta_1, y) for y in logL200])\n", " g += point((beta_1, mxLogL), size=40, color=\"red\")\n", " g += point((beta_1, mean(logL200)), size=40, color=\"green\")\n", "# 結果を表示\n", "g.show(figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

バイアスbの計算

\n", "\t

\n", "\t\t非常に時間がかかったので試しに10セットのデータに対してバイアスbを計算し、その分布を表示してみました。\n", "\t

\n", "\t

\n", "\t\tバイアスbは、2.7と久保本の1.01とはずれています。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# バイアスbの計算\n", "def bias(lambda_true, N):\n", " sample = genSample(lambda_true, N)\n", " sampleDf = pd.DataFrame({'y' : np.array(sample)})\n", " # k=1でのglm回帰を実行\n", " fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit() \n", " # 最大対数尤度\n", " mxLogL = fit.llf\n", " # 切片\n", " beta_1 = fit.params[0]\n", " logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)]\n", " return (mxLogL - sageobj(mean(logL200)) )" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 時間がかかるので、50個だけサンプルを生成しました\n", "bias_sample = [bias(lambda_true, N) for i in range(50)]" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2.7066865786446273" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mean(bias_sample)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAECCAYAAAD+VKAWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAE11JREFUeJzt3X1wVPW9x/HPObsk2SzZkI0EamsEQsACsTzItdKOeKcz\ntEUvjrWtik/XSxEf2tFBrU6d6fW2FzpI7x2ZYYYRREd5cnzqw7XJ4FWqgwiXiFrBaIBRCQJGTIIb\nNkuS3T33D09it3k6ZLM52bPv1z/RPSe73/yy57zZXdg1LMuyBADIeabbAwAARgaCAACQRBAAADaC\nAACQRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBCQYdu3b9f48eO1ePHiHttee+01zZs3T8XF\nxaqoqNCKFSscX69lWYpEIuKdV4Ch43e648mTrQPuY5qGwuGgmpujSiY5UJ3y6rqtXbtG27Zt0sSJ\nFWpv70y5Dx079okuv/wK/eY3K/XCCzfq3Xff0U9/epVKS8fr6qt/OuB1R6Otmjjx6/roo2MKBosy\n+WN4ilfva5nmhXUbO3bg42RIHyGYpiHDMGSaxlBered5dd0CgQJt3/5XTZgwsce2kyc/0w033Kwb\nb/xX+Xw+zZo1R5deepl2737D0XUbhpHyFc549b6Wabmybo4fIQBna8mSZX1umzlztmbOnJ1y2fHj\nn2jatOmZHgtAHxwHwTQHrqPPZ6Z8hTNeXzfD+PJPV35/3z/f+vXrdOTIx1qy5Gf97tfl79fMyf7Z\npKOjQwcO7M/IdZumodGjC3T69BnHT33MmFGlvLy8jMyTLbx+jHZxHIRwOOj44XkoFBj0QLnMq+uW\nn++XlFBJSbDX7WvXrtWqVStVXV2tysoJjq7T50tI+nLNQqHerzdb1dbW6Z7Vz6motNztUdTa1KAN\nvw1o7ty5bo8yInj1GO3iOAjNzVFHjxBCoYAikZgSiWTaw+UKr69be3tcHR1xtbREe2xbseI/tHXr\nFv35zzWaOnVGr/v0JhqNSZK9Zr4hnddtkUhMRaXlGjO+0u1RJH05j9Pfi1d54Rjt6w9kf89xEJJJ\ny/FDzEQiqXg8OxfNTV5dN8uyZFlWj59t3bq1ev7551RT84rOPffrZ/Wzdx2UXlyzkXbC8eIaD5bX\n14IXleGKjz/+SKtX/647BgDcRxCQMeXlZTIMQ52dnZKk6uoXZRiGjhxp1AsvPKtYrE0LFszv3t+y\nLJ13Xrl27XrTrZGBnEYQkDENDZ/1uW358l9q+fJfDuM0AAbi7b9DBQBwjCAAACQRBACAjSAAACQR\nBACAjSAAACQRBACAjSAAACQRBACAjSAAACQRBACAjfcyQs7r6OjQe+9l5hPKBqO+/gO3R0COIgjI\nee+9t1+//O8XRsQnlElS44e1GjeJTyjD8CMIgDSiPqGstemo2yMgR/EaAgBAEkEAANgIAgBAEkEA\nANgIAgBAEkEAANgIAgBAEkEAANgIAgBAEkEAANgIAjJqx46XNX36ZN1227/12LZz52v6wQ/+WRUV\n39D8+d/W888/48KEALrwXkbImLVr12jbtk2qqJjcY1tjY6Nuuuk6/e53q/WjH/1Ee/a8oZtuulaV\nlVN04YUzXZgWAI8QkDGBQIG2b/+rJkyY2GPb888/o8mTK3XttdcrLy9Pl156mb7//R9q8+YnXZgU\ngEQQkEFLlizT6NFFvW579923deGF30q5rKpqpt55563hGA1ALxw/ZWSahkzT6Hcfn89M+QpnvL5u\nhmHIMAz5/V/9fC0tLfrGN85Luay0NKzm5uaUy/ry92vmZH8n14XeDcUaZzuvH6NdHAchHA7KMPoP\nQpdQKDDogXKZV9ctP98vKaGSkmD3ZaNG+ZSf70+5LBjMl2kaKZf1xedLSPpyzUKhgffvj1fXfaiE\nQgFHv5Nc4PX7iuMgNDdHHT1CCIUCikRiSiSSaQ+XK7y+bu3tcXV0xNXSEu2+rLi4RMePN6Zc9skn\nJxQOn5NyWV+i0Zgk2WvmS2u+SCSW1vd7XSQSc/Q78TIvHKNOou44CMmkpWTScrRvIpFUPJ6di+Ym\nr66bZVmyLCvlZ7vwwpl6+umtKZft27dPs2fPcbQGXQflUKxZth7gw8Wr98vB8PpaePsJMYxYV199\njY4ebdDWrZvU3t6ul1/erlde+V/ddFPPf68AYHjw7xCQMeXlZTIMQ52dnZKk6uoXZRiGjhxp1Dnn\nnKPNm5/Rr351nx544B6dd1651q17TBdc8E2XpwZyF0FAxjQ0fNbv9m9/+xLt2PH6ME0DYCA8ZQQA\nkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQ\nAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQ4LL9+9/V1Vf/iyor\ny1VVNUV33LFUTU1Nbo8F5CSCANckEgldf/1PNHfuP+n99z/Uzp3/p88/P6kHHrjH7dGAnEQQ4JrG\nxk/V2Pipfvzja+X3+zVmTIkuv3yR9u//m9ujATmJIMA1X/vauaqq+paeeuoJRaNRnTx5Ui+++Cct\nWPBDt0cDchJBgGsMw9DGjU+ppuZFVVR8XVVVlUomk3rwwX93ezQgJ/md7miahkzT6Hcfn89M+Qpn\ncnHdOjo69Pbbb+m2236mefO+o8WL1ykWO6M1a/5L1133Iz300H/2+/2xWJskaf/+vykQKExrlkOH\n6tP6fq/z+Uz5/blz3+xNrhyjjoMQDgdlGP0HoUsoFBj0QLksl9attrZODzzyJ533nTt1VNKqbQck\nScb5V0qSHnp8b7/fH+88I0la8dSb8o8qSGuWxg9rNW7S3LSuw8tCoYBKSoJujzEieP0YdRyE5uao\no0cIoVBAkUhMiUQy7eFyRS6uWyQSU1FpucaMrxzU93e2f/kIobisQqPy03uE0Np0NK3v97pIJKaW\nlqjbY7jKC8eok6g7DkIyaSmZtBztm0gkFY9n56K5KZfWLVsPqlyUS/fLgXh9Lbz9hBgAwDGCAACQ\nRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBAA\nADaCAACQRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBAAADaCAACQRBAAADaC\nAACQRBAAADaCAACQJPmd7miahkzT6Hcfn89M+QpnhmPdOjo6dODA/oxd/9k6dKje7RHgkM9nyu/P\n7WM6V85tjoMQDgdlGP0HoUsoFBj0QLksk+tWW1une1Y/p6LS8ozdxtlo/LBW4ybNdXsMOBAKBVRS\nEnR7jBHB6+c2x0Fobo46eoQQCgUUicSUSCTTHi5XDMe6RSIxFZWWa8z4yoxc/9lqbTrq9ghwKBKJ\nqaUl6vYYrvLCuc1J1B0HIZm0lExajvZNJJKKx7Nz0dyUyXXL1jsx3Mfx/BWvr4W3nxADADhGEAAA\nkggCAMBGEAAAkggCAMBGEAAAkggCAMBGEAAAkggCAMBGEAAAkggCAMBGEAAAkggCAMBGEAAAks7i\n7a8B5J5kIq76+g/cHiPF9OlVysvLc3sMTyIIAPoUPXVCG/9yXEV7Trs9iiSptalBDy+XZs2a4/Yo\nnkQQAPRrJH3SHjKL1xAAAJIIAgDARhAAAJIIAgDARhAAAJIIAgDARhAAAJIIAgDARhAAAJIIAgDA\nRhAAAJIIAgDARhAAAJIIAgDA5vjtr03TkGka/e7j85kpX+HMcKwbvxN4hc9nyu8f3vtzrpzbHAch\nHA7KMPoPQpdQKDDogXJZJteN3wm8IhQKqKQk6Npte5njIDQ3Rx09QgiFAopEYkokkmkPlyuGY90i\nkVhGrhcYbpFITC0t0WG9TS+c25xE1HEQkklLyaTlaN9EIql4PDsXzU2ZXLdsvRMD/8jN84vXz23e\nfkIMAOAYQQAASCIIAAAbQQAASCIIAAAbQQAASCIIAAAbQQAASCIIAAAbQQAASCIIAAAbQQAASCII\nAACb43c7BQC3JRNx1dd/MOy329fbX0+fXqW8vLxhnydTCAKArBE9dUIb/3JcRXtOuz2KWpsa9PBy\nadasOW6PMmQIAoCsUlRarjHjK90ew5N4DQEAIIkgAABsBAEAIIkgAABsBAEAIIkgAABsBAEAIIkg\nAABsBAEAIIkgAABsBAEAIIkgAABsBAEAIIkgAABsBAEAIOksPg/BNA2ZptHvPj6fmfI113V0dOjA\ngf0D7meahkaPLtDp02eUTFoZmeXQofqMXC+Qq5KJuA4dqh8x57sZM9L/9DbHQQiHgzKM/oPQJRQK\nDHogL6mtrdM9q59TUWm526Oo8cNajZs01+0xAM+InjqhDf9zXEVvtLo9ilqbGrThtwHNnZveMe44\nCM3NUUePEHr73NFcFYnERsynO7U2HXV7BMBzRsrxLX15vmlpifa5vaQkOOB1OA5CMmk5fjojkUgq\nHicIRBHAcBmK8+7IePILAOA6ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAAkEQQAAA2ggAA\nkEQQAAA2ggAAkEQQAAA2x+922pd4PK7Xd+2SLMnn++qDXhKJzHzQy0DKysZq2rRprtw2AGSztINw\n4sRxrdpYrdFl3xyKedI2Jr5dGx5Z6fYYAJB10g6CJAWLy1R0jvufCiZJBV8cd3sEAMhKvIYAAJBE\nEAAANoIAAJBEEAAANoIAAJBEEAAANoIAAJBEEAAANoIAAJBEEAAANoIAAJBEEAAANoIAAJBEEAAA\nNkdvf21ZlqLRVhmG0WPbmTNtine2q7O9bciHG4x4Z7t2737d7TEkSQcP1utU42HFO8+4PYpam44q\nEW8fEbNI6c/T9X1Nx96Tf1SBq7MMtZE0z0iaRRpZ84ykWU43f6IzZ6arre10n/tEIpaKiop6PY93\nMSzLGvCjzSKRiIqLiwc3KQBgRPjiiy8UCoX63O4oCJZl6ciRE/2WRZJ8PlOhUECRSEyJRPLsp81R\nrNvZa2uLatq0StXVHVJhYdDtcbIG97XB8cK6lZQEB3yE4OgpI8MwFAwWDbif328qFAoqkfApHs/O\nRXMD6zZ4hYVBFRaOdnuMrMF9bXC8sG6h0MDncF5UBgBIcviUETDSdL2uNdBzogCcIwjISpZlqbW1\ndcDnRAE4RxAAAJJ4DQEAYCMIAABJBAEAYCMIAABJBAEAYEs7CPF4XPfee698Pp9eeumllG2WZenB\nBx9URUWFSktLtXDhQn300Ufp3qSnPPnkk/L5fCosLFRhYaECgYAKCwv15ptvuj0aPMY0ze77V9fX\nu+66y+2xRpzt27dr/PjxWrx4cY9tO3bs0MUXX6zi4mJVVVVp69atLkyYOWkFoa2tTd/97nfV0tLS\n6/a1a9fq6aefVk1NjRoaGjR58mRdddVV6dykJ82fP19tbW1qa2tTLBZTW1ubLrroIrfHgscYhqGD\nBw+m3M/WrFnj9lgjyurVq3X33XdrypQpPbZ9+umnuvLKK3XHHXfo5MmTeuSRR7R06VK99dZbLkya\nGWkF4fTp01qyZIk2btyo3v45w/r167V8+XJNmTJFwWBQK1euVF1dnfbu3ZvOzQIYBMuyej1O8ZVA\nIKC9e/eqoqKix7YtW7Zo6tSpuvnmm5WXl6fvfe97WrRokR577DEXJs2MtIJQVlampUuX9rrtzJkz\nqqur06xZs7ovGz16tCorK1VbW5vOzXpOQ0ODFixYoHA4rMmTJ2vLli1ujwSPuv/++3X++ecrHA5r\n2bJlikajbo80ovz85z9XUVHvbwK3b98+zZ49O+Wy2bNne+p8lrEXlVtaWmRZlkpKSlIuD4fD+vzz\nzzN1s1ln7Nixmjp1qn7/+9+rsbFRK1as0C233KJXX33V7dHgMZdccokWLFigw4cPa/fu3dqzZ4/u\nvPNOt8fKGk1NTZ4/nzl6++t08BC1fwsXLtTChQu7//+aa67RH/7wBz3xxBO67LLL3BsMnrNr167u\n/546dapWrVqlRYsWacOGDRo1apSLk2UPr5/PzuoRwubNmxUIBLr/hkJ/wuGwTNNUU1NTyuVNTU0q\nKys7+0k9wskaTpgwQcePHx/myZBrJkyYoEQioc8++8ztUbLC2LFjPX8+O6sg3HDDDYrFYt1/Q6E/\n+fn5mjFjhvbt29d92alTp3T48GFdfPHFg5vWA/5xDR999FE9++yzKfu8//77mjRpkksTwoveeecd\n3XvvvSmX1dXVKT8/X+eee65LU2WXiy66KOV8Jkm1tbWeOp9l9B+m3X777VqzZo3q6+vV2tqq+++/\nX3PmzOnxwkwua29v1y9+8Qvt27dP8Xhc27ZtU01NjW6//Xa3R4OHlJWVaf369Xr44YfV0dGhgwcP\n6te//rWWLVvG24c7dP311+vjjz/W448/rvb2dlVXV6umpkbLli1ze7ShY6Vh06ZNVkFBgRUIBCzT\nNK38/HwrEAhYt956a/c+Dz30kDVu3DgrGAxaV1xxhXXs2LF0btKTVqxYYU2cONEKBALWtGnTrOrq\nardHggft3LnTmjdvnlVUVGSNHTvWuu+++6z29na3xxpRus5nfr/f8vv93f/fZefOndbMmTOtgoIC\n64ILLrD++Mc/ujjt0OPzEAAAkngvIwCAjSAAACQRBACAjSAAACQRBACAjSAAACQRBACAjSAAACQR\nBACAjSAAACQRBACA7f8Ba2ubNw2FjiMAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# biasのヒストグラムを表示する\n", "histogram(bias_sample, figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

Rで計算

\n", "\t

\n", "\t\t久保本のサポートページにあるRの関数を使ってサンプル1000セットでバイアスbを計算してみました。\n", "\t\tこれでも久保本の1.01とはずれており、その分布をみると-20から20の区間でかなりぶれていることがわかります。\n", "\t

\n", "\t

\n", "\t\tバイアスの定義を変形すると、平均対数尤度$E(logL)$は以下の様になります。\n", "$$\n", "\tE(log L) = log L^* - b\n", "$$\t\t\n", "\t\tここで、一定モデルのk=1であることからバイアスbとkを入れ替えるとAICの定義となり、\n", "\t\tAICが予測の「悪さ」を表す指標となっていることがうなずけます。\n", "$$\n", "\tAIC = -2 \\times (log L^* -1)\n", "$$\t\t\n", "\t

\n", "\n", "\n", "Rの関数をdata/bias.txtに定義します。" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting data/bias.txt\n" ] } ], "source": [ "%%writefile data/bias.txt\n", "bias <- function(lambda.true, sample.size){\n", " sample.rpois <- rpois(sample.size, lambda.true)\n", " fit <- glm(sample.rpois~1, family=poisson) \n", " lambda.estimated <- exp(coef(fit))\n", " likelihood.mean <- mean(sapply(1:200, function(i){sum(log(dpois(rpois(N, lambda.true), lambda.estimated)))}))\n", " logLik(fit) - likelihood.mean\n", "}" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[1] 0.7929698" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Rでbiasを計算する関数を定義すると非常に速く計算できる\n", "r('source(\"data/bias.txt\") ')\n", "# サンプリングが1000くらいでもbiasの平均は、結構ぶれる\n", "r('N <- 50')\n", "r('lambda.true <- 8')\n", "r('bias.sampled <- sapply(1:1000, function(i) bias(lambda.true, N))')\n", "r('mean(bias.sampled)')" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC\nAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAmFQTFRF\n////9/f39PT0+Pj49fX1+vr68PDw6+vr9vb277y88Kys7szM+j8//wAA9m5u+/v7/f39/UJC+3Nz\nzMzMPz8/zBERpR0dnZ2drKysLy8v6S4u2FBQXl5ebm5u/0RE/3d3vLy8/8zM/7u7/93dX19fOzs7\nAAAAjY2N29vb3d3dZmZmDw8P0NDvEBAQ7u7ud3d3REREUVFRbGyLMzMzERERqqqq4OD/MjIyBAQE\nioqKmZmZGRkZCwsLR0dlOTk5IiIioaGhVVVVnp69u7u75+fnuLi40NDQoKCgFRUVfHx8iIiIfn5+\nxMTE8/PzFxcXHR0dZ2eFVVV0OTlYcXFxycnJkpKSwsLCc3NzDAwMPj4+9zxE8Wl3QEBeICA/SUlo\nGBg2fn6dtLS0QkJC1dXVbGxsIyNBTk5OPT1cIiJBeXmYlJSUHx8fRUVFLi4u5rPM6KS75MLdXFxc\nDQ0NCAgIfX19Kioqq6vJqqrJPDxbFBQzS0tpMDBOlZW0LCxKFhYWxcXFGxs5LCwsJSUlXl59cHCP\ndXV17IaZZ2dnBwcHZGRkTU1NNjY2SUlJJiYmVFRUa2trYGB/WFh2UFBQkJCQQUFggoKhISFAJSVE\n5ubmPT09Nzc3ISEhg4ODcnKRQEBfW1t5LS1MlZWVWFhY2dn4RERiQ0NhMDAwt7e3YGBgQ0ND1NTz\n2dnZ3Nz7YmJiMTExpKTD29v539/94ODgHh4eSEhmSkppoaHAp6fFeHiXjIyqHx89FBQUAwMDU1Nx\nBgYGNDQ0Dg4OPj5d842Nqamp+HBw+0BA9pCQaWlpra2t5OTkZWVldNvb5QAAAAFiS0dEJloImLUA\nAAAJcEhZcwAAAEgAAABIAEbJaz4AAB1BSURBVHja7d2Nn9vmQQdwdVmaUJauIINlsOAYyeOX6sU9\nyW7r9ur5fGfPTi4HS7ODJCtdoNdmXpZtddmyraXpStuxdUC6sVFgoyswxoAB42Uw3sfYxl+FXvwi\n39mnd+l55N+v+ZwcV7pH0jePbEmPnofjEARBEARBkJXNHW85lAVv+cixt4ZZOlzZx0MtHarstx4L\ntXSYnXY8EPydJw7l5IkwOX5nmKXDlc2FWjpU2XceD7V0mJ32FsADHvBJlw14wAdZGvChygb8isP/\nyF133fWjyZYNeBrg33bq1Km7ky0b8IAPsjTgQ5UNeMADHvCJlA14wAdZGvChyga895w8nLefDJN7\njoVZOlzZk635MQP+x5Mt+9g9oZYOs9PuQI1HjQc84AGfSNmAB3yQpQEfqmzAAx7wqwfP537CgP/J\n3Dj5RMoGfPrwQuGnDPifFuwUxUTKBjwF8LmfMeB/dryBIuBXB965gYAHfJxlAx7wQZYGfKiyAQ94\nwAM+kbIBD/ggSwM+VNmABzzgAZ9I2YCnG35NcGYturIBTzd8rpCbpZCLrmzAUw4vOP6XAPjodgDg\nAy0N+FBlAx7wgF9F+J97xzve8fOAXz345S1wAA94wB8M4AME8IAPsjTgQ5UNeMADHvDOOQAPeMAf\nDOADBPCAD7I04EOVDXjAAx7wzjkAD3jAHwzgAwTwgA+yNOBDlQ34FYd3biDgVxT+9JkcMZMrlU8A\nfmXgK1VyrygIkpAXZVkE/IrAKyqp3bcujVOSC2uAXwV4XtPrjfUpvFTOyfcDPvvwvF5TGk54SSo+\n8CDgsw7fJNVG4wC89NDDPOCzDT92PwC/8UhLAXyW4Xnddj8I/85WG/AZhlcm7gfhN3lSB3x24TVN\nWQLfUHUF8FmFV/VKYxm80lIBn1H4Duk0lsI36qQC+EzCn9HVxhHwjVp1ZeC3trcmk26v19vmjB9b\nmYV/V0s5Er4zrvLZh9/qb/e3xhMDftDb6mcY/qzjQL8QflLlsw/fO8dtC9MJ1zeqfXcnq4f68vld\nh3vjF35x05F3W/DjKr8C8F2D2jkRLvQulDnu0YsX3xPx4LrpDyq890u/fOnUqVOXx/BX3vuYI7/y\nuPWm/cX+fVejK5vOQYXn4S9Y753rZRP+V+Vfe2wO/oqz/j9hwzetc/nswwvb5p/xxDzan9uy4LN4\nqN/PbW444Z980gn/1DV7qjdX4lC/0+/1d7p9a8K93zjGd41X2fxyJxDBC7yqrQS8IT03Mf4pTF6x\nA8/n5pJfMIsBnytKXuArhF8N+KVhB14oODueXjiAJHciT3hP8NYZHeA973wPiQ/evRt67sRQlLzB\nN3XA+9j5HpIqfEkue4RXSB3w3ne+h6QKb1Z4b/CNahXw3ne+h6QJPzArvEf4OlEA73nne0iK8GXd\nrPAe4Y1TecB73vkekiK8qJd9wFfbgPe88z0kPfiyrEo+4OsE8J53voekBy/Kih/4BtkDfHQ7ID34\nsiw2fMG3PwD46HZAevCiXPYH33wY8NHtgNTgjQov+YNXyPXothvwacEbFX4K/8EbN258yBW+cX03\nuu0GfErwZoWfwjupj4D/8Eei227ApwRvVni/8A8SPrLtBnw68FaF9wu/+fQosu0GfDrwI+sqvV/4\nh/Yj227ApwNv3ZbzDf+MHNl2Az4VeLvC+4Y/Q4SothvwacDbn/D+4TcKorey3QP4NOBFu8L7hxcL\nUW034FOAn1R4//ACKUe03YBPAX5S4f3DSyQf0XYDPnn4aYUPAL8veirbPYBPHl4clgPDjwqeynYP\n4BOHL8slKTD8muNDHvBhNiF5eHEoBYeX5Hw02w34pOH5WYUPAu/4kAd8mE1IHL44q/BB4B0f8oAP\nswlJw/MkfxjeWwscC97xIQ/4MJuQNHwxZ1r/ut3DzUetn9c+5gPe8SEP+DCbkDD8GhFM+JsbVj5u\n/fzEU37gc6KHst0D+GThc1aFl246ca/5gp9drgd8mE1IFl4ga2HhZ5frAR9mE5KFLxSlsPDS9J48\n4MNsQqLwJcKHh59+yAM+zCYkCj+c3J0JAy/m3Mt2D+AThB9NbseGgs+TKLYb8MnBz27HhoIvk7UI\nthvwycE77s4Eg//ks+tmnjsjWvmNMCsO+MTgHbdjA8Jfe96Cv/WC5Z7bC7PigE8M3nl3JiD8U9ak\n2bJ/416YFQd8UvCnnXdnQsHzdu/1gGcD/lZOigi+odcB77LzPSQh+NP23ZlI4Gsq4F12vockBH9r\nV4oMXq0B3mXne0gy8BVy3xHwHnvEmMB3COBddr6HJANf3V0/At75Nw/wCukA/gQT8EaFjxC+oQ0A\nf4IJ+GptPUr4ahvwJ1iAr5BOpPDWJRzA0w9frTUihbcu4QCeevgKqUcL3yB1wDMAr7YaEcObl3AA\nTzu8ojejhjcv4QCedviB3ogavk4ATz98S40cXiH8qsEvGKSW0kGF33fVXuM6UQyqT1150ZE5+Mef\ncP5t4aDCh2ZsNbm93wyz4nQOKpypGl+rWhXZbjY1Ttga32hXV63GMwdfIfxhz9DwqgZ4yuFVrRED\nfIcAnm5461wuevgG6QCeavimrsQCrw0ATzW8Vm24w1++dOnSS/7gq23A0ww//mrnAu+z6ZWZQQvw\nNMOPv9pFD8+Tl/fCrDjg44Uff7WLHr5B7tsLs+KAjxXevmoXC3xN2wuz4oCPFb5dbcQFr+7uhVlx\nwMcJr5B6bPDNV/bCrDjg44QftBqxwVfIp8OsOODjhNfU+OAbr/xWmBUHfIzw05P4WOB3PxNmxQEf\nI/z0JD4WeG0YZsUBHyN8axAn/GeJ99U8HMDHB38/qcQJ//KkE6RAAXx88LuOI30M8Ov3lkKsOODj\ng396EC/8q6L39TwUwMcG/+D0cm1M8J8peF/PQwF8bPAfeFcjXvjPhfl2B/jY4B9+1TP8S5cvX77h\nG/6358aW9hnAxwWfJ2c8wy+f8Sj4Fwuj4CsO+Ljgi48cMbhUNPDFYvAVB3xc8PLLscOPQny7A3xM\n8AJ5Z+zwc4OK+wzgY4IX9zdjh5dCfLsDfEzww1EC8IXgl3AAHw/8GuETgBf3A6844OOBFwtSAvCl\n4HdmAR8P/HCUBLxxXAm64oCPBd4QSQJeIvmgKw74WOCNI30i8LnA3+4AHwt8QUwGfjoGne8APg54\n3hxCNgn4vBx0xQEfB/zIHHcoCXg+cPMrwMcBbx7pE4GXhqWAKw74GODL1vgzicDvB/12B/gY4Euy\n5A8+WEMME14MeoMO8DHA7xd9wgdremXCC0GbXwE+Bnh7aMFE4APfoAN89PB5IiUHH7T5FeCjhy/u\nJwhfDHiDDvDRw49HjU4GPugNOsBHDm/eoEkOPugNOsBHDm/eoEkOPugNOsBHDl8YJQof8AYd4KOG\nt27QJAH/O7+7aeb2dcGKXzrARwxv3aBJAv7J926YeYbkzMg+5QEfNfy+mBS8PYaJOTyNkRzgfSZi\n+LJ92S45+EZrAPggiRg+L0sJw1ergA+SiOHHl+0ShB+0AB8kEcOPL9slCM8TBfABEi385LJdgvDW\nCMOA951o4SeX7RJoiDGFN0cYBrzvRAtfEIPAL5/RC7w5wnBW4Le2t6aTXq+3NXmDdvjpZbsk4Tsk\nM/Bb/e3+1niy1Tfgx29QD1+answlCK+QTlbge+e4bWE86fa6O5M3qIe3W9slDN/QBpmB73Ld3ngi\nXOhdKI/fuKhprwX6ZpBUpmMPGfnox51Mn//8UvgvPLF8xie+4GHGapvjap20t91zji//X05446/n\nepNXlNd4gZTTqPFNPSs1Xtg2/9iTc1sG/PgN2uHFnJQGfIVUMgK/0+/1d7r98aTX37Je0Q8/aYOR\nMHxDb2YEnuO6s8lO1/EG1fDOk7lE4dvVzMAvDc3wpaGUDvygBXjfiRDeeTKXKDxPFMD7TYTwcj4l\n+IZeB7zfRAc/dzKXLHxNBbzfRAfvuDOXNLyqAd5vooMviL9305kvJgffIWcA7zORwRsnc549o4Zv\nkGcA7zORwZdkKSh8uIYYZmq7gPeZyOCLxcDw4ZpemVF1wPtMZPByKUX4jt+eMQAfFfyacTKXHnyD\nnPVJB/iI4M2TuRTh9ds+6SiCV5iGz4mpwu9e90lHEXxVG7jb0wpvdWqYIvx1n/2e0QTP8Wqr3WQU\n3urqKkX4DZ/f7qiC57i6Rloqk/Difsrw1/31jEETfKeqt5uKcnSlpxV+OEoZ/ra/rutpgleb1kd8\np8IgvN34Jk34s/4+5KmCN380eSY/4+3GN2nCb/rr/ooeeL6m14zoFSbh7cY3qcL76/6KHvhxjXcL\npfD2Y/GpwvvrwZwmeE+hE37NbnyTKry/IYbpga/xqnmor7EF/yW7k7nbLavPuVThJTnvh44a+I7C\nd8ywBb83tHqZ03etPufShd8v+qGjBt68cMe1SZUxeNFaK1Jf5OkZ/vKlS5deCg0/8tORNU3wWofX\nOI2t0zkbvk4WenqGXz6jH3hfHVnTBE+4pspVOwzCqzUa4H2NRUYTfFXV+HpLYRBeG1AB72e0Cprg\nlUGda7L25c6ErxCeCviSjwFHqYJn8Vu9Cd/U/XrGA+9nwFGa4Ntt1QiD8NU2HfB+hqSiCZ5wHkIj\nfKvp1zMmeB+jitMEX+XZhOdJhRL4vPertlTB6+xdsjXhBy3fnjHB+xiZiCZ4Ji/ZGvDtKjXw3geX\npgmea6rNuls7Wwrh9To18N6v2tIEr9ZUdaAxB88ThRp471dtaYLXlQ6Ll2xVzb9nXPDS0OsJHU3w\nLRO+xtxNmppKEbznq7Y0wTe1aq1dY+5QT6Yf8RTA571etaUJnuNV1e1BGvrgJ7dk6YAve32ghh54\ndRzW4Ce3ZMPAf/DGjRsfigTe81CzNMG3W6paY64FzuSWbBj4aJpeWRl5bGtLDzzHaeY5vMZWu/q9\nlye3ZCmB93pCRxO8bsK3GTudu6ov90wD3mszHJrg1Vq9o7J2AedWe7lnKvAeT+hoguea7Zpr1wi0\nwd/bXO6ZCrzHEzqq4L2EMvhPT2/J0gLv8YQO8OHgf//eIzxTgfd4Qgf4cPCv3qIE/mPXNsf58HlB\ncP9mD/hw8E9/lhL4pz6xMc795PWCexMswIeCF8jLtMDPZmwNBMC7JxS8eH6dPvhqG/AeEgq+8BkK\n4esE8B4SBr5MPkchvEL2AO+eMPB5+UUK4RttD12fAT4MfHGfSvjB+QzCnzyct58Mk3uOBV/2D/7w\nj15c6vn43NjPV64cMaMTfn7GJx5f/hu/vGzGCnnEddWP3RNmpx0LsdNO3sF8jecJT2WNb5xvuddZ\n1mo8TfCjoRQN/PIZg8Hf1t3pAB98G/aLlMLvuT9DB/gQ8CRPKbyH7k0BHxzeHFKUUvh3ufZ8Bvjg\n8GJOohX+Vddn6AAfHL4wohb+dddOUQAfGN7qop5S+A3XTlEAHxi+JEv0wrt2igL4wPBWF/W0wrt2\nZQ34wPBynmJ4166sAR8UXrC6qKcW3q0ra8AHhTdP5iiGL7mc0AE+KLw13hi98G69XAI+ILw93hi9\n8G69XAI+IPzIGm+MYniXEzrAB4TfFymHdzmhA3ww+DLJRwgfQY8Yh+BdTugAHww+L0sRwkfdAueT\nz66vr79wSxxn0dc8wAeDL+7TDH/teQP+6r1j94VPUQI+GLw9siS18OaMFTLuMlAE/KIEgrcv29EN\n39AGgI8a3r5sRzm8WgN81PD2ZTvK4etEAXy08GuEZwC+QeqAjxZeLEgswLergI8Wfnqkpxt+0AJ8\ndPBlUTxNTq+Pc+UrFMNXSAXwkcELQ/HMcxP39T8Odgk+GfhGawD46OBznDYbmuAK1fDVNuAjhK84\n+i2mG75OAB8h/GycOdrhFdIBfHTwjiM95fCNmgr4yODvd3ZfSzm81Rk44KOB39Ucu5hyeN48oQN8\nNPBPO8YiiQh++Ywh4Rt6E/ARwT8411E57fDVKuAjgr/9Eef+ph2+qQM+IviHX2UJ3myGA/go4PPk\nDEvwZjMcwEcBX3xkgyl4tQb4KODLZI8t+DoBfBTwJXmTLfgGqQM+AviCyBp8WwV8eHierLEGP2gB\nPjy8WJBYg+fJacCHhpdLzME3WlcBHxa+RMrswVdfAHxY+FxRYg+++QrgA8KvCXbeIGc3N6/N7W8G\n4BXyVcAHgyc5O7tPb2xsfGJu7zMA33hOA3xA+HHZevPQTmUBXtsFfCj4pq7EBB9bCxwr9y3qDQfw\n3uG16uGdygL8+qLxKgDvGZ63m94wCP/Cgu5NAe8ZvlpbsFOZgD+zoHtTwHuFr5A6q/CnCc8K/Nb2\n1nSytV3muF6vt5UuvNpatFOZgF8fHu7elE74rf52f2s86V7o9s9t9dOGV6xzOUbhxX1G4HvnuG1h\nPDHqe7dn/LeT7qHePpdjFD5PWIHvGtjTyc7gnHChd8E44D968eJ7Dg9SG/ugwiZ8a/LA3PzIvl8+\naqzguX8hCQ4qfGDGF1/8EnkfG4MKz8Fv97vme+d6acI3icIw/JtX99iAF7bNP/bk3PtN9S0LPr1D\nvX3xxgyLh3ppVGDjUL/T7/V3un1r8v4L/f5213iV5pe7zuy5KSbh1w6d0NEJz3HduYnxT2HyKh34\nWnXJTmUEXhqWGIFfmlTgHRWeUfjiwRM6wHuBr7WX7VRW4Esy4APAkw7r8GUiAN4/fG3pTmUFXioc\naHgHeHd4wdG9WTzwUY9JswD+4Akd4N3hc2S5JwtNryz4gyd0gHeFF0gW4A+e0AHeFT5XPIKJHfgD\nJ3SAd4PPEz4T8Afu0AHeDX5YlDIBX55vcgl4F/iSXM4GvLQ/d0IH+KPhy7IoZQR+NNfkEvBHw4tG\nhc8I/PyI8oA/Ep63RpTMBvz8iPKAPxK+aA03lRF40XnxDvBHwQtEyBD83MU7wC+EH5FxblphG/5P\nbo5jb9Ia4JfDi6JZVnPc/oJx+OlvtMYmygmAd4FXJk2qswJvDTYLeFd4dfIMRVbgrWELAO8GXyFN\ndya24M3BZgHvBl+btruJHT7eFjiz32gOWwB4F/jOrN1NZuAVUge8G3xrNrpcZuDNYz3gj4affrPL\nFLxxrAf8kfCnJ4/DZwveONYD/kj4F5wtqrMDbxzrAX8U/H1zLaozBN/UAX8U/L3OYUSzBK+QZwC/\nHF585eWMwjfaDwF+KTwvX13PKnydvBPwy+D3C+uZhW888DLgl8ALRMgw/EOPAH4J/LAoZRj+ut0O\nB/Bj+DXiyM2bf5pZ+I3zI8A74IXcuABFV33cRGUQ/nYB8IvgrT5rswz/htXqDvAH4O1OqrMMv5kr\nAv4wvN38ImH4eHvEOAhfksuAPwg/bn6RMPzyGeOAl+QS4A/Cj7suzTZ8sQD4A/BNvbIC8Objk4B3\nwk9b0mcbXjK+3gHeCT9tb5VxeOPrHeAd8Io+8MvEJrw0HAHeAT8Zbyj78OIQ8DP4yqyBZdbhy/Jr\ngJ/CV1tL9lX24KXinwF+Aj97Vm4F4HnyNcCP4avasn2VQXjpz78OeBv+dWen9NmH7x7swH5l4R+q\nLd1XWYT/i0/nAG/mDWeFXwX4b4Sp8hmC/4AeiIld+L8shqjy2YHnyfVVg+dDVPnswBevb6QIf/nS\npUsvJQ0vFYeBd1pm4HlyNk34ZJteTeDLZoMMZuAXDFIbwaDCf/XX33xsztPrWMEsDCo8/xsf++ab\nZv7mb998c0/+EmWDCidd43lS2ly9Gi9JQzHgTsvKoV4cSisJL8z1ZL568MaH3WrCS/uFYDstI/BG\nhV9ReF4OdrDPBrxZ4VcUXhoFO9hnA96s8KsKL+UCHewzAf935nhDKwsf7GDPNvx+zszXv/XAmY2N\njXc/uyLwz757w8zf/4M1WZfyQa7csg1PBDNfO39708ingl1ZZw/+qU+Zm7v5j/9k/vz2hiQVzUfp\nVgve+oVNuy19wFsqDMI7Z9w04MsF/7fpsgA/fnhmdeGlNf8f8xmAH1f4VYY3PuZLqwffqvpmyhy8\nJMo+z+bZh58MNLXa8L6/4LEPP34ePmX4ly5fvnwjTfhyoeBLnnn4zqTCpwvvkSk+eKks768UfK3q\nZV9lF/6fn9oY5/4HHso980zOGX75rmMdflbhVxT+2r9sTnKWfPjWLcGRwhFX9FiHn1X4VYV3zNgk\nQ825p3PZhXdUeMCb8i+sCLyjwgPeyHf+lV8NeEeFB7w543d0fiXgqx731crAP6lpSvbhS0TxuK9W\nB15xyGcVvjwkETBlDL6h6O2sw4sy4BfMyOvVbMOX5VEUTJmDb/BkkGn44lAC/MIZm6SZYXiBCIBf\nMqNqn9RlE76wLwF+2YxVvZJV+JHMA375jNZJXRbheXkk0QSfZI8YXmZUtFo24XMFiSr41JpeLZvR\nPKnLIPxIXgP80TN2SJM5eNk5ZiRZ0KJojYwkwLvM2CTXWYMnztmFw0+JlM1v9IB3m7H6wIMZg98f\nlgHvYUa9tbzhLYvwRfsDHvBuM55pLX90nkH4EREkwHuZceMNuZgd+BEpSYD3Br+5/GlK5uCLDnfA\nu3adUVr2NCVj8OWcLEiA9w5vyC8+p2MLviQXeAnwvjrLKS5+jpYl+NKQiN++ORfAu/eStFieGfg1\ncUiKvLS+HgXTSsEvlmcA/nWhJOZkUhiZV20A7x9eKiyQpxmer6vtmk4IGebE/PhaHeADwJcXyNMK\nzw8Mcr2mqvWzrzu/zgE+ALxU3j8kTye82iKaWrcfj5o++U81fGo9YniCNz/n89TDlwpEG8yeiWMD\nPhhTcvCSSEZUw5dFWS4S59oDPhJ4qSTvO+/V0QVvsA/F8vx1GcBHAy+tFWTHRTyq4PPysHToghzg\nA8N/ce5i17+VRbI/7RWHKnixuOBKLOADwx/qHksoyGKZRngR8LHCG5/0Q7m4BvjVgzfoC6QwWqMF\nXhDNaNq6lbmd+sln1x158iuADzbjbD+e1p4j//6t//jsV8VpnN/3y+Jc1mKFF3NL4a8974R/NnKm\nVYGf24+n9f/8FiHkld3dFzQj+tlxz3imsjB0uudEn/Bb21vTyezHUnjR/D+qGtfeB/zBGf+r0VA6\ndVVVa0b083NPMhDNeE+1tUSf8Fv97f7WeDL7AXia4GdxnPALwtnrnabxD4IEg++d47aF8WT2A/DU\nw09H8woK3+W6vfFk9oPjLmraawtmHxCE0gwWcB33D8/FNX584IQrmwu1dKiy7zweaumYTueEbfOP\nPZn9AHyEZdMJv9Pv9Xe6fWsy+wH4CMumE57juo5J1/FGNuH/++67774r2bJphV+aTMK/7dSpU3cn\nWzbgAR9kacCHKhvwgAc84BMpG/CAD7I04EOVDXjAAx7wiZQNeMAHWRrwocoGPOABD/hEygY84IMs\nDfhQZQMe8IB3yVvuPJSTd4bJd78RZulwZb82nv7P9773vf9NtuxvfDfM0js7IRa+IxB85Pn+o+mV\nraVX9KPfT69sOgL4Fc0Pfphe2RfTK/qHP0ivbARJMV3nI5nJJqVi091oWtK90OVmj2QmmpSKTXej\nqUnv/7rc7JHMZItOp9h0N5qe9Lvc7Mm8RJNSseluNAXZ7hkB/Oql3DVi7YPZI5mJJqVizaS30fTE\n2AezRzITTUrFprvRlKUb/lcwVCwVpSMIgiBI5qJafcB1at6XcM6rqmmvPxIwNp3Ce18C8JmIOTiW\nyvEqp+rGVKkSzfo3YL+w3uOrNaKqek1pquasJnyzZb5Qda0NeFaj6hWlPejUOhrHafxA5TqWpfXC\nfq+jKwoZcO2mMSvXHhjwvKYotWazpigtwLMa1arDRi1WOs1Wp6417YO+/cJ6zzy0E3NGe1bjj1rt\ndFS1XcehnuFM4PmW2tQ6xjFfa1vvmy/s9xbA14zXzVoH8AxH1Uw+A9Mg1DrNuoXc4awX9nszeGtW\nA75p/NsY1I1PA64GeFajalpbU4w6r6s1rdppqVrVdLZejN+rTeG1dksxK32tVjU+5o0FNcCzm0pn\nMlHMP+MTO+uF/d4kqlqZ/IW3XnQqaa87kkTwgb6i6XTC/w4EQRAEQRAEiT3/D29I+WI9k+lIAAAA\nJXRFWHRkYXRlOmNyZWF0ZQAyMDE2LTA3LTE3VDA0OjMyOjQ4KzAwOjAw0uIhBQAAACV0RVh0ZGF0\nZTptb2RpZnkAMjAxNi0wNy0xN1QwNDozMjo0OCswMDowMKO/mbkAAAAgdEVYdHBkZjpIaVJlc0Jv\ndW5kaW5nQm94ADUwNHg1MDQrMCswpXe8owAAABR0RVh0cGRmOlZlcnNpb24AUERGLTEuNCAcRzp4\nAAAASnRFWHRzaWduYXR1cmUAZTE0Nzc3MjYwOTExNzdlMGIwNzIwMTc3OTZhNWMyMjEwYWY0ZmM2\nNDQwNjE3OTc2MmZjNThiODBiYTdjMWI0ZEc/FykAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "graph = preGraph(\"images/fig-4.10-R.pdf\")\n", "r('p <- qplot(bias.sampled, geom = \"blank\") + geom_histogram(aes(y = ..density..), colour = \"black\", fill = \"white\") + geom_density(alpha = 0.2, fill = \"#6666FF\") + geom_vline(aes(xintercept = mean(bias.sampled)), color = \"red\", linetype = \"dashed\", size = 2) ')\n", "r('plot(p)')\n", "postGraph(graph)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 7.2", "language": "", "name": "sagemath" }, "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 }