{ "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\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": false }, "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": { "collapsed": true }, "source": [ "\n", "\t

例題のデータ

\n", "\t

\n", "\t\t３章の例題に使われているデータは、架空の植物の第i番目の個体のサイズ$x_i$と肥料の有無$f_i$が種子数$y_i$\n", "\t\tにどう影響するかを調べます。\n", "\t\t久保本の図３．１を以下に引用します。\n", "\t

\n", "\t

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

\n", "\t

\n", "\t\t３章のデータは、ネット公開されており、以下の例でもネットのデータを利用します。\t\n", "\t

\n", "\t

\n", "\t\tデータをロードしたら、summary関数とvar関数を使ってデータの素性を確認します。\n", "\t\tここで、確認することは、以下の3項目です。\n", "\t\t

データの個数と欠損値の有無
データの平均
データの分散

\n", "\t

\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 3章のデータをネットから取り込む\n", "d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

カテゴリ

\n", "\t

\n", "\t\t肥料の有無（施肥$f_i$）は、CとTという文字列で表されており、Rでは因子（factor）と呼んでいますが、\n", "\t\tアンケートなどのカテゴリと言った方分かりやすいかも知れません。\n", "\t

\n", "\t

\n", "\t\tPandasのデータフレームもRのData.Frameに劣らず、$y_i$は整数型、$x_i$は実数型、$f_i$は因子型と扱ってくれます。\n", "\t

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

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

	y	x	f
0	6	8.31	C
1	6	9.44	C
2	6	9.50	C
3	12	9.07	C
4	10	10.16	C

\n", "

" ], "text/plain": [ " y x f\n", "0 6 8.31 C\n", "1 6 9.44 C\n", "2 6 9.50 C\n", "3 12 9.07 C\n", "4 10 10.16 C" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# どんなデータか調べる\n", "d.head()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "

	y	x
count	100.000000	100.000000
mean	7.830000	10.089100
std	2.624881	1.008049
min	2.000000	7.190000
25%	6.000000	9.427500
50%	8.000000	10.155000
75%	10.000000	10.685000
max	15.000000	12.400000

\n", "

" ], "text/plain": [ " y x\n", "count 100.000000 100.000000\n", "mean 7.830000 10.089100\n", "std 2.624881 1.008049\n", "min 2.000000 7.190000\n", "25% 6.000000 9.427500\n", "50% 8.000000 10.155000\n", "75% 10.000000 10.685000\n", "max 15.000000 12.400000" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 個数と平均、最大、最小を、標準偏差をみる\n", "d.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

データの可視化

\n", "\t

\n", "\t\t施肥の有無（C, T）別に種子種$y_i$とサイズ$x_i$の関係を散布図で表示してみます。\n", "\t\tseabornを使うと、このようなグラフも簡単に表示することができます。（施肥別の指定は、hue='f'で行っています）\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAFhCAYAAACBP4ZvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wJHd95/H3zErr9Wple70rY2kXr+1Afr4YTLBDWGwS\njCFxHOJzjicn5xCeKxByqVxBUiSVmABVuRwXU4TkgLscJoRLcubC07kSnhwID147gA3YSeyf1w/Y\n1kq2tWJ3LWsfpNmZ+2M0WkkrzW8kTU/3SO9XlWunW/3w7e6xPtPdo/6WarUakiQ1U867AElS8RkW\nkqQkw0KSlGRYSJKSDAtJUpJhIUlK6sl6BSGEZwGfBd4fY/xQCKEH+DjwDOBJ4JUxxkNZ1yFJWrlM\nzyxCCJuBDwK3zBn9ZuCJGOPzgZuAn8qyBknS6mV9ZnEUuAp455xxVwPXA8QY/1fG65cktUGmYRFj\nrALHQghzR58L/HwI4b8Bo8CvxxgPZlmHJGl18rjBXQLuiTG+GPhX4PdyqEGStAyZ3+BexGPA12de\nfxH4w2YT12q1WqlUyromScpS1/8SyyMsPk/9PsZfApcAsdnEpVKJsbGJDpTVHgMD/dabkW6qFaw3\nS91UK9Tr7XaZhkUI4WLgBmAXMB1CeCXwH4EPhhDeCEwAr82yBknS6mV9g/tO4MWL/OjVWa5XktRe\n/gW3JCnJsJAkJRkWkqQkw0KSlGRYSJKSDAtJUpJhIUlKMiwkSUmGhSQpybCQJCUZFpKkJMNCkpRk\nWEiSkvLoZyF1jWqtyu2jdzAyOcpQ3yC7By+hXPIzltYfw0Jq4vbRO/jGvj0APHDwIQAuHXpeniVJ\nufAjktTEyORo02FpvTAspCaG+gabDkvrhZehpCZ2D14CMO+ehbQeGRZSE+VS2XsUEl6GkiS1wLCQ\nJCUZFpKkJMNCkpRkWEiSkgwLSVKSYSFJSjIsJElJhoUkKcmwkCQlZR4WIYRnhRDuDyH8+oLxV4YQ\nqlmvX5K0epmGRQhhM/BB4JYF408B3gmMZLl+SVJ7ZH1mcRS4CljYBOD3gD8HpjJev1QY1VqVPSPf\n5u/2/j/2jHybas0Ta3WPTMMixliNMR6bOy6E8KPARTHGTwGlLNcvFUmj694DBx/iG/v2cPvoHXmX\nJLUsj0eUvx/4T8uZYWCgP6NSsmG92emmWmF+vQeGx+np2XBiuDpeuO0pWj3NdFOta0FHwyKEMAQE\n4K9DCCVgMITw1Rjji5vNNzY20ZH62mFgoN96M9JNtcLJ9W4tb6NSifOGi7Q93bR/u6lWWBvB1smw\nKMUYR4BnNkaEEB5KBYW0Vth1T90s07AIIVwM3ADsAqZDCK8AXh5jPDgzSS3L9UtFYtc9dbNMwyLG\neCew5JlDjPH8LNcvSWoP/4JbkpRkWEiSkgwLSVKSYSFJSjIsJElJhoUkKcmwkCQlGRaSpCTDQpKU\nZFhIkpIMC0lSkmEhSUrKo/mR1DHVWpXbR++Y91jwcsnPSN3AY1cshoXWtEYrU4AHDj4E4GPCu4TH\nrliMaa1pI5OjTYdVXB67YjEstKYN9Q02HVZxeeyKxctQWtNsZdq9PHbFYlhoTbOVaffy2BWLl6Ek\nSUmGhSQpybCQJCUZFpKkJMNCkpRkWEiSkgwLSVKSYSFJSjIsJElJhoUkKcmwkCQlZf5sqBDCs4DP\nAu+PMX4ohPB04EagF5gCfiXG+ETWdUiSVi7TM4sQwmbgg8Atc0a/F/hIjPFy6iHy9ixrkCStXtZn\nFkeBq4B3zhn31pnxAGPAczOuQWtMtVbj1rtGGR6bZOdAH5ddNEi5VMq7rK6zsG3p1dsvX+b87T8O\nHtviyjQsYoxV4FgIYe64IwAhhDLwNuDdWdagtefWu0b5ynf3AXDf8EEAfuo5Q3mW1JUWti3t79/E\ns/uf3fL8WRwHj21x5dLPYiYoPgH8Y4zxq6npBwb6sy+qjaw3OwMD/YxPTtHbc+IK6vjkVGG3oah1\nARwYHqenZ8Ps8COH9nHF+Ze2PH8Wx2E5yyzyvl2L8mp+9DEgxhjf28rEY2MTGZfTPgMD/dabkUat\n2/o2Ml2pzo7f1rexkNtQ9H27tbyNSiXODp9z+o5l1ZvFcWh1mUXftwuthWDreFiEEK4DjsUY39Pp\ndWttuOyiei/mude1tXwL25Zeft5uxvdPtjx/FsfBY1tcpVqtltnCQwgXAzcAu4BpYB9wFvUb3BNA\nDfi3GONvNFlMrds+QVhvNrqpVrDeLHVTrQADA/1df5c+6xvcdwIvznIdkqTs+RfckqQkw0KSlGRY\nSJKSDAtJUpJhIUlKMiwkSUmGhSQpybCQJCUZFpKkJMNCkpRkWEiSkgwLSVJSXv0spI6r1mp8865R\nvnXP4wD85AVn8cLnDCXbdi5sP7p78BLKpc58zlpNm9Fqrcpto9/hOz94gOrhfi7Z/tyWtreV5ea1\nP5Qfw0Lrxq13jXLzrT9g4vAUAI//8AilUinZtnNh+1GAS4eel22xM1bTZvT20Tv4/N6vM3F4GoDH\n9x6mVHrhqtuU5rk/lB8/DmjdGB6bZKpyfHZ4qnKc4bF0s5+RydGmw1laWF8r9TaMTI7O6zpX6T20\nrPmbLbfZsNYmw0Lrxs6BPjbO6Tm9sWcDOwf6kvMN9Q02Hc7SwvpaqbdhqG9wXj/rnunTlzV/s+U2\nG9ba5GUorRuXXTRIDebds2ilbefC9qON4U5YTZvR3YOXUKN24p7FM5/bljalee4P5SfTtqptYlvV\nDHVTvd1UK1hvlrqpVlgbbVW9DCVJSjIsJElJhoUkKcmwkCQlGRaSpCTDQpKUZFhIkpIMC0lSkmEh\nSUoyLCRJSZk/GyqE8Czgs8D7Y4wfCiHsBD5BPahGgdfEGKezrkOStHKZnlmEEDYDHwRumTP6PcCf\nxRhfBDwAvCHLGiRJq5f1ZaijwFXUzyAaLgdunnl9M/DSjGuQJK1SppehYoxV4FgIYe7ovjmXnZ4A\nfBi+gJNbiL7g2Wez5+7HZh8p/pLnncNzzj9z0bagWbf6XF170xPzDm7fzP2PHmR4bJKnn7WF1/78\nBfSUy4tOO3c9c8fv2L6ZDWftY3TyscK3NbUF69qRdz+Lrn9sr9pnYQvR+x49yL2PHJxtg7r/0H08\ntXvXom1Bs271uZr2pnPn3fMvoxydOs6GconHfngYgDf+wo8l1zN3/PcPfJeeQ4/Qv7m38G1NbcG6\nduQRFhMhhFNijMeAHcBIaoaBgf7sq2oj612Z8cmpeZ3dRg8cplKtUpr5BH9s+jjjk1OL1ntgeJye\nOV3wDlTH27pdC2tbqo6FBgb65807fbze5rSxTaMHDs9bzlLrmTt+8pQnqVRrs9vbzm1t93shy+NS\nlPftepFHWNwCvAL4m5l/v5CaocuanFjvCm3r2zivZ/TO7X0cmpiiVqsAcErvBrb1bVy03q3lbVQq\ncd5wO7drYW1L1TFXY9/Onbd3Q5njx4/TaDo2uHXzvOUstZ6548vHTqPn9INUZvqJt2tbs3gvZHVc\nivS+bcVaCLZMwyKEcDFwA7ALmA4hvBK4Dvh4COHXgIeBj2dZg7rHwhaiS92zWEzWrT5X09507rwv\n+vGhk+5ZtLKeueN3bL+MDWedO++eRVHZgnXtsK1qm3XjJ55uqbebagXrzVI31Qq2VZUkrROGhSQp\nybCQJCUZFpKkJMNCkpRkWEiSkgwLSVKSYSFJSjIsJElJhoUkKcmwkCQlGRaSpCTDQpKUlHenPK1T\n1VqNb35/hDv2f5fy5gl+4twf4QWDP9Hxlpsrbfu5cL6fPPtibrv7cR4de4rp/ofZvPUIF0ycx7/r\n+zFuu/vxOY9dfxrfeuzOk+YbHptkx0Af1GoMj01y5FiFTZs2MN3/CIcqY1SPnMYl25/LC58zNK+d\n62ravRadLVmLxbBQLm69a5R/2HsrR7Y8AIdgbO8+SpQ63nJzpW0/F86399GD/OCeMzjS9xBHph+g\n/2gvj0w+wp1Tj/GDe84A6m1SHzz6rzxW/reT5gO4876x2eVPHJ7i1MFRKsfqNZXLJR7fe5hS6YXz\n2rmupt1r0dmStViMaeVieGySSu+h2eHpSpWRydGO17Fwna3WsHC64afqw41tanS1a4xfON1iw1OV\n47P/AUz3HqRKjSq12WUPj03Onz8x3M1WemyUDcNCudg50EfP9Omzw709ZYb6Wu8+1y4L19lqDQun\n27mlPtzYpka/7Mb4hdMtNryxZ8PsfwC902dQpkSZ0uyydw70zZ8/MdzNVnpslA0vQykXl100SK12\nGXfs3zx7zyKPlpsrbfu5cL6fPPtibtv0OI+Onc50/9b6PYvBmXsWm+bes7iYbz12xknzLX7PYpDp\n/u0n7lk887kntXNdTbvXorMla7HYVrXNurDdY9fU2021gvVmqZtqBduqSpLWCcNCkpRkWEiSkgwL\nSVKSYSFJSjIsJElJhoUkKcmwkCQlGRaSpCQf9yFJXSqE0AN8HfhUjPGGLNfV8bAIIfQBfwVsBTYC\n74kxfqnTdUjSGjAIHMk6KCCfy1CvA+6NMV4BvAr40xxqkKS14NeBZ4YQ3pL1ivIIi/3AtpnXZwJj\nTaaVJC3tI9Q/fH8k6xUlL0OFEH4uxviFdq0wxnhTCOF1IYS9wBnAy9q1bHWPdrUD7fbWm9ValdtG\nv8N3fvAA1cP9i7ZOzWa99f2/sA3shVsupFwqJ49Pq/u9HdMtVsvqt7+73zd5aOWexW+GEP4c+Gvg\nxhjjw6tZYQjhOuDhGONVIYSLgI8C9kpcZ9rVDrTbW2/ePnoHn9/7dSYOTwMs2jo1C439v7AN7MTT\njnLp0POSx6fV/d6O6Rar5eUvPW01m9/175sFOvL482RYxBh/PoSwFfgPwIdDCAAfAz4dYzy+gnVe\nBnxxZtl3hRCGQgilGOOSjTUGBvpXsJr8WG/a+OTUbDe5xnArdSyc5sDwOD0zneUADlTHC7X/U7Uc\nGB6nUq3R+NBePeXJlvfFajT2/+QpT1IqQaVa/9+vsf9Sx6fV/d6O6RarBVb3vi36+2aZOtKUqKVv\nQ8UYD4QQ/g8wRf2GyjuAd4UQ3hRjvH2Z67wf2A18JoSwC5hoFhRAtzU5sd4WbOvbONunujGcqmOx\nWreWt1GpxHnDRdn/rezbreVt9JRLNHqQlY+d1tK+WK3G/i8fO41a7xg95dJsPWNjE8nj0+p+b8d0\ni9UCq/u90On3TVZBNHOl52czWfgCrdyz+Gng9cCLgU8Db4wx3hNCOBf4DPDcZa7zfwA3hhD+CdgA\n/Noy59ca0K52oN3eenP34CXUqJ24Z7FI69QsNNaxsA3shVsunPfzpY5Pq/u9HdNl0Tq22983eUi2\nVQ0hfJP6Hff/G2M8tuBnvxtj/C8Z1ge2Vc1UN9XbTbWC9Wapm2qFtdFWtZV7Fi9s8rOsg0KSVAB+\nV0ySlGRYSJKSDAtJUpJPnZWkdSSE8AzgA8B26t9I3QP8doxxqtl8nllI0joRQigDnwL+OMa4O8bY\n+LP1P0jN65mFJBXU1W//XBnoB568+YZr2vGX2j8D3BNj/Oaccb8DVJeYfpZnFpJUQFe//XM7qf/h\n8z8CN1399s8NtGGxFwDfmzsixngsxjidmtGwkKRi+g1gx8zr84E3t2GZNer3KZbNsJCkYupLDK/E\nvcDz544IIWwMIVyYmtGwkKRi+mvqD28FOALc1IZlfhk4J4TwMpi94f1fgVenZjQsJKmAbr7hmtuB\nX6J+A/ram2+45q7VLnPmCd9XAr8WQvgW8HXgYIzxXal5/TaUJBXUzTdc8wjwSDuXGWN8HPj3y53P\nsNC6066WrvVlVblt5Dvc+cT3gBIXP+0iXjD4vKYtOivVCn9776cZfmqEnVuG+OULXk65VF5Vm89W\n2oTaSlSrYVho3WlXS1eot+f8wsO38NTUJABPHNlPiXLTFp1/e++nufOJ79enPzwGwI+ccd6q2ny2\n0iZ0jbUSVYcZFlp3hscmmw4vx8jkKNPHT3xFffr4NCOTo83X/9TIScOn9m46abnLrSM1fyvTSEvx\nHFTrzs6BvqbDyzHUN0jvht7Z4d4NvQz1Ne/ktnPL0EnDC+dJLWOxOlLzr3YdWt88s9C60842nbsH\nL6FWq827Z5Fq0fnLF7y8vv4F9yxg5W0+W2kTaitRrUayrWoB2FY1Q91UbzfVCtabpW6qFdZJW1VJ\n0toQQvgT4BLgbOp/EX4/8MMY4ytT8xoWkrROxBjfARBCeC1wYYzxd1qd17CQpIJ69U1vnX1E+Sev\n/XCu9wz8NpQkFdCrb3rrvEeUv/qmt7bjEeUrZlhIUjFl8YjyFTMsJKmYsnhE+YoZFpJUTFk8onzF\nDAtJKqBPXvvheY8o/+S1H171I8pXw29DSVJBffLaD7f9EeUAMcaPL3cezywkSUm5nFmEEK4DfhuY\nBq6PMX4+jzokSa3p+JlFCOFM4HrgUuAXgGs6XYMkaXnyOLN4KfDlGONh4DDwlhxqkCQtQx5hcS7Q\nF0L4HHAG8O4Y41dyqEPL1Go70izbd6aWvVSNjfGPjj3FVP/DHKqMUT3ST//R89m8qZenD2xZVXvV\nxWq7evvlTaY9Uc90/8Ns3nqEHVsajw0vzduGFzz7bG67+7HZaU/dephjlWNs6jmFHVuGVtCCdflt\nZdt1TNvV0rZaq7Jn5Nsrqsf2siuTR1iUgDOBXwTOA74K7MqhDi1Tq+1Is2zfmVr2UjU2xh/pe4jJ\nY/fPTl8ZO8iWIz/C3uFDS27PSmvr79/Es/ufvei0c+s5Mv0A/Ud7eXBzfXuOj+2Ytw33PXqQ4f2T\ns9NuPFKhwlG2bOzjwUM/OGkfpKykrWy7jmm7Wtr+00O3r7ge28uuTB5h8TiwJ8ZYAx4MIUyEELbH\nGPcvNcPAQH/nqmuDtVrv+OQUvT3lecOLzXtgeJyeng0nhqvjbdsnB6rNl71UjY3xk6c8SY0Tz2Mr\nnfoUlckqvT3lJben5doWbPcjh/ZxxfmXLjrt3HpKJahUa/T0bOBAdZzpyYF52zB64PCCaSuwoUSl\ndnx2nuXU3ew4LrWcdh3TVt9DKY8M71txPVm+P9eyPMLiS8DHQgjvo36G0dcsKIBua3KyZuvd1reR\n6Up13vBi824tb6NSifOG27FPBgb6k8teqsbG+PKx0yhtfmL257UjW+gpl5muVJfcnlYtrO2c03cs\nuby59dR6x+gpl6hUjrO1vI3jC7Zh5/Y+hvdPzk7bW+6hUqvQU9owO89y6l5qHzV7L7TrmLb6Hko5\n5/Qd/MvoyurJ6v3ZzFoIo46HRYxxJITwd8DtQI36w7LUBVptR5pl+87UspeqsfHvo2OnM9W/9cQ9\ni4H59yzaWdvl5+1mfP/kotPOrWe6f+v8exaDpXnbcOKeRX3axe5ZLMdK2sq265i2q6Xt5eftZmLi\n6Irqsb3sythWtc3W8plF3rqpVrDeLHVTrbA22qr6FQBJUpJhIUlKMiwkSUmGhSQpybCQJCUZFpKk\nJMNCkpRkWEiSkgwLSVKSYSFJSjIsJElJhoUkKSmPR5SvS+3qEFZ0czvAHTla4dRTenj6Wa11oWt0\nMNv31CiHD5xK78SutnSwS9e8ss5pzear1qp85cE9xMceWqKjn93a1F0Miw5pV4ewomts51OHp5k4\nPEX/5o3s3ddaF7pGB7OJw9NMHJ7m1KcOsHf4vJbmXY2Vdk5rNt/to3dw2+P/TKVyfNFl2q1N3caP\nMh0yPDbZdHitaGzXVOX4vH9b2d6RyVGA2eY4ld5DLc+7Go31LjW8kvlSy1zpOqW8GBYdsnOgr+nw\nWtHYro0zbSsb/7ayvUN99UY4jbabPdOntzzvajTWu9TwSuZLLXOl65Ty4mWoDmlXh7CiO9EB7uR7\nFimNjmWz9yx6d/H0H119B7tW17vczmnN5ts9eAn9/Zvm3bNoxzqlvNgpr826sINX19TbTbWC9Wap\nm2oFO+VJktYJw0KSlGRYSJKSDAtJUpJhIUlKMiwkSUmGhSQpybCQJCUZFpKkJMNCkpRkWEiSknIL\nixDCphDC/SGEX82rBklSa/I8s/gDYDzH9UuSWpTLI8pDCAG4APj7PNavlVvr7WHntjsd7Dub40/s\nYN/+wydta5H3gy1blYW8+lncALwNeF1O69cKrfX2sHPbnX5v5D4qT5zDqZPnnbStRd4PtmxVFjoe\nFiGE1wB7YowP108wSH4cGxjoz7yudlrL9Y5PTs12smsMd3J7s17XgeFxema6+1WqU1RPeZLeY/Xt\nnbutre6HPN4Lc7cB4EB1vOU6uum92021rgV5nFm8DDgvhHA1sBM4GkJ4NMb4laVm6LImJ2u63m19\nG2d7ZDeGO7W9ndi3W8vbqFQiAD3lEpVjp81u79xtbWU/5PVemLsNjeFW6uim92431QprI9g6HhYx\nxl9qvA4hvAt4qFlQqFjWenvYue1OB3eczfHT59+zaCjyfrBlq7JgD24tS7lUKsy1+SyUS+X51/eX\n2NQi74eTtkFqg1zDIsb47jzXL0lqjd+nkyQlGRaSpCTDQpKUZFhIkpIMC0lSkmEhSUoyLCRJSYaF\nJCnJsJAkJRkWkqQkw0KSlGRYSJKSfOrsOlat1fjyPz/MPQ+OF641aJaK3BK1E2y7qpUwLNaxW+8a\n5Rt3jzJdqRauNWiWitwStRNsu6qV8OPEOjY8Ntl0eK1ar9vdMDI52nRYWoxhsY7tHOhrOrxWrdft\nbhjqG2w6LC3Gy1Dr2GUXDdLfv2nePYv1oMgtUTvBtqtaCcNiHSuXSvzM83fx4+efmXcpHVXklqid\nYNtVrYSXoSRJSYaFJCnJsJAkJRkWkqQkw0KSlGRYSJKSDAtJUpJhIUlKMiwkSUmGhSQpybCQJCXl\n8myoEML7gBcCG4A/jjF+Jo86JEmt6fiZRQjhcuDHYoyXAlcBH+h0DZKk5cnjMtTXgFfNvD4IbA4h\nrJ+elspVtVbjG98f4W9v2cs3vj9CtVbLuySpK3T8MlSMsQYcmRl8E/APM+OkzK33lqrSSuXWzyKE\ncA3weuBnU9MODPRnX1AbWW92Vlvr+OQUvT3lecNZbn837Vvornq7qda1IK8b3FcCvwtcGWOcSE0/\nNpacpDAGBvqtNyPtqHVb30amK9V5w1ltfzftW+iuerupVlgbwdbxsAghnAa8D3hJjPFQp9ev9W29\nt1SVViqPM4trgW3AJ2dubNeAX40xDudQi9aZ9d5SVVqpPG5w/wXwF51eryRp5fwLbklSkmEhSUoy\nLCRJSYaFJCnJsJAkJRkWkqQkw0KSlGRYSJKSDAtJUpJhIUlKMiwkSUmGhSQpybCQJCUZFpKkJMNC\nkpRkWEiSkgwLSVKSYSFJSjIsJElJhoUkKcmwkCQlGRaSpCTDQpKUZFhIkpIMC0lSkmEhSUoyLCRJ\nSYaFJCmpJ4+VhhDeD+wGqsBvxRi/k0cdkqTWdPzMIoTw08AzYoyXAm8CPtjpGiRJy5PHZaiXAJ8F\niDHeC5wRQtiSQx2SpBblERZnA2NzhvfPjJMkFVQRbnCX8i5AktRcHje4R5h/JjEEjDaZvjQw0J9t\nRW1mvdnpplrBerPUTbWuBXmcWXwJeCVACOFiYF+McTKHOiRJLSrVarWOrzSE8EfAi4DjwNtijHd3\nvAhJUstyCQtJUncpwg1uSVLBGRaSpCTDQpKUlMuzoVoVQrgO+G1gGrg+xvj5nEtaUgjhDcBrgBr1\nvx25JMZ4Wr5VLS6E0Af8FbAV2Ai8J8b4pXyrWloIoQR8BHgWcAx4S4zxvnyrOlkI4VnUn07w/hjj\nh0IIO4FPUP9QNgq8JsY4nWeNcy2sd2bcbwJ/ApwRYzycZ30LLbJ/nw7cCPQCU8CvxBifyLPGhkVq\nfQHwPuq/y45Sfy+M51njchX2zCKEcCZwPXAp8AvANflW1FyM8cYY44tjjFcA7wI+nndNTbwOuHem\n1lcBf5pvOUnXAKfFGC+j/jyxG3Ku5yQhhM3Un3N2y5zR7wH+LMb4IuAB4A151LaYxeoNIbwGOAvY\nl1ddS1li/74X+EiM8XLqv5jfnkNpJ1mi1t+iHmZXALcDb86jttUobFgALwW+HGM8HGN8PMb4lrwL\nWobrqb+Ri2o/sG3m9ZnMf/xKET0T+BZAjPFBYNfM2UaRHAWuYv4fmF4O3Dzz+mbq7+miWKzeT8cY\nfz+nelIWq/etwKdnXo9Rfy8XwUm1xhivjTE+PPO+3QEM51XcShU5LM4F+kIInwshfC2EcEXeBbUi\nhPATwCNFOR1eTIzxJuq/cPcC/wS8I9+Kku4GrgwhlEMIATgP2J5zTfPEGKsxxmMLRvfNuez0BDDY\n4bKWtFi9Rf7j2CXqPRJjrIUQysDbgL/Jp7r5lngvEEK4ErgXOCvG+L87X9nqFDksStQ/Kfwi8Hrg\nY/mW07I3AX+ZdxHNzNwLejjG+EzqTwH+7zmX1FSM8QvUzyy+BvwmcA/d90yxbqu3K8wExSeAf4wx\nfjXvepqJMX4xxhiAGEL43bzrWa4ih8XjwJ4YY23m0sNECKFQnyaXcDmwJ+8iEi4DvggQY7wLGCrg\nZZ15YozXxxh/Ksb4NuDMIp+5zTERQjhl5vUO6s9F6wbd9Je6HwNijLHIl30JIfzinMFPUf9/sKsU\nOSy+BFwRQiiFELZRP6Xfn3dRzYQQBoGJGGMl71oS7qfeqZAQwi7qNRf2F0QI4aIQwkdnXv8ccEfO\nJbXqFuAVM69fAXwhx1qaWfhBodAfHBpmzpCPxRjfk3ctLfjDEMJFM6+fD8Q8i1mJQj/uI4TwZuqX\ndWrAe2OMf59zSU3NPBjxvTHGl+VdSzMzX529EXgasAH4/Rjj1/KtamkzZz0fBS4EjgDXxRgL9Y2d\nmWN/A7CL+tcj9wHXUf9W3CnAw8DrY4zHcytyjiXq/TLws9R/mX0buC3G+M7cipxjiXrPon4zeYL6\n74h/izH+Rm5Fzlii1t+h/q3Daerv4dcU/cPvQoUOC0lSMRT5MpQkqSAMC0lSkmEhSUoyLCRJSYaF\nJCnJsJAkJRkWkqQkw0KSlGRYaN0KIfznEML/nHkdQgj3zPx1u6QFDAutZx8AfjSEcCn1J+++uciP\n6ZbyZFgJecjNAAAAfUlEQVRo3Zp5eOIbgU8Cd8UYv5lzSVJhGRZa77ZRfxDdOXkXIhWZYaF1K4Sw\nCfgwcDUwFUL4lZxLkgrLsNB69m7qfafvB36Les+BoZxrkgrJR5RLkpI8s5AkJRkWkqQkw0KSlGRY\nSJKSDAtJUpJhIUlKMiwkSUmGhSQp6f8DsGKx9GAMW/sAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# F別の分布をみる\n", "sns.lmplot('x', 'y', data=d, hue='f', fit_reg=False )\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t私は箱ひげ図には馴染みがなく、最初にみたときには驚きました。\n", "\t\t慣れてくるとばらつきを知るには良いプロット方法だと感じるようになりました。\n", "\t

\n", "\t

\n", "\t\tこの分布の違いは、後で示す施肥別のヒストグラムをみると分かりやすいと思います。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAESCAYAAAD67L7dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADtBJREFUeJzt3W+MZuVZx/HvzDMs7gzbdJYOSkVoDXohJX1BQ2hqVai1\nsFCgFtga17bEP6G2lRfqNoSkFGiCVIS02NAaFMQGYze2Ak0tbbFoNY2VpP4Jwl5xTQUVDLMwlGX2\n3+zM+GJml1nYmYXdPc89O9f38+rsmbPP/cvu2d9z7znPc5+B2dlZJEl1DLYOIEnqL4tfkoqx+CWp\nGItfkoqx+CWpGItfkooZ6nqAiDgDuBe4NTNvj4gh4G7gVOB54LLM/EHXOSRJczqd8UfEMHAb8OCC\n3b8BPJ2ZZwNfBH6mywySpP11PePfCawDrl6w7yLgWoDM/OOOx5ckvUSnxZ+ZM8CuiFi4+w3ABRFx\nM/AU8OHMfK7LHJKkF7W4uTsAPJaZ5wL/DlzTIIMkldX5zd0D+D/g2/PbXweuW+rgPXumZ4eGel1n\nkqSVZmCxH7Qo/q8xd93/T4G3ALnUwRMT2/sQSZJWlrGxNYv+bKDL1Tkj4kzgFuAUYAr4X+CXmfuk\nz4nANuCDmTm+2GuMj29z+VBJepXGxtYsOuPvtPiPBItfkl69pYrfb+5KUjEWvyQVY/FLUjEWvyQV\nY/FLUjEWvyQVY/FLUjEWvyQVY/FLWhY2b36UzZsfbR2jhBZr9UjSy9x335cAOO200xsnWfmc8Utq\nbvPmR8l8jMzHnPX3gcUvqbm9s/2XbqsbFr8kFWPxS2rukksuPeC2uuHNXUnNnXba6UT81L5tdcvi\nl7QsONPvHx/EIkkrkA9ikSTtY/FLUjEWvyQVY/FLUjEWvyQVY/FLUjGdF39EnBERWyLiwy/Zf15E\nzHQ9vqSjg8sy90+nX+CKiGHgNuDBl+w/FrgaeLLL8SUdPVyWuX+6nvHvBNYBT71k/zXAZ4HdHY8v\n6Sjgssz91WnxZ+ZMZu5auC8ifhJ4c2Z+CVj0m2WS6nBZ5v5qsVbPrcBvvdKDR0eHGRrqdRhHUmvH\nHNPbb3tsbE3DNCtfX4s/Il4PBHBPRAwAJ0bEQ5l57mK/Z2Jie9/ySWrjggvewyOPPLJve3x8W+NE\nR7+l3jz7WfwDmfkk8BN7d0TE95cqfUk1uCxzf3X9qZ4zgVuAU4CpiLgUeG9mPjd/iCtvSgJclrmf\nXJZZklYgl2WWJO1j8UtSMRa/JBVj8UtSMRa/JBVj8UtSMRa/JBVj8UtSMRa/JBVj8UtSMS2WZZa0\njGzadA8PP/zd1jGYnJwEYGRkpGmOs846m/XrNzTN0DVn/JKWhd27d7F7966DH6jD5iJtkpaFjRuv\nAuDmm29rnGRlcJE2SdI+Fr8kFWPxS1IxFr8kFWPxS1IxFr8kFWPxS1IxFr8kFWPxS1IxFr8kFdP5\nIm0RcQZwL3BrZt4eET8G3AkcA+wGfiUzn+46hyRpTqcz/ogYBm4DHlyw+5PA5zPzHObeEH6nywyS\npP11falnJ7AOeGrBvt8Evjy/PQ6s7TiDJGmBTi/1ZOYMsCsiFu7bARARg8BHgOu7zLBcLIc1z5fL\neudQY81zablq8iCW+dL/AvA3mfnQUseOjg4zNNTrT7AOrV69il6v7b30vWudv+Y1a5rmgLk/j7Gx\n9jm0fOz99+F50b1WT+C6C8jM/OTBDpyY2N6HON276KLLueiiy5tm2Lve+U03fbppjr3Gx7e1jqBl\nZHp6BvC8OFKWegPt+xQ0IjYAuzLzhn6PLUnqeMYfEWcCtwCnAFMRcRlwArAzIh4CZoFHM/OjXeaQ\nJL2o65u73wPO7XIMSdKr4zd3JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakY\ni1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+S\nirH4JamYoa4HiIgzgHuBWzPz9og4CfgCc286TwHvz8yprnNIkuZ0OuOPiGHgNuDBBbtvAP4wM38O\n+E/gV7vMIEnaX9eXenYC65ib2e91DvCV+e2vAO/sOIMkaYFOL/Vk5gywKyIW7h5ZcGnnaeDELjPc\neON1TEw82+UQR429fw4bN17VOMnyMDq6lmuuua51DKnvOr/GfxADBztgdHSYoaHeIQ/w/PPP8cwz\nzzBwzOpDfo2VYnb+P3jPPr+9cZL2Zqd20OsNMja2pnUUzev15s5P/06616L4t0XEsZm5C/hR4Mml\nDp6YOLySmp6eYeCY1Rx36sWH9TpaWV7Ycj/T0zOMj29rHUXzpqdnAPw7OUKWegNt8XHOB4FL57cv\nBR5okEGSyup0xh8RZwK3AKcAUxFxGbABuDsirgQeB+7uMoMkaX9d39z9HnDuAX70ri7HlSQtzm/u\nSlIxFr8kFWPxS1IxFr8kFWPxS1IxFr8kFWPxS1IxFr8kFWPxS1IxFr8kFXPQJRsi4vzMdCE1qQM+\nL+JFPi9if10+L+KVrNVzVUR8FrgHuDMzH+8kiVTQxMSzPPPsVgZXt340Rnszg7MATOx4rnGS9mZ2\n7On09Q96tmXmBRExCvwi8Ln5p2ndBXw5M6c7TScVMLh6iNHzT24dQ8vIxANPdPr6r+gaf2ZOAH8B\n/DnwWuB3gX+NiLd2mE2S1IGDFn9E/GxE3AU8CpwJ/Fpmng28G/hcx/kkSUfYK7mweCPweeBD849L\nBCAz/ysiNnWWTJLUiVdyjf/tS/zs945sHElS1/wcvyQVY/FLUjEWvyQVY/FLUjEWvyQVY/FLUjF9\nXyAkIkaAPwNGgVXADZn5jX7nkKSqWsz4rwA2Z+Y7gMuBzzTIIElltSj+rcDx89trgfEGGSSprL5f\n6snML0bEFRHxH8wt+HZhl+NNTk4yO7WTF7bc3+UwOsrMTu1gcnK2dQypiRbX+DcAj2fmuoh4M/An\nwFmLHT86OszQUO+QxxscHDjk36uVbXBwgLGxNU0z9Hp+vkIH1usNdnZ+tnj6w08DXwfIzH+LiNdH\nxEBmHnD6NTGx/bAGW716mB1TcNypFx/W62hleWHL/axePcz4+LamOaanZ5qOr+VrenrmsM7Ppd40\nWkw3tgBvBYiIU4Bti5W+JOnIazHj/yPgzoj4W6AHXNkggySV1eLm7iTwvn6PK0ma450lSSrG4pek\nYlpc45c0b3Jykplde5h44InWUbSMzOzYw+TMZGev74xfkopxxi81NDIywu7BKUbPP7l1FC0jEw88\nwcjqkc5e3xm/JBVj8UtSMRa/JBVj8UtSMRa/JBVj8UtSMRa/JBVj8UtSMRa/JBVj8UtSMRa/JBVj\n8UtSMRa/JBVj8UtSMRa/JBVj8UtSMRa/JBVj8UtSMU0evRgRG4CNwBRwbWZ+rUUOSaqo7zP+iFgL\nXAu8DXg3cEm/M0hSZS1m/O8EvpmZ24HtwIcaZJCksloU/xuAkYi4D3gtcH1mfqvLAWendvDClvu7\nHOKoMDu9G4CB3qrGSdqbndoBDLeOAcDMjj1MPPBE6xjNzeyeBmBwVa9xkvZmduyB1d29foviHwDW\nAu8B3gg8BJyy2MGjo8MMDR36iXDCCWP0et7DBti6dSsArxs9rnGS5eA4jj/+eMbG1jRN4fn5on3n\n53FrGydZBo6j0/NzYHZ2tpMXXkxEXAH8cGZ+av7XjwDnZObWAx0/Pr6tvwFXsI0brwLg5ptva5xE\nejnPzyNrbGzNwGI/azHV+AbwjogYiIjjgZHFSl+SdOT1vfgz80ngL4F/BL4KfLTfGSSpsiaf48/M\nO4A7WowtSdV5V0mSirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+S\nirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SirH4JakYi1+SimlW\n/BHxQxGxJSI+0CqDJFXUcsb/ceCZhuNLUklNij8iAjgN+GqL8SWpslYz/luA3wYGGo0vSWUN9XvA\niHg/8J3MfHxu4r90+Y+ODjM01OtLtpWu15t7nx8bW9M4ifRynp/90/fiBy4E3hgRFwEnATsj4r8z\n81sHOnhiYntfw61k09MzAIyPb2ucRHo5z88ja6k30L4Xf2b+0t7tiPgE8P3FSl+SdOT5OX5JKqbF\npZ59MvP6luNLUkXO+CWpGItfkoqx+CWpGItfkoqx+CWpGItfkoqx+CWpGItfkoqx+CWpGItfkopp\numSDpPY2bbqHhx/+busYTEw8C8DGjVc1zXHWWWezfv2Gphm6ZvFLWhZWrTq2dYQyLH6puPXrN6z4\nGa725zV+SSrG4pekYix+SSrG4pekYix+SSrG4pekYix+SSrG4pekYix+SSqmyTd3I+L3gbcDPeCm\nzPyrFjkkqaK+z/gj4hzg9Mx8G7AO+HS/M0hSZS0u9fwdcPn89nPAcEQMNMghSSX1/VJPZs4CO+Z/\n+evAX8/vkyT1wcDsbJvOjYhLgKuBd2XmtsWOGx/ftiLeFJbDmud71zsfHV3bNAfUWPNcamlsbM2i\nV1KaFH9EnAdcD5yXmT/oewBJKqzvxR8RrwH+Hvj5zNza18ElSU0+zvk+4Hhg0/xN3VngA5n5Pw2y\nSFI5za7xS5La8Ju7klSMxS9JxVj8klRMk7V61H8RcSpzy2O8jrk1kr4DbMzM3U2DqbyI+APgLcCP\nACPAFuDZzLysabAVzJu7BUTEIPDPwEcy8x/m930GeD4zP940nDQvIj4IvCkzP9Y6y0rnjL+GXwAe\n21v68z4GzDTKI6khi7+G04B/WbgjM3c1yiKpMW/u1jDL3HV9SbL4i9gMnL1wR0Ssiog3NcojqSGL\nv4ZvAidHxIWw72bvp4D1TVNJasLiL2D+eQfnAVdGxD8B3waey8xPtE0mqQU/zilJxTjjl6RiLH5J\nKsbil6RiLH5JKsbil6RiLH5JKsa1eqRDFBGbgB8HLs7MJ1vnkV4pi186dO8FRlzwTkcbL/VIhyAi\n7mDu388DEXFS6zzSq+E3d6VDFBHTwND8khjSUcMZv3R4BloHkF4ti1+SirH4pUPnbF9HJYtfOnRe\n29dRyZu7klSMM35JKsbil6RiLH5JKsbil6RiLH5JKsbil6RiLH5JKsbil6Ri/h/Dz2/ENpr6xAAA\nAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# F別の箱ひげ図\n", "sns.boxplot(x='f', y='y', data=d)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

可視化から見えてくるもの

\n", "\t

\n", "\t\t施肥別の散布図からは、久保本でも整理してあるとおり、\n", "\t\t

サイズxが増加すると共に種子数yが増えているように思える
施肥fの効果はないように思われる

\n", "\t\tのように感じます。\n", "\t

\n", "\t

\n", "\t\tこのようにデータを整理することが大切です。Sageでは処理と文章が1つのドキュメントにまとめられており、\n", "\t\tどのようにして結果を導いたのかを後から簡単にフォローすることができます。\n", "\t

\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

施肥の違いをもう少しみてみる

\n", "\t

\n", "\t\tPandasのデータフレームもRのData.Frame同様にデータの絞り込みができます。\n", "\t

\n", "\t

\n", "\t\t施肥のあるもの$f_T$を、以下の様にfのフィールドが'T\"のものとして、d_Tにセットすることが1行で記述できます。\n", "\t

\n",
    "\td_T = d[d['f'] == 'T']\t\t\n",
    "\t

\t\t\n", "\t\tまた、Pandasでは、列（Series）をドット(.)で指定することもできるので、以下のようにも記述できます。\n", "\t

\n",
    "\td_T = d[d.f == 'T']\t\n",
    "\t

\t\t\n", "\t

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

	y	x
count	50.00000	50.000000
mean	7.88000	10.370600
std	2.65453	0.944307
min	3.00000	8.520000
25%	6.00000	9.757500
50%	7.50000	10.225000
75%	9.00000	10.947500
max	15.00000	12.400000

\n", "

" ], "text/plain": [ " y x\n", "count 50.00000 50.000000\n", "mean 7.88000 10.370600\n", "std 2.65453 0.944307\n", "min 3.00000 8.520000\n", "25% 6.00000 9.757500\n", "50% 7.50000 10.225000\n", "75% 9.00000 10.947500\n", "max 15.00000 12.400000" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 施肥の有無の違いをみる\n", "d_T = d[d['f'] == 'T']\n", "d_T.describe()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "y 7.046531\n", "x 0.891716\n", "dtype: float64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d_T.var()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t同様に施肥のないもの$f_C$を、抽出します。\n", "\t

\n", "\t

\n", "\t\t平均と分散がほぼ同じであることから、どちらもポアソン分布と仮定しても良さそうです。\n", "\t

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

	y	x
count	50.000000	50.000000
mean	7.780000	9.807600
std	2.620874	0.999813
min	2.000000	7.190000
25%	6.000000	9.115000
50%	8.000000	9.880000
75%	10.000000	10.412500
max	12.000000	11.930000

\n", "

" ], "text/plain": [ " y x\n", "count 50.000000 50.000000\n", "mean 7.780000 9.807600\n", "std 2.620874 0.999813\n", "min 2.000000 7.190000\n", "25% 6.000000 9.115000\n", "50% 8.000000 9.880000\n", "75% 10.000000 10.412500\n", "max 12.000000 11.930000" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d_C = d[d['f'] == 'C']\n", "d_C.describe()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "y 6.868980\n", "x 0.999627\n", "dtype: float64" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d_C.var()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t列の度数分布は、以下の様にgroupbyでy毎の頻度を求めることができます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "y\n", "2 1\n", "3 2\n", "4 3\n", "5 4\n", "6 10\n", "7 1\n", "8 5\n", "9 8\n", "10 9\n", "11 4\n", "12 3\n", "dtype: int64" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 施肥なし（C）の度数を調べる\n", "d_C.groupby('y').size()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

施肥別のヒストグラム

\n", "\t

\n", "\t\tpandasのplot機能を使えば、施肥別のヒストグラムも以下の1行でプロットできます。\n", "\t

\n", "\t

\n", "\t\tbyでfを指定することで施肥別のヒストグラムを図化できます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEGCAYAAACNaZVuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFPFJREFUeJzt3X2QXXV9x/H3skBkEwxJTCqQmQSi+QrjqLVjB58REJuK\nYouMDhRUhhlsZYqDUqtWEGyrhUYFlVbjwyDqFEcriqggFbQd9A/HdrQOfiOGhMfJAwlrYBMMu9s/\n7o1ZaLK79+w55z6c92uG4e7N3t/5/vbufu65557f+Q5NTk4iSWqWg7pdgCSpfoa/JDWQ4S9JDWT4\nS1IDGf6S1ECGvyQ10MHdLkDqNRFxMfA2Wn8fBwO3AO/LzN92tTCpRO75S1NExD8BZwKvzszjgOcB\n84CbulqYVLIhF3lJLRGxCHgAeH5m/nrK/YfSejG4uWvFSSXzsI+0zwnAfVODHyAzfwcY/BooHvaR\n9lkMbO52EVIdDH9pn23A0d0uQqqD4S/t8xPgDyLiBVPvjIiDI+LvI+JpXapLKp3hL7Vl5ihwFfDF\niFgFEBEjwGeAF2Tm7m7WJ5XJs32kp4iIC4G/pLVzNAF8E/hg+4NfaSDMKvwj4rnAjcBHM/PaiFgO\nXE/rj+Mh4JzM3FNppZKk0sx42Kf9tvca4LYpd18BfCIzXwn8BjivmvIkSVWYzTH/3cAaWnv4e53I\nvhWPNwGnlFuWJKlKM4Z/Zk5k5uNPuXv+lMM8W4AjS69MklSZMs72GSphDElSjYpe3mFnRMxrvyM4\nGnhwum+enJycHBryNaJK69ev55z3foWRhctKH3tsdAvXf/gsVq9eXfrYkkrTUcgWDf/bgDOAr7T/\n/71pKxoaYuvWnQU31TuWLj28Z+exffujjCxcxoJF1SxQ3b790Z6aey8/F7PlHHrHIMxj6dLDO/r+\nGcM/Il4IrAVWAHsi4o3A2cB1EXEBsAm4rvNSJUndMmP4Z+bPgFft559OLb8cSVIdvLyDJDWQ4S9J\nDWT4S1ID2clL6sD4+DgbN24oZawdOxawffujv/965cpjGR4eLmVsaSaGv9SBjRs3cNFV3yp9PcXY\n6BauvuT1rFr17FLHlQ7E8Jc6VOV6CqkuHvOXpAYy/CWpgQod9omI+cAXgUXAocAVmXlrmYVJkqpT\ndM//rcCvMvMk4Ezg6tIqkiRVrmj4bwOWtG8vBraWU44kqQ6Fwj8zbwBWRMSvgTuAd5dZlCSpWkWP\n+Z8NbMrMNRHxPOBzwItKrUyagzIXY011772bSh9T6oai5/m/FLgFIDN/HhFHRcRQZk4e6AGdXmu6\nV/XqPHbsWFDp+IsXL+i5uU9Xz/r16ytZjPXw/XexZPlxpY65Vy/+jGejH2ven0GZx2wVDf+7gROA\nb0TECmDndMEP9H2jBOjthg9TLxNQ1fi9NPeZnouqmtuMjW4udbypeu1nPBu9/DfRiUGYR+nNXA7g\n08DnI+IOYBi4oOA4kqQuKBT+mfkY8KaSa5Ek1cQVvpLUQIa/JDWQ4S9JDWT4S1IDGf6S1ECGvyQ1\nkOEvSQ1k+EtSAxn+ktRAhRu4t6/seQmwB7g0M79bWlWSpEoV2vOPiMXApcBLgNOA08ssSpJUraJ7\n/qcA38/MMWAMeHt5JUmSqlY0/FcC8yPim8ARwOWZ+YPSqpIkVapo+A/R6t37BuAY4HZgxXQPGJRG\nCb06D5u5PFnVP48q9OLPeDb6seb9GZR5zFbR8N8M3Nlu4LIhInZGxDMyc9uBHtDvjRKgtxs+2Mzl\nyar+eVSh137Gs9HLfxOdGIR5dPriVfRUz1uBkyJiKCKWAPOnC35JUm8pFP6Z+SDwNeAnwM3AhWUW\nJUmqVuHz/DNzHbCuxFokSTVxha8kNZDhL0kNZPhLUgMVPuYvzdX4+DgbN24o9NgdOxZMezrnvfdu\nKlqW1AiGv7pm48YNXHTVtxhZuKz0sR++/y6WLD+u9HGlQWH4q6tGFi5jwaKjSx93bHRz6WNKg8Rj\n/pLUQIa/JDXQnMI/Ip4WEXdHxLllFSRJqt5c9/w/ADxcRiGSpPoUDv+ICOA5tK7tI0nqI3M522ct\n8A7greWU0gxzObd9Op7X3t8mJyYqfQ5XrjyW4eHhysZX/ykU/hFxDq3r+W9qvQFgaKbHDEqjhLnO\nY/369ZWc2171ee1VNBrpx4YrVdm1cytrb9jGyMKHSh97bHQL13/4LFavXl362ODfdr8quuf/WuCY\niHgdsBzYHRH3TdfKsd8bJUA5DR+2b3+0knPbqz6vvYpGI/3YcKVKVa15gOoaxQxCExQYjHl0+uJV\nKPwz8817b0fEZcA99vCVpP7hef6S1EBzvrxDZl5eRiGSpPq45y9JDWT4S1IDGf6S1EBe0nk/DrQQ\na6YGIrPhYizVrcoFZIsXP7+ScVU9w38/bDKiQVLVArLW4rEFLFp0ZKnjqh6G/wHYZESDpMoFZOpP\nHvOXpAYy/CWpgQof9omIK4GXAcPARzLzG6VVJUmqVKE9/4g4ETg+M18CrAE+XmZRkqRqFT3s80Pg\nzPbtR4CRiJjxss6SpN5Q9Kqek8Cu9pfnA99p3ydJ6gNzOtUzIk4H3gacWk45nbErliQVM5cPfF8D\nvBd4TWbO2AWhii45/doVqx/ZyUsHMigdsAZlHrNVtI3j04ErgZMzc3Q2j6miS06/dsXqR3by0oH0\newcssJNXJ94ELAG+2v6gdxI4NzPvLzieJKlGRT/wXQesK7kWSVJNXOErSQ1k+EtSAxn+ktRAXtJZ\nM6qqGYjrKaTuqSX83/PBj7Fx8+9KH3fLfXfBsj8ufVw9WVXNQFxPIXVPLeE/9sQhjI0cW/q4jx+y\nhXmlj6r9cT2FNFg85i9JDWT4S1IDGf6S1EBzubDbR4ETgAngnZn509KqkiRVqmgnr1cAz2p38jof\nuKbUqiRJlSp62Odk4EaAzPwVcEREeH1eSeoTRQ/7PBOYephnW/u+u/f3zbseG2XskScKburAHv/t\nVsYPKv8a3Lt2bgeq6UpZ1djWXM/Y/TZulWOPjW4pfUzVp6zz/Kf9zfrsx95vf19JPa1pzVyKHvZ5\nkNae/l5HAeUu/5QkVaZo+N8KvBEgIl4IPJCZj5VWlSSpUkOTk5OFHhgR/wi8EhgH3pGZvyizMElS\ndQqHvySpf7nCV5IayPCXpAYy/CWpgQx/SWogw1+SGsjwl6QGMvwlqYEMf0lqIMNfkhrI8JekBjL8\nJamBDH9JaiDDX5IayPCXpAYy/CWpgQx/SWogw1+SGsjwl6QGMvwlqYEMf0lqIMNfkhrI8JekBjL8\nJamBDH9JaiDDX5Ia6OBuFyD1ioi4FnhV+8tVwAPAbmASeFFmPtat2qSyDU1OTna7BqnnRMQG4OzM\n/HG3a5Gq4GEfaf+G2v9JA8nwl6QGMvwlqYEMf0lqIMNfkhrI8JekBjL8pf3zHGgNtFmd5x8RzwVu\nBD6amddGxHLgelovHg8B52TmnkorlSSVZsY9/4gYAa4Bbpty9xXAJzLzlcBvgPOqKU+SVIXZHPbZ\nDayhtYe/14nATe3bNwGnlFuWJKlKM4Z/Zk5k5uNPuXv+lMM8W4AjS69MklSZMi7sNuMS+MnJycmh\nIVfKd2r9+vWc896vMLJwWeXbGhvdwvUfPovVq1dXvi1JlegoZIuG/86ImNd+R3A08OC0FQ0NsXXr\nzoKb6h1Llx5e6zy2b3+UkYXLWLDo6Nq21y/PU93PRRWcQ+8YhHksXXp4R99f9FTP24Az2rfPAL5X\ncBxJUhfMuOcfES8E1gIrgD0R8UbgbOC6iLgA2ARcV2mVkqRSzRj+mfkz9jW4mOrU8suRJNXBFb6S\n1ECGvyQ1kOEvSQ1k+EtSAxn+ktRAhr8kNZDhL0kNZPhLUgMVurZPRMwHvggsAg4FrsjMW8ssTJJU\nnaJ7/m8FfpWZJwFnAleXVpEkqXJFw38bsKR9ezGwtZxyJEl1KBT+mXkDsCIifg3cAby7zKIkSdUq\nesz/bGBTZq6JiOcBnwNeNN1jOr3WdK+qcx47diyobVsAixcv6KvnqZ9qPRDn0DsGZR6zVbSZy0uB\nWwAy8+cRcVREDGXm5IEe0O+NEqA7zVzqZDOXejmH3jEI86irmcvdwAkAEbEC2Dld8EuSekvRPf9P\nA5+PiDuAYeCC0iqSJFWuUPhn5mPAm0quRZJUE1f4SlIDGf6S1ECGvyQ1UNEPfDVgJicmuPfeTbVs\na3x8HBhieLj4vseOHQtmdSpsGdvqxMqVxzI8PFzLtqS5MPwFwK6dW1l7wzZGFj5U+bYevv8uDjt8\nCSMLlw3UtsZGt3D1Ja9n1apnV74taa4Mf/3eyMJlLFh0dOXbGRvdPJDbkvqJx/wlqYEMf0lqoMKH\nfdoXd7sE2ANcmpnfLa0qSVKlCu35R8Ri4FLgJcBpwOllFiVJqlbRPf9TgO9n5hgwBry9vJIkSVUr\nGv4rgfkR8U3gCODyzPxBaVVJkipVNPyHaLVvfANwDHA7sKKsoiRJ1Soa/puBO9vX8N8QETsj4hmZ\nue1ADxiULjmD3MlLc9dpN7RB+LsYhDnA4MxjtoqG/63AFyLiSlrvAOZPF/xgJ68i6u7kpbnrpBva\noHSP6vc5wGDMo5ZOXpn5IPA14CfAzcCFRcaRJHVH4fP8M3MdsK7EWiRJNXGFryQ1kOEvSQ1k+EtS\nA3lJ5w6Nj4+zceOGWrZVV3MVSc1j+Hdo48YNXHTVt2prRLJk+XGVb0dS8xj+BdTZiESSquAxf0lq\nIMNfkhrI8JekBppT+EfE0yLi7og4t6yCJEnVm+ue/weAh8soRJJUn8LhHxEBPIfWhd0kSX1kLnv+\na4GLaTV2kST1kULn+UfEObSauWxqvQGY+QVgUBolLF5sgxUdmM1c+tegzGO2ii7yei1wTES8DlgO\n7I6I+6br49vvjRKg9cthgxVNx2Yu/WkQ5tHpi1eh8M/MN++9HRGXAffYwF2S+ofn+UtSA8352j6Z\neXkZhUiS6uOevyQ1kOEvSQ1k+EtSAw3E9fwnJiZY+6nPccihI5Vu57CRQ9l0z0bgmZVuR/1pcmKi\no+5rO3YsKHzq8Pj4ODDE8HA9+28rVx7L8PBwLdtSPQYm/P/nnjHmLVtV7YZ2wPZtWzn0sGo3o/60\na+dW1t6wjZGFD1W+rYfvv4vDDl9SS0e5sdEtXH3J61m16tmVb0v1GYjwl3pFnV3e6tqWBpPH/CWp\ngQx/SWqgwod9IuJK4GXAMPCRzPxGaVVJkipVaM8/Ik4Ejs/MlwBrgI+XWZQkqVpFD/v8EDizffsR\nYCQivK6/JPWJolf1nAR2tb88H/hO+z5JKmx8fJyNGzfUsp2p6yTmsuZiNnpxncScTvWMiNOBtwGn\nzvS9VTZKeOKJJzjoIN94SFWZrklNmX/b69ev56KrvlX5+oW610lc/+GzWL16deXb6sRcPvB9DfBe\n4DWZOWMXhCobJTzxxBNMTPjGQ6rKgZrUlN0EZfv2R2tZv1D3OolOmvwUVUszl4h4OnAlcHJmjhYZ\nQ5LUPUX3/N8ELAG+2v6gdxI4NzPvL60ySVJlin7guw5YV3ItkqSauMJXkhrI8JekBjL8JamBvKSz\npGlN16Sm7MVRnTTD0dwY/pKmVXeTmiXLj6t8OzL8Jc1CnU1qVA+P+UtSAxn+ktRAc7m2z0eBE4AJ\n4J2Z+dPSqpIkVapoM5dXAM9qN3M5H7im1KokSZUqetjnZOBGgMz8FXBERCworSpJUqWKhv8zga1T\nvt7Wvk+S1AfKOtWzq51UhoaGOHLB48w76J5KtzNv3sEccsgWHhjdNfM3l2DXzu3U9aN1W26rSduq\nc05jo1tq2U6niob/gzx5T/8oYLoVIENVdvICuO5fPlTp+JI0SIoe9rkVeCNARLwQeCAzHyutKklS\npYYmJ4u1P4yIfwReCYwD78jMX5RZmCSpOoXDX5LUv1zhK0kNZPhLUgMZ/pLUQIa/JDVQZeEfEf9v\nBUVELK9qe1WLiJO6XUMZIuIZ3a5hrvb3uyWpM6Wf7RMRfwZ8HBgBvgNcmJk72//2g8zs+RCNiHOf\nctcQ8HfAhwAy84u1F1VARKwBTs/Mt7dfvL4A7ATm03pebu5qgbMQEacCV9O6nMi7gU/RWlS4E7gg\nM3/YxfJmLSIOAc4DTgGObN/9IPA94LrMHO9WbbMVEcuAdwGLga9k5u1T/u2TmXlh14qbg4i4PTNf\n1e06OhERC4GXZ+a3I+II4H3A8UACH8nMrdMOQDWdvP4W+EPgEVpX/Lw1Iv4kM0fp8mUgOnAp8DBw\nM/tqfhpwTNcqKuYK4LT27cuAV2XmhohYQmtuPR/+tJ6Lk2gFzh3AyZn584hYAXwJeHkXa+vE9cBv\ngLXAFlq/V0cDZ9B6UX7qDkcv+hKtCzr+FLgsIl6WmXuX1h/fvbJmLyImaL3o/o59f9tHRsQ9wGRm\nHtu14jrzNeCG9u1rgV8CHwT+CLgO+NOZBqgi/Mczc3v79mciYjNwS0ScBvTLooLnAh8Ang9cnJmb\n2i9gl3e5rk4dQmsPGVovxnsvflTfhU3m7neZ+RDwUEQ8kpk/B2g/Jz2/tzzFkZn55qfc9xvgRxHR\nF+9egEMz81qAiPg6cH1EXJqZV9A/v09raO2gfjIzvw4QET/OzBd3t6yOPT0zP9u+fWRmntW+/dOI\n+IvZDFBF+P9XRHwbODMzd2XmNyNiN/AfwJIKtle6zNwNvD8iAvhURNxJf344fhXw3xHxfVqBf2N7\nLicBn532kb1jR0T8A63fnbsj4l+BW2g1Euqnhq8TEfHnwE2ZuQcgIubR2vN/vKuVzd6eiDgD+PfM\nnIiIc4AvRMRngGov3lWSzLwlIu4A3tcOyYvpn53Sqe6OiI8BXwZuj4gzgR/RenGb7jprv1d6oGXm\n3wD/DOyect8ttN6e99Wec7acBtzHvr3mvpGZXwZeTOuX4i7gTlqB+bbMXNfN2jpwLq236bdn5hrg\nP4FX05rHed0srEPn0DoElxGxuf2O+JfAK4C3dLWy2TsPeB2tQ6Bk5kRmvgX44d77+kFmPp6ZlwHv\nAT4BLAVoHzvvF28B1tM6tPvm9v//DVgG/PVsBvDyDlKX9cuJENPp5zlExFGZ+WA/z2Gq2c6jisM+\nkp4iIv5qmn8+urZC5mCQ59A+fbgv5gDlPBeGv1SPi4Hb2P/x2ENqrqUo59A75jwPw1+qxxuAa4CL\nMvNJH/BGxIldqahzzqF3zHke/XgGi9R3MvN/aX3gu2c///yumsspxDn0jjLm4Qe+ktRA7vlLUgMZ\n/pLUQIa/JDWQ4S9JDWT4S1ID/R9apL1ek21CsAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# F別のヒストグラムをみる\n", "d.hist(column='y', by='f', sharex=True, layout=(2,1), bins=11)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

ポアソン回帰の統計モデル

\n", "\t

\n", "\t\t施肥別の平均、分散からも$個体_i$の種子数$y_i$の確率$p(y_i | \\lambda_i)$は、ポアソン分布に従うとして、\n", "\t\t以下の式で表します。\n", "$$\n", "\t\tp(y_i | \\lambda_i) = \\frac{\\lambda_i^{y_i}exp(-\\lambda_i)}{y_i!}\n", "$$\t\t\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

$個体_iの平均種子数\\lambda_i$の式

\n", "\t

\n", "\t\t久保本の3.4.1では、$個体_iの平均種子数\\lambda_i$を以下の様な指数関数で表しています。\n", "$$\n", "\t\t\\lambda_i = exp( \\beta_1 + \\beta_2 x_i)\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tここで、$\\beta_1, \\beta_2$は、式を決定づけるパラメータです。図3.4では、このパラメータを変えることで、\n", "\t\tどうような曲線になるか示しています。これをSageを使って表示すると以下のようになります。\n", "\t

\n", "\t

\n", "\t\t直線よりも指数関数の表現力があるようにみえます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEBCAYAAACXArmGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4E+XawOHfTJK26ZJSKEsrtMguAgKioCKgqKAIKsgi\nIqvHDcEFRRSrgB5R2VQQFxA8B3FDccEFFxBFyhH9BBdEQARaFoHSJW3apsnMfH8ECghIaNNOlue+\nrlxpwmTy5GUyT951FMMwDIQQQghANTsAIYQQwUOSghBCiHKSFIQQQpSTpCCEEKKcJAUhhBDlJCkI\nIYQoJ0lBCCFEOUkKIiQZhoHT6USm2QgRWFZ/NzxwoPCEz6uqQs2aceTmutB1+YKeDim7inO5Cjnz\nzDPYvn03cXEJZocTUuS4q7hQL7vatU/9Xal0TUFVFRRFQVWVyu4q4kjZVZyiKMfcC//JcVdxkVB2\n0nwkhBCiXKWTQlaWwg8/BCIUIYQQZqt0Uli82ErfviD9fUIIEfoqnRTattXJzobs7PBtYxNCiEhR\n6aTQqZMGQGampdLBCCGEMFelk0JSErRuDZmZ0mcthBChLiBn8osvhrVrpaYghBChLiBJoUsX2LZN\nZd8+6VcQQohQFrCaAsB330ltQQghQllAkkJqKjRqpEsTkhBChLiA9Q5fcIEmSUEIIUJcwJLChRfq\nbNqkkpcXqD0KIYSobgFMChqGobBundQWhBAiVAUsKaSlGaSm6qxZ4/dq3EIIIYJMwJKCokDnzhpr\n1khNQQghQlVApyF37uzl11+lX0EIIUJVQJPCRRf5+hUyM6UJSQghQlFAk0KDBgYNG+p8+600IYlT\n+/XXX+jXrw9NmjSgVaum3HLLcPbv3292WEJEtICvYnfxxV5JCuKUysrKGDjwOi6+uAubNv3JN9/8\njwMHDvDAA/eaHZoQES3gSeGiizQ2b7awf7+sgyROrqSkmIkTH2Xs2Hux2WzUrFmLXr168/vvv5kd\nmhARrUqSAiCjkMQ/SkysweDBN6GqvkPwjz+28uabr3Pttf1MjkyIyOZ3j7CqKqjq8b/+LRb1mPsz\nzoDmzX3zFfr31wMUZnj6e9lFol27sunQ4Rw0TWPYsBE89NDDKMqpa5lHl53VGrnlVxFy3FVcJJSd\nYhj+XV3ZMAy/vqwAY8bAp5/CH39UKjYRQbZt28Ytt9xCvXr1WLx48Sm3dzqdJCYmUlBQgMPhqIYI\nhYgMfieFgweLTlpTcDjsOJ0laJqvZvDRRxaGDo3h55+LqV/fr91HpBOVXST7/vt19OzZna1bd1Cz\nZq1/3NblKqJBg3pkZ/9FXFx8NUUYHuS4q7hQL7ukpLhTbuN385GuG+j6yU/wmqbj9foKqWNHHUWJ\nZtUqhUGDNH/fImIdXXaR4ttvv2H8+HvIzPy/8ud03VcbVVXrKcvj8BcyEssuUKTsKi6cy65KGsaS\nkqBVK51vv5VJbOLEzjmnLYWFhUyZ8gglJSXk5OQwffqTXHDBRcTHJ5gdnhARq8p6Szp31li92oJ/\njVMi0iQkOFiy5APWr/8/zjqrEV27diIxsQYvvviK2aEJEdGq7Kd8165eXnghiq1bVZo1C89qlqic\nFi3O4r33PjY7DCHEUaqsptCpk0Z0tMGqVTJfQQghQkWVJYXYWOjYUWPVKulXEEKIUFGlMzC6dfOS\nmWnB7a7KdxFCCBEoVZoULrlEo7hYLtEphBChokqTQsuWOrVr63z1lSQFIYQIBVWaFBQFunWTfgUh\nhAgVVb6qU7duXn79VZbSFkKIUFDlSaFrV98yF19/LU1IQggR7Ko8KdSpY9CqlTQhCSFEKKiWRcG7\ndfOyapUseSGEEMGuWpLCJZdoHDig8uuv4XthCiGECAfVcpbu2FEjPt7gyy+lCUkIIaqEYWBd9x14\nPJXaTbUkhago3wJ5n38uSUEIIapC1PJPSLr6cix/bK3UfqqtPeeKK7z8+KNKTo4MTRVCiIAqLSU+\n40Hc3S9Ha3FWpXZVbUnh0ks1DENh5UoZmiqEEIEU+8Js1L27cT32pG/WcCVUW1KoW9egbVuNL76Q\nJiQhhAgUdfcuYp+dQcm/bkdr0rTy+wtATH67/HIvX31lrWw/iBBCiEPipmRgxMVTPG58QPZX7UnB\n6VT4/ntpQhJCiMqyrV1DzHvvUpQxGSPBEZB9VmtSaNPGt2qqjEISQohK8nqJf/B+POd2wD3ghoDt\ntlqTgqr6agtffik1BSGEqAz7/BexbNpI0dTpvpNrgFT7FOPLLtPYssXCjh0yNFUIISpC3bOb2Kee\noHTEzXjbtg/svgO6Nz906+bFZpPZzeFq165shg+/kRYtGtKqVVPGjr2dwkLncdtNmzaVlJQk0tPr\nkp5el7S0OqSn1yUnJ8eEqIUILfEZD0JsLK4HMwK+72pPCvHxcOGFGsuXS1IIR0OGDCQpKYn16zfx\n5ZffsHnzJiZNeviE2w4YcAM7d+5j5859ZGXtZ+fOfSQnJ1dzxEKElqgVnxO97H2KpjyBkVgj4Ps3\nZYW6K6/0kplpIT/fjHcXVcXpLKBdu/ZMnDgJu91OvXopDBgwmLVr15gdmhDhoaSE+An3UXZxN9x9\n+1fJW5iWFLxeRUYhhRmHI5FZs+Yc82t/9+5dpKSknnD7jRt/pVevy2ncuD5du3Zi1aqV1RWqECEp\n9tnpqHv3UPT0jErPXD4Zv8/KqqqgqscHYbGox9z7o0EDOPdcjeXLbQwerPv9unBTkbILJevX/8iC\nBS/zxhtLsFqP/Yz169enUaNGPProFOrWrcfCha8wZMgA1qxZR+PGTU6576PL7u/7Fv8s3I+7qmRm\n2albNhM7+xlK77oXpXlz/0/ep0kxDP8ufWMYBkoAM9NTT8GUKZCTA3Z7wHYrgsSaNWvo06cPU6ZM\nYfTo0X69plOnTvTo0YPJkyefclun00liYiIFBQU4HIGZtCNE0DIM6N4dsrLgl1+q9KTpd7LJzXWd\ntKbgcNhxOkvQNP9/9V9yicKECbEsXVrKVVdpfr8unFS07ILd8uWfcNtt/2LatJn07z+QvDyXX69L\nTa3P9u1Zfm3vcpUAHCo7mfdyOsL1uKsOZpVd1JK3iPvqKwqXvI+3VIdS/75Tf5eUFHfKbfxOCrpu\noOsnr1Romo7X638hnXkmNGumsWyZhSuuiOzFkE637ILZunXfMXr0rSxYsIguXbqd9HPNmjWN887r\nSOfOXcqf27x5M9dd18+vsjj8hQynsqtuUnYVV51lp+TkYH9oPKXX9KW066VQxe9raqPiVVf5Lrzj\n9ZoZhQgUTdMYN24MGRlT6NKl23H/ftFFHVi37jsAcnNzmTBhHNu2bcXtdjN37mx27NjOwIGDqzlq\nIYJb/MPjQdcp+vfT1fJ+pieFvDyFtWul+h8Ovv9+HVu3bmHixPHlk9EO3+/alc22bX/gchUBkJEx\nmUsvvZx+/frQvHk6H3zwLkuXLqNevRSTP4UQwSPqs0+JWfoORY89iVGnTrW8p98dzQcOFJ7weatV\nJSkpjrw812lXpwwD2rePo2dPL1Onuk/rteGgMmUX6YqLi2jYMJUdO/YQGxtvdjghRY67iqvOslMK\n8km6uCPes1vhfP2dgAxBrV074ZTbmFpTUBTfnIVPP7XiX2oSQojIEDc5A6WoiKLpz1bZnIQTMX2g\n8lVXedmzR2X9etNDEUKIoGD7ZhX21/6D65Ep6GfUr9b3Nv1M3KmTRnKyzgcf2MwORQghzOdykXDv\nWMou7Ezp0BHV/vamJwWrFXr39vLBB1Z0ad4UQkS4uKlTUA/so3Dm7IBeJ8FfpicFgOuu8zUhyWU6\nhRCRzLruO+zzXsT1wMPojRqbEkNQJIXzz9dISdF5/31ZIE8IEaGKi0m463a87dpTcusdpoURFElB\nVaFPHy8ffmhFi8wVL4QQES5+SgaW3bsonP0SWMxrNQmKpABw7bUeDhxQycyUJiQhRGSxfbUC+4J5\nFD0yBa1pM1NjCZqk0L69TlqaNCEJISKLkpdLwl13UNblEkpH3mJ2OMGTFBQFrrnGw8cfW/FE9vp4\nQogIEj9hHEpJCYXPzTVltNHfmR/BUa691ktursrq1dKEJIQIf9HvvUPMe+9S9OR09NQzzA4HCLKk\n0KqVTuPGOu+/LxPZhBDhTd27h/jx91J6Td8qu95yRQRVUjjchPTJJ1bckbc+nhAiUhgGCXfdgRET\nQ9HTM6t1baNTCaqkAL6JbE6nwooV0uEshAhP9nkvELVqJYXPzsVIqml2OMcIuqTQvLlOq1YaS5ZI\nUhBChB/rLz8RN+URim+9A8+ll5kdznGCLikADBzo4fPPreTmmh2JEEIEUFERCbeMwNusBa6HJ5sd\nzQkFZVLo29eLYcB770mHsxAifMQ//ACWvXspfHkhREebHc4JBWVSqF3boHt3jbfflqQghAgP0e+9\ng/31RRQ+OR2tSVOzwzmpoEwKAAMGeFi/3sKWLUEbohBC+EXduYP4++6mtO/1uAcONjucfxS0Z9wr\nrvBSo4bB229Lh7MQIoSVleG4dQRGUk2Knp4VVMNPTyRok0J0tG+RvCVLbLJyqhAiZMVNmoj1119w\nzluI4Ug0O5xTCtqkAL5RSHv3yrIX4WrXrmyGD7+RFi0a0qpVU8aOvZ3CQqfZYQkRMNHvvUPs/Jco\neuxJvO3ONTscvwR1UmjfXqdJE4233pIO53A0ZMhAkpKSWL9+E19++Q2bN29i0qSHzQ5LiICwbNlM\nwj1jKO3bn9Lho8wOx29BnRQUBQYO9PLJJ1YKC82ORgSS01lAu3btmThxEna7nXr1UhgwYDBr164x\nOzQhKs/lwjHqJrT69Smc/mzQ9yMcLaiTAkD//h7cbmSRvDDjcCQya9YckpOTy5/bvXsXKSmpJkYl\nRAAYBgn33YUlOxvngtcgPt7siE6L30N7VFVBVY/PdhaLesx9oKWlwWWXabz2mo0RI8Krx7mqyy6U\nrF//IwsWvMwbbyzBaj11eRxddv5sL46Q467i/Cm7qIXziXn3bYrmLURpeZb/J9kgoRiGYfizoWEY\nKCZVgT78EK65Bn78Edq1MyUEUYXWrFlDnz59mDJlCqNHj/brNU6nk8TERAoKCnA4HFUcoRB++uEH\nuOgiuOUWmD3b7GgqxO+kcPBg0UlrCg6HHaezBE3TAx4ggNcL55xj58orNaZPL6uS9zBDdZRdsFu+\n/BNuu+1fTJs2k/79B/r9OperiAYN6pGd/RdxcaFVPTebHHcV909lp+zfh+PSLuj1Uij8+LOgXMYi\nKSnulNv4XbPRdQNdP3n+0DQdr7fqDrAbbvDw8stRZGSUEnfqzxVSqrrsgtW6dd8xevStLFiwiC5d\nup1WGRz+QkZq2QWClF3FHVd2bjc1bhqM4fVSsPA1dIsNQrRsQ6ZRcfBgDy4XfPBBqLXQiRPRNI1x\n48aQkTGFLl26mR2OEBVnGMRPGIf1p/U4X12MHuKDJUImKaSlGXTrprFoUZTZoYgA+P77dWzduoWJ\nE8eTllaH9PS65fe7d+8yOzwh/BbzykvYF/+XwunP4u1wvtnhVFpI/ey+6SYPI0fa+eUXldatQ7Nq\nJnw6dbqAv/7KNzsMISrF9s0q4jMepPi2O3EPutHscAIiZGoKAD17eklN1Vm4UOYsCCHMpW7/E8fN\nQ/Fc3BXXI1PMDidgQiopWK0wbJiHd9+1kZdndjRCiIjldJI4dBB6zVo4X17oOzmFiZBKCgBDhnjQ\nNFi8WGoLQggT6Dpxt92MumcPzkVvYdRIMjuigAq5pFC7tsG113p59dUoWVJbCFH9HnwQ22efUvjS\nK2hNm5kdTcCFXFIAGDWqjKwslS++kCW1hRDVJ+rVV+Dppyl5bCpll/UwO5wqEZJJoV07nXPP1Xjl\nFRmeKoSoHlErPif2/nthzBjct/u3HEsoCsmkADByZBlff21l69aQ/QhCiBBh/eUnHKOG4bm8B8wK\n/ktqVkbInlH79PGSnKzzyivS4SyEqDrq7l04bhyAt2kzXPMWgiW8m61DNilER/uGp775po18mQMl\nhKgCSqGTxMH9wWaj4LW3CbuF104gZJMCwMiRHnQdXn1V+haEEAFWWopj2GDU3bsoWLwEo25dsyOq\nFiGdFGrXNhgwwMO8eTZKS82ORggRNjQNx+03Y/thHc7X3kJrcZbZEVWbkE4KAHfcUUZOjsKSJdK3\nIIQIAMMgfvy9RC3/GOfLr+LpdKHZEVWrkE8KjRoZXHWVl7lzo9BljTwhRCXFPvU49kULKZw5m7Ke\nV5kdTrUL+aQAMHp0Gdu2qSxfHj7rjwghqp993gvEzZxGUcYU3DcMMTscU4RFUujQQadjRy/PPy8d\nzkKIiol+923iJz5A8e1jKLnzLrPDMU1YJAWAO+8s4/vvLXz3XXiPIRZCBF7Usg9IuPNWSgcOxvXo\nY2E9Oe1UwiYpXH65RtOmGs8/Lx3OQgj/RX3+KY7bRuLucy2FzzwPaticFiskbD69qsKYMWUsX27j\nl1/C5mMJIaqQbdVKHCNvouzynhTOeTnsZyv7I6zOntdf7yU9XWfmTOlbEEL8M1vmtyQOu4GyLt1w\nvrQAbNLKAGGWFKxWuOceNx9/bOO338LqowkhAsj6/Xc4bhyAp0NHnAte862bI4AwSwoA/ft7SUuT\n2oJZVq78krPPbsJtt438x+2mTZtKSkoS6el1SU+vS1paHdLT65KTk1NNkYpIZf1pPYk3XI+3dRsK\n/vsGxMSYHVJQCbukYLPBXXeVsWyZld9/D7uPF9TmzHmWjIwJNG7cxK/tBwy4gZ0797Fz5z6ysvaz\nc+c+kpOTqzhKEcmsP/5AYr8+aE2b4Xx9SUQscHe6wvKsOXCghzPOMJg1S2oL1cluj+Gzz76iYcMz\nzQ5FiONYv/+OxP7XorU4i4K338OITzA7pKAUlkkhKspXW3j/fSubN4flRwxKo0bdSvxpfNE2bvyV\nXr0up3Hj+nTt2olVq1ZWYXQiktn+l0nigOvwtmpN/ptLMRIcZocUtPxeF0JVFVT1+AkdFot6zH2w\nuOkmjdmzDZ56Kpr//tdtdjgnFKxlV1mKoqAoClbryT9X/fr1adSoEY8+OoW6deuxcOErDBkygDVr\n1vnV/HR02f3T+4jjhetxdzLWb78hftD1eM/tgOv1JVgr0WQUCWWnGIZh+LOhYRgoITbLb9EiGDoU\nvvsOzj/f7Ggix4gRI3C73bz++uun9bpOnTrRo0cPJk+efMptnU4niYmJFBQU4HDIrz5xEp98Av36\nQefO8MEHEBtrdkRBz++aQm6u66Q1BYfDjtNZgqYF1zKlV14JZ51l5/77Dd5/P/guuBDMZVcZbreX\nsjIveXmu03pdamp9tm/P8ut1LlcJwKGykwlHpyNcj7u/s737NnG334Lnip645r8KbgPcp3dM/l2o\nl11S0qlrSX4nBV030PWTVyo0TcfrDb5CevDBUoYOjWXFCoWuXTWzwzmhYC27ijIMA8Mw/vEzzZo1\njfPO60jnzl3Kn9u8eTPXXdfPr7I4/IUMt7KrTuFcdjEL5xM3YRzu/oN8S1dYrRDAzxrOZRe+DWOH\n9Oih0aGDxuOPR+NfQ5moKhdd1IF1674DIDc3lwkTxrFt21bcbjdz585mx47tDBw42OQoRUgzDOzP\nziDhgXspuflWCp97wZcQhN/CvrQUBTIy3FxzTSwffWSld2+v2SGFrbS0OiiKgsfjAeCTTz5CURR2\n7twHwLZtf+ByFQGQkTEZRVHo168P+fl5NG/egqVLl1GvXopp8YsQZxjETXmE2OefxTX+IYrHPRDR\nq51WlN8dzQcOFJ7weatVJSkpjrw8V1BXp264wc727SrffOMiKkimL4RK2QWj4uIiGjZMZceOPcTG\nxpsdTkgJy+PO6yV+/D3YX/sPRf9+ipJ/3V4lbxPqZVe79qmHjId989FhGRluduxQWLBAFr0SIqwU\nFeEYOoiYNxfjnPNSlSWESBExSaFlS52hQz1Mnx7NwYNSpRQiHCj79lHj2quw/W8tBYuX4B5wg9kh\nhbyISQoA48eXAfD000HSfiSEqDDLls0kXdUddf8+8j9cjueS7maHFBYiKikkJxuMG+fmP/+xsWlT\nRH10IcKKbe0aavS6HCM+nvxPV6C1am12SGEj4s6Mo0Z5aNjQICNDhqgKEYqily4hsf81eNu0JX/Z\nZ+hn1Dc7pLAScUkhKgomTy7lm2+sfP65zIQVImToOrFPTMFx2yjc1/aj4I13MByJZkcVdiIuKQBc\ncYVGt25eJk6MobjY7GiEEKeiFBXiGH4jsc/OoChjCoWzXyRoxpaHmYhMCooCTz5Zyr59ilyhTYgg\np+7cQY1el2P79huci96kZMzdMimtCkVkUgBo1Mjg7rvLmDs3SjqdhQhStsxvSerRDaW4mPxPvqTs\niivNDinsRfTZ8M47y2jYUOf++6PRQ29yohBhLeY/C0i8vg/elq3I++wrtBZnmR1SRIjopBAdDU8/\n7WbdOitvvCEznYUICqWlxN93Nwn3303pTcMpeOs9jJq1zI4qYkR0UgDo3FljwAAPU6ZEk5Mj7ZRC\nmEnN2kmN3j2IefM1Cmc8R9FTM8EmP9iqU8QnBYBJk9wYBjz6aLTZoQgRsaK+/Iykyy5Gzcsl/+Mv\nKL1puNkhRSRJCvhmOk+eXMqSJTaZuyBEddM0Yp98nMTB/fGc15G8L77Ge047s6OKWJIUDhk0yMtl\nl3m5994YcnPNjkaIyKDk5JA4qC+xz0zH9dAjOBe9hZFU0+ywIpokhUMUBWbOLKWsTOGhh2LMDkeI\nsGdbs5qk7p2xbvyFgrffp/ju+0CVU5LZ5H/gKPXqGUydWsrSpTaWLQv7i9IJYQ6Ph9gnppDY92q0\nRo3JW/Etni7dzI5KHCJJ4W/69vXSq5eH8eOj2b9fRiMJEUjqzh3U6NOT2NmzKH4wg4J3PkRPSTU7\nLHEUSQp/oyi+uQuKAvffLyupChEo0UuXkHRpZ9QDB8hf9pmvucgiAzuCjSSFE6hd2+Dpp918+qmN\nt96SZiQhKkMpyCfhzltx3DaKssuvIG/larwdzjc7LHESkhRO4uqrvQwa5GHChBj++EOakYSoCNvK\nL0nqegFRn36M87kXKHzhFVnuOshJUvgHTzxRSmqqzi232HG7zY4mPK1c+SVnn92E224baXYoIoCU\nokLix91FjUF90Zo0I+/rtbgH3Sirm4YASQr/ID4eXnqplC1bVKZMkdnOgTZnzrNkZEygceMmZoci\nAsi2ZjVJ3S4k5t23KXxqJgVL3kev38DssISfJCmcQuvWOpMnu5k3L4pPPpH+hUCy22P47LOvaNjw\nTLNDEYFQXEzcxPHUuK4X2hn1yV2VSemIm6V2EGLkLOeHkSM9fPuthTFjYmje3EXjxjIkKRBGjbrV\n7BBEgNgyvyX+3jFY9uym6LGplPzrdpmIFqL8TgqqqqCqx2d8i0U95j5cPf98GZdfbmH48Fi++KKE\n+PjK7zNSyu5UFEVBURSsVv/L4eiyO53XicAed0peLvZHHiZ68X/xnt8J5xtL0Js2C9tfm5HwnfX7\n/65mzTiUf6gGOhz2gAQUrJKS4IMP4Pzz4b774njzzcDVisO97E4lOtoKaCQlxfn9GotFA3xl53D4\n/zpxRKWOO8OAN96Ae+6B0lJ44QWst9xCYoTUDsL5O+t3UsjNdZ20puBw2HE6S9C08L58Wb16MHu2\nhREjYmjd2s3o0d5K7S+Syu6fuN1eysq85OW5/H6Ny1UCcKjsZALU6ajscafu3EHsuLuxrfySsmuu\no3jqNIx69aCgpAqiDS6h/p3154eX30lB1w10/eRt6Zqm4/WGXiGdrl69dEaPVpk0KYqWLTUuvlir\n9D4jpexOxjAMDMM4rTI4/IWM9LKrjNMuO7cb+0tziZvxJHrNWhS89taRayZH2P9BOB93kVHXC7CJ\nE91cfLHGyJF2tm6VIhThz7byC5K6diJu6hRKho4kd/W6IwlBhJVw7Q+qUlYrzJ9fwtVXxzJ4sJ1P\nPy0mOVlGJJ2utLQ6KIqCx+MB4JNPPkJRFHbu3GdyZOIwdcd24h95kOjln1DWuQvOhYvRzmppdlii\nCklSqCCHA157rYSePWMZNszOu+8WEyOXYTgtWVn7zQ5BnExxMbHPzSD2+efQayXjnPcq7j7XyZyD\nCCBtH5WQlmawaFEJv/yicvfdMbKiqgh9uk70kjepeVEHYuc8S/EdY8hd8wPua/pKQogQkhQq6dxz\ndebM8V2Y56mnoswOR4gKs61ZTY0ruuEYfQvetu3J/eY7ih98BOJkyG8kkaQQAH36eMnIcDNzZjTz\n59vMDkeI02LZugXHTQOpcV0vsFrI+/AznAtfQ2/U2OzQhAmkTyFA7ryzjJwc3/Wda9QwuP76ys1h\nEKKqKX/9Rfy0J4n570L0M+rjfHmhNBMJSQqBoigwaZKb3FyFsWNjqFGjhMsuq/wcBiECTck9CFMn\nkzhnDkZUNK6HJ1Ny860QLSsBC0kKAaUoMGtWKQUFMYwYYec//ynh0kslMYjgoDgLsL8wh9iX5gIG\npaPH4rp1NEZiDbNDE0FE+hQCzGqFefNK6dJFY9gwOytXyhIMwmQuF/bnZlKzQ2tin38W9/CRsH07\npQ9lSEIQx5GkUAWio2HBgpLyxPDVV5IYhAlcLuwvzqHWeW2Ie+rfuPv2J/f7nymZ8m9ITjY7OhGk\nJClUkaMTw9ChkhhE9VGcBcQ+M51aHVoRNzkD9xU9yf3feoqenIFet57Z4YkgJ0mhCh1ODBdfLDUG\nUfWUnBxin5hCzXZnEzvjKdy9ryX3uw0UPfM8eoM0s8MTIUI6mqtYdDQsXFjCyJF2brrJzosvlnL1\n1TJcVQSOuncP9rnPYV/0KqBQMmwkJXeMkVqBqBCpKVSDw4mhVy8vN98cw+LFMsFNVJ7ll59JGH0L\nNTu0JubN1ym+7U4O/vgrrsn/loQgKkxqCtUkKgrmzi0lMTGae+6J4eBBhXvukRqDOE26TtTKL7C/\nMIeo1V+j1W+A6+HJlA4ZipHgMDs6EQYkKVQjiwWeespNrVoGjz8eTVaWyvz5ZkclQkJJCTHvvIX9\nxTlYt26L1/rVAAAWFklEQVTB0669bwby1df4xkELESByNFUzRYEHHigjPV3n3ntj+OsvePlliI01\nOzIRjNQd27H/dyExr/8XJS+Psp69KJwxG2/HTrIchagSkhRMMmiQlwYNShk+3M5VV9lZtKiYBg1k\n7W0BaBpRK78gZuF8olZ8gZHgoPSGGykZ8S9ZpE5UOeloNlHXrjqZmVBUBFdcEcvq1TJkNZIpOTm+\nmcfnn0PijQNQ9++naNYcDv68GddjT0pCENVCkoLJWraElStLaNVKp39/O88/b5OL9UQSw8CW+S0J\nt99MrbYtiJs2Fc+FnclbvpL8L76mdPBN0rYoqpU0HwWBmjXhzTdLmDo1ismTY9iwwcKsWaXEx5sd\nmagq6p7dxLz1OjFvvIZlx3a0hmfieuhRSgcNxqhZy+zwRASTpBAkLBZ4+OEyzjlHZ+zYGK66KpZX\nXimlaVPd7NBEoLjdRH32CfbXF2FbtRJiYnD3vpbCZ+fi6XShdByLoCDNR0Gmd28vy5cXo2lw2WWx\nLFokzUkhTdex/S+T+PvuplabZiTePAzF6aRo+rMc/GULhbNfxHPBRZIQRNCQpBCEmjfX+fzzYq6/\n3sO4cTGMHBlDXp7ZUfln165sbryxPy1aNKRDh9Y89tijJ9xu2rSppKQkkZ5el/T0uqSl1SE9vS45\nOTnVHHHVsPy2kbjHHqVmh9bU6NOTqJVfUDpkOLmr15H/yZeUDhkmk81EUJLmoyAVFwczZri55BKN\ne++NoVu3OObOLeWii4L7oj0jRgyhbdv2vPTSQg4c2M/gwddTp04dbr119HHbDhhwA88+O9eEKKuG\nuiub6KVLiHl3CdZNG9GTknD36UtpvwF4z+8IqvwGE8FPkkKQu/pqL+3buxg9Ooa+fe2MGuXhoYfc\nQdkJvWHDj/z2268sXbqM+Ph44uPjue22O5k374UTJoVwoGZnEf3xh0R/9CG2df/DsNtx97wK18RH\nKOvW3be+iRAhRH66hIDUVIN33inhscfcLF5so1u3OL75JvjmNPz88080aJBGwlHNIm3anMMff2zF\n5XIdt/3Gjb/Sq9flNG5cn65dO7Fq1crqDLfCLNu2Yn92BjUu70qtc1sR9+/J6ElJOOe8xMGNf1D4\n0kLKrrhSEoIISX7XFFRVQVWP7wyzWNRj7oX/TqfsrFa44w6NK68s4a67orn++liGDvUwaVIZNYLk\niooFBXnUqJGE1Xrk8yQn+4ZXOp15JCYmlD9fv359GjVqxKOPTqFu3XosXPgKQ4YMYM2adTRu3OSU\n73V02R39flXCMFA3bSTqww+IWvYBlk2/YcTG4rnsCoruHIvnip6Q4PtswZeqjyff2YqLhLJTDMO/\nsS2GYaDICImgoOswbx7cfz/ExMCTT8Lw4eY3WU+dOpX33nuPdevWlT+3bds2mjVrxp9//kl6evo/\nvr5Tp0706NGDyZMnn/K9nE4niYmJFBQU4HBUQYdtWRmsXg3LlsFHH8G2beBwQO/e0K8f9Oghk8pE\nWPK7ppCb6zppTcHhsON0lqBpMqb+dFSm7AYMgC5dFB59NIpRo6zMnasxbVoZbdua938QG+vgwIEc\n8vKONBVt374LRVGwWmOPef5EUlPrs3171im3A3C5SgAOlV1gfp8rOQewffE5ts8+xbZyBUpRIXrq\nGXh69KTsiafxdunmuzgGgNsA96njDEbyna24UC+7pKS4U27jd1LQdQNdP3mlQtN0vN7QK6RgUNGy\nS06G558vYcgQCxMmRNO9ewxDh3p44IEykpOrf3JD69Zt2bUrmwMHckhKqgnADz98T7NmLYiKijnm\nM86aNY3zzutI585dyp/bvHkz113Xz6+yOPyFrNRx5/Vi3fAjUatWErXyS6z/9z2KYeBpfy7Fd96F\n+/KeaK1aHzuHIIyOcfnOVlw4l134NoxFkAsu0FixopjHHnOzdKmN88+PY+bMKE7Qt1ulWrduQ9u2\n7Xn88UkUFRWydesWXnzxeUaMuBmACy88l3XrvgMgNzeXCRPGsW3bVtxuN3PnzmbHju0MHDi4SmNU\nd+4g5j8LcIwYQq2zGpF01WXYX3wevW49Cp95npxf/yB/+VcU3zserXUbmVQmIo4MSQ0TVivccouH\nfv28PPNMFDNmRLFwoY377y9j8GBPtV2HZcGCRdx77xhatWpKQoKD4cNHMXz4KAD+/HMbLlcRABkZ\nk1EUhX79+pCfn0fz5i1YunQZ9eqlBDQepdCJ7dvVRK1agW3VSqzb/8SwWPCeex4lt95BWbdL8bZt\nLxeqEeIQvzuaDxwoPOHzVqtKUlIceXmusK1OVZWqLLudOxWmTo1m6VIbZ56pc889bvr182ILk8tD\nFxcX0bBhKjt27CE29sikDaUgH9v/1mLL/Bbb2m+x/vwTiq6jNTyTsm6XUtatO57OF2M4Ek2M3lzy\nna24UC+72rUTTrmN/DwKU+npBi++WMro0WXMmBHF2LF2pk/XufvuMgYM8ITNEHolP4+or7/Glrka\nW+YarL/+jGIYaKln4LmwM6VDR1LWuQt6wzPNDlWIkCBJIcy1bq3z6qulbNxYxqxZUYwbF83MmVGM\nGVPGoEEe7HazIzwNhoG6KxvbD+sw1n4LQGLbs0kEtAZpeC64iJKbb8VzwUXo6Q2lP0CICpCkECHO\nPltn/vxSfv9d5ZlnopgwIZqnnopi6FAPI0d6qFcvCJdiLS3F+tMGbD+sw/bDOqw/rMOy7y8AtPSG\nABRPfwZvt+7oaf88B0II4R/pUzCRmWW3Y4fCK69EsXixjdJSuOYaL7feauI8B8NA3bEd20/rsR5O\nAr/8jOLx+GYPt22Pt8P5eDqcj+fc83DF2U/YpyBOTb6zFRfqZSd9CuKkGjY0eOwxN+PHu3n9dRvz\n5kXx7rtxdOigMXRoGX36eKtuwq5hoO7cgfXnDdg2rMf60wasP29ALcgHQEtriOe88yntfwPe887H\n27LV8aODiouqKDghIpvUFEwUTGWnabB8uZVXX7Xx9ddWHA6D/v093Hijh1atKhGbrqNm7cT6y0/Y\nftqAdcN6rD+vR80/lADOqI+3TVu857TF07Yd3jbtMJKTT7nbk40+EqcWTMddqAn1spOagvCbxQK9\nennp1cvL9u0Kr79u4/XXbbzyShRnnaVx/fVe+vXzkJp68t8QSqETy6ZNWDf+gvW3jVh/+xXLpt9Q\ni3w/KLSUVLzntKPk1tF427bD06YdRu3a1fURhRB+kJqCiYK97Dwe+OorC++8Y2P5cituN3TurNH/\nulKubr6Zmns3Yt30G9aNv2L9bSOWrB0AGFYrWtPmeFuejbdlK7xnn4337DYYdesGLDapKVRcsB93\nwSzUy05qCqJSbO4irqqzhd6XbcZVL4tla5J54/uO3LW6I+NoTXf2c228h6tb20i6uk95EtCaNjuy\ncJwQIqRIUoh0hoH6114s2//Esu0PLFs2Y93yO5atW7Dsyi7fLC71DIY3a86QC7aQXWcjH+y/iA83\ndOOOH3pyx3dwgapx1Rleusd6aRQdhMNbhRB+kaQQCQ6f+P/c5jv5H32/cztKcbFvM1VFS2+I1rwF\n7uuux9u0GVrzFmhNm2HEH6l21gRGACPwkpPjYvlyKx9/bGXSpGgmTozhzDN1unf30r27lwsv1EJr\ngpwQEU6SQphQigpRs7KwZGehZu/EkpWFJWun7+S/40+UEt/1BwxFQW+QhtawEZ6OnSgddCNao8Zo\nZzbyTQg7zWaf5GSDIUM8DBnioagIVq+2smKFheXLrcyfH0VMjMGFF2p07epLEK1a6VhC4fJkQkQo\nSQoh4oQn/ews1OwsLNk7UfPyyrc1oqPR6jdAb5CGp9MFlN5w+MTfGC0tvcra++Pj4corvVx5pRfD\ncLNli8rKlRZWrLAydWo0paUKDodBx44aF17oSxKtW+uyQKkQQUS+jmYzDJSCfCzZu1H37Maydw/q\nUTfLnj2oe3cfe9KPikJrkIbeIA3vOe1w977G9+u/QRp6Wjp67TqmX5tTUaB5c53mzXVuv92D2w3r\n11tYu9bCmjUWpk2LprhYIT7eoEMHjfbtNc49V6NdO92UCwQJIXxkSGpVcrtRcw6g7t+H+tdfvpP8\n3j2oe3b72vgPPT76ajiGomAk10ZLPQM9JRU9JQU9JRWtQRpag3T0tDT0OnVNP+lXVlkZ/PSTSmam\nlXXrLKxfr5KT4/tMaWk67dtrtGun0b69TuvW2nGzq2VIasXJd7biQr3sZEhqVfB4fCf6A/tR9+9D\nOXDopH/osXr040Ozdg8zrFbfib5eiu+k37oNlsYNKaqRjKfuoQRQtx5hs671P4iKgvPO0znvvDIA\nDAOysxV+/NFy6KayfLmvyUlVDRo10mnZ8vBNo3FjWQFViKoQ2TUFw4DiYtTcg6i5B1EOHkTNyz3y\nd+5BlLw81IM5vpP8gf2oubnH7UZPSkKvUxe9dh30OnV897XroNepi1G7NnrtOmj1Un3LNxz1Cz+k\ny64aeDzw++8qP/1kYdMmlY0bVX77zUJ+vgI4gUQ6dMjl7LPjaNJEL7/Vr29IZ/Y/kOOu4kK97CKr\nplBSguosQMnPRykoQC3IQykoQMnPQy0/weeiHjx00j+cCNzu43ZlxMSg16yFXrMWRs1a6MnJeM9q\neegk/7cTf3LtiPhlbwabzXc9iNatj3z5DAP27lVYv76UESMgJcVg3ToLb71lo7TUV3uIjvbVLBo3\n9iWJxo11GjY0SEvTqVvXCPWWNyGqVHAkBcOAkhKUwkJUVyFKYSGK04mSn49a4DvJKwWH/z500s/P\nR3EeuT/RyR18nbJ6rWTfyb1mTfSatdCaNPGd9GvVwkjyPWfUqlWeCKpueVBRWYoCqakGNWpoAMye\n7SY2thhdh127FP74Q2XbNpU//vDd3nrLxt69R7JAVJRB/fq+BNGggU56ukGDBvqhxwbJyZI0RGSr\nfFIwDNi2DcuufSj5BShFh07qRUWH7gvLn1OLilAKnb5/K/rbNpp20rfQExwYNWpgOBLRD917mzXH\ncCRi1KiBnpiIkVgDIzERPTEJI/HQ845EsNvlClwRQFUhLc0gLU3j0kuPPZZcLsjOVsnKUsjOVtm5\nUyU7W2HDBgvLlqmHmqN8bDaDunUN6tUzSEnRSUkxqFfPd5+SYlC3rm90lMMhh5UIT5VOCtEvPA8P\nT8Dxt+cNVcWIT8BIOHSLi/fdxyegpaRixB96HHfUNoe3j49HdyT6Tu6ORKSBWFRGXBy0aKHTogXA\n8T8+CgogK0slO1tl716FffsU9u71/b15s8revTYKC4/NADabQa1avprF4VutWga1axskJ+vlj5OT\nDZKSDBISJImI0FDppOAePpLYizriNKx47bHo8Q6M+HhfE4x8C0QISEw8vu/i74qK4K+/FPbvV8nJ\nUThwQOHgQYWcHN9t1y6FDRtUDh5Uycs7/rhXVYPEREhMNKhRw8Dh8N0ffnz0vyUmGiQkGMTHQ3y8\nQVyc72+Z5CeqQ+UPs9hY6NYNLc+FFoK98UL4Iz4emjQxaNLk5M2ch3k8kJvrSxw5OQr5+b6b06mQ\nnw8FBb7HeXkKO3f6mq8KChQKCsAwTv5DKibGOJQkfMni73/Hx4PdbhATc+y93e57re8e4uMV6tWD\nsjIFm00p/zepkAsIlo5mIcKIzQZ16/r6Jk6HrvtqJPn5CoWFCkVFCi4XuFwKRUVQVKQcuh1+Tinf\nfvdulaIiKClRKC09cu/x/FNt/dgBFVFRBtHRvtFbNptvUF1UlHHo3tdkFh3NoX878vyJ/j769Tab\nr5Zjs/mGClutvpvvb+Oovw8/b/xtG992p3qtxRLyczqDgl9JwTAMXK5ClBM0B1ksKhaLhstVgqZJ\nTeF0SNlVXHGx65j7cGG1QnKy7xYImgYlJeB2H0kUZWUqqhrDwYNuXC6D0lIoLVUO3YPXq+B2+2o8\nZWUKHo/v9V4vuN0KZWW+GelFRcqhbTi0jXJom2NfW1bm22d1URTfCLLDN4vF15J99HO+m287RTmS\nUA7fjt3eOOrfFaxWJ6CddJ+KwilvvjiP3A6/5/G3I6Ph/HltXJzB6NEeEk4yHcHpNEhISDjhuby8\n/PyZvOZ0OklMTDyd/xchhBBBqKCgAIfj70ODjvArKRiGwc6de09aU3A47Did8mv3dEnZVVxxsYuW\nLZvy229biY2NMzuckCLHXcWFetklJcWdsqbgV/ORoijExZ24PmK1qjgccWiaJSSnfZtJyq7yYmPj\nZEG80yTHXcWFetk5HKde5kK6ZYQQQpTze0E8IYLJ4X6uU7WPCiFOjyQFEZIMw6CwsPCU7aNCiNMj\nSUEIIUQ56VMQQghRTpKCEEKIcpIUhBBClJOkIIQQopwkBSGEEOUCmhRcLhdpaWmMHDkykLsNa0uX\nLqVt27YkJCRw1llnMX/+fLNDEmEsKyuLvn37kpycTEpKCiNGjMDpdJodVsi55557UMN0SdaAfqpH\nHnmEoqKiQO4yrH3//fcMGTKExx9/nIKCAmbOnMno0aPJzMw0OzQRpnr37k3NmjXJzs7m//7v/9i4\ncSP33Xef2WGFlA0bNrBo0aKwnR8TsKTw888/8+abbzJ8+PBA7TLs5ebmMnHiRK6++mpUVeXKK6+k\nTZs2fPPNN2aHJsJQQUEB5513HlOnTsVut5OamsqwYcPkeDsNhmFw++23M27cOLNDqTIBSwq33347\nTzzxhCyxfRp69OjBxIkTyx9rmsbevXs544wzTIxKhKvExETmz59P7dq1y5/LysqS4+00vPjii9jt\ndgYPHmx2KFUmIEnhpZdewmKxMGzYsEDsLmKNHz+e+Ph4Bg4caHYoIgL88MMPzJkzh4cfftjsUELC\nvn37mDRpEi+88ILZoVSpSl+Oc//+/Tz66KOsXLkyEPFErAceeIC33nqLVatWERUVZXY4IsytWbOG\nPn368PTTT3PJJZeYHU5IGDduHKNGjaJ58+bs3LnT7HCqzGnXFF577TXsdjuxsbHExsZy3333MWzY\nMFq2bFkV8YWVw2V3uPzA10Y5bNgwPvroIzIzM2nSpInJUYpwt2zZMnr16sVzzz3H6NGjzQ4nJKxY\nsYLMzEwyMjIA3/c2XFV6QTxVVUlKSiofnlVcXIyu6yQkJLB///6ABBnOxo4dy//+9z+++OIL6Y8R\nVS4zM5PevXvz9ttv0717d7PDCRkjR47k7bffxm63A6DrOnl5eSQnJzNnzhwGDBhgcoSBU+mksGfP\nnmMez5gxg927dzNr1ixSUlIqFVy4O1yF//3334/p/BOiKmiaRps2bbjnnnu4+eabzQ4npBQUFOBy\nucofZ2dnc8EFF7B7926SkpKIiYkxMbrAqnSfQmpq6jGPHQ4HeXl5khD8sHDhQpxOJ+np6cc836VL\nF5YvX25SVCJcrV27lt9//52xY8cyZswYFEXBMAwURWHz5s00aNDA7BCDVmJi4jE1eY/Hg6IoYXme\nk+spCCGEKBee87SFEEJUiCQFIYQQ5SQpCCGEKCdJQQghRDlJCkIIIcpJUhBCCFFOkoIQQohykhSE\nEEKUk6QghBCinCQFIYQQ5SQpCCGEKPf/RkKRDjcr/0IAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# fig3.4\n", "p1 = plot(lambda x: exp(-1+0.4*x), [x, -4, 5], rgbcolor='red')\n", "p2 = plot(lambda x: exp(-2-0.8*x), [x, -4, 5], rgbcolor='blue')\n", "(p1 + p2).show(figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

線形予測子とリンク関数

\n", "\t

\n", "\t\tglmで使われるリンク関数と線形予測子の関係を整理すると以下の様になります。\n", "\t\t線形予測子は、パラメータが線形結合しているもので、ターゲットとなる変数\n", "\t\t（ここでは$\\lambda_i$）を線形予測子で表したときに、以下の関係がある場合、\n", "$$\n", "\t\t\\lambda_i = f(線形予測子)\n", "$$\t\t\n", "\t\tリンク関数は、fの逆関数を使って以下の関係にあります。\n", "$$\n", "\t\tf^{-1}(\\lambda_i) = (線形予測子)\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tですから指数関数で表された$\\lambda_i$のリンク関数は、対数を使って以下の様になります。\n", "$$\n", "\t\tlog \\lambda_i = \\beta_1 + \\beta_2 x_i\n", "$$\t\t\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

ポアソン回帰の実行

\n", "\t

\n", "\t\t久保本では、モデルから推定される予測の良さを重視しており、その指標として対数尤度が最大になることを目安としています。\n", "\t

\n", "\t

\n", "\t\t今回のデータに対する対数尤度(Log-Likelihood)は、以下の様になります。この対数尤度の最大値を追求することになります。\n", "$$\n", "\t\tlog L(\\beta_1, \\beta_2) = \\sum_i log \\frac{\\lambda_i^{y_i} exp(\\lambda_i)}{y_i !}\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tここでは、Pythonのstatsmodelsパッケージを使ってGLM（一般化線型モデル）を使ってポアソン回帰を行います。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 15, "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": "markdown", "metadata": {}, "source": [ "\n", "\t

glmの引数

\n", "\t

\n", "\t\tglmの引数について説明すると、\n", "\t\t

fit: 回帰の結果をセットする変数
y ~ x: モデル式を線形予測子で表現します
data: dでデータフレームを指定します
family: sm.families.Poissonで確率分布としてポアソン分布を指定します
link: sm.families.links.logで対数を指定します
fit(): 最後にfit関数で回帰を計算します

\n", "\t

\n", "\t\t回帰の結果は、summary2関数で表示します。結果は、\n", "\t\t

対数尤度(Log-Likelihood): -235.39
逸脱度(deviance): 84.993
切片$\\beta_1$(Intercept): 1.2917、標準誤差は0.3637
xの係数$\\beta_2$: 0.0757、標準誤差は0.0356

\n", "\t\tと求まります。\n", "\t

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

Model:	GLM	AIC:	474.7725
Link Function:	log	BIC:	-366.3137
Dependent Variable:	y	Log-Likelihood:	-235.39
Date:	2016-07-17 01:38	LL-Null:	-237.64
No. Observations:	100	Deviance:	84.993
Df Model:	1	Pearson chi2:	83.8
Df Residuals:	98	Scale:	1.0000
Method:	IRLS

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

	Coef.	Std.Err.	z	P>\|z\|	[0.025	0.975]
Intercept	1.2917	0.3637	3.5517	0.0004	0.5789	2.0045
x	0.0757	0.0356	2.1251	0.0336	0.0059	0.1454

" ], "text/plain": [ "\n", "\"\"\"\n", " Results: Generalized linear model\n", "==============================================================\n", "Model: GLM AIC: 474.7725 \n", "Link Function: log BIC: -366.3137\n", "Dependent Variable: y Log-Likelihood: -235.39 \n", "Date: 2016-07-17 01:38 LL-Null: -237.64 \n", "No. Observations: 100 Deviance: 84.993 \n", "Df Model: 1 Pearson chi2: 83.8 \n", "Df Residuals: 98 Scale: 1.0000 \n", "Method: IRLS \n", "---------------------------------------------------------------\n", " Coef. Std.Err. z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------\n", "Intercept 1.2917 0.3637 3.5517 0.0004 0.5789 2.0045\n", "x 0.0757 0.0356 2.1251 0.0336 0.0059 0.1454\n", "==============================================================\n", "\n", "\"\"\"" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fit = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "fit.summary2()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "84.9929964907 -235.38625077 474.77250154\n" ] } ], "source": [ "# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットされている\n", "print fit.deviance, fit.llf, fit.aic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

回帰結果のプロット

\n", "\t

\n", "\t\tλの予測値は、predict関数で取得することができます（fitのnuにセットされているみたいです）。\n", "\t

\n", "\t

\n", "\t\t施肥別の分布図に予測値を重ね書きします。\n", "\t\t予測値をyhatにセットし、pandasでソートした後、lineplotします。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# λの予測値nuをyhatにセット\n", "d['yhat'] = fit.predict()" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAESCAYAAAD67L7dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8XHWd//HXpGluzW2STC+hZVoEvspiBXwgW7ErIA+w\nroqIuD4W0AJWuemy3UuFZUEB6Zrl0h+y1v1x68pjXdTlpmjBBbuKC0ildMEf8OXWxsAATTqTS5t7\nMr8/ziSdpJkkk8w5ZzLn/Xw8+shk5sz3fL7fc/rpt9855zOhZDKJiIgER5HfAYiIiLeU+EVEAkaJ\nX0QkYJT4RUQCRolfRCRglPhFRAKm2O0dGGOOBh4EbrbWfs8YUwz8G3A40Al8zlrb4XYcIiLicHXG\nb4ypAG4FHkt7eh2wx1p7AvAjYLWbMYiIyFhuz/h7gTXAN9Ke+xRwNYC19g6X9y8iIuO4mvittcNA\nnzEm/enlwCeMMf8MvA1cYq1tdzMOERE5wI8Pd0PAS9bak4H/B1zpQwwiIoHl+oe7E3gH+E3q8aPA\nNyfbeHBwKFlcPM/tmERECk0o0wt+JP6tOOv+W4APAnayjROJbg9C8lYkUkVra5ffYfhG/Vf/g9x/\n8GYMIpGqjK+5mviNMccBNwFRYMAY8zngL4FbjTEXAl3Al9yMQURExnL7w90dwMkTvPR5N/crIiKZ\n6c5dEZGAUeIXEQkYJX4RkYBR4hcRCRglfhGRgFHiz4HnnnuWq67aMO3tH330URejERGZnBJ/joRC\nGW+SG2NgYIC7777b5WhERDJT4p+Br3xlLbHYWwC0tu7htts20dPTzXXX/SNf+tIX2LLFKTr6+98/\nw0UXXcDXvvZVrrzy7xgcHOS7372ZV199lZtv/o6fXRCRAPOjZENO/fhXr7H95T05bfP49y7k86cc\nnvH1j3/8z3n88V9y3nnn89vf/oZTTz2d++77ET/84X0MDQ1x9tmfZu3aL9PV1ck3v/ltFi9ewvXX\nX8MzzzzNX/7lF3n11ZdZv376S0Pij454nCc2rKemeTcd0Sirm26hJlznd1gis6YZ/wyceupp/OY3\n2wB48sknaGho4Mgj30tJSQnl5eWj29XWhtm48Touu+wrPPfcs3R0qPr0XPLEhvWsfeh+PrtzB2sf\neoAn/n693yGJ5MScn/F//pTDJ52du6G6uoZIZBEvv/wiyWSSSGQh8+YdXEF048ZrufHGWzn00Ci3\n3NLkaYwyezXNu0fLG4ZSv4sUAs34Z+j009dw003f4aSTPkYymZxwm+7u/SxatIiuri527Pg9g4OD\nhEIhBgcHPY5WZqIjGmXkyCaBjuhyH6MRyR0l/hk68cQ/IxZ7k5NP/ljGbc4882wuuugCbrzxBs45\n50vcc8/dFBUVMTAwwNVXX+FhtDITq5tuYcsZn+X+Y45jyxmfZXXTzX6HJJIToUyz1XzR2tqVlwHu\n2PF7Hnnk51x55TVZvzfo9cjVf/U/yP0Hz+rx59UXscx5d975r2zf/juuv17r9iIy9yjxz8CFF36V\nCy/8qt9hiIjMiNb4RUQCRolfRCRglPhFRALG9cRvjDnaGPOaMeaScc+fbowZdnv/IiIylqsf7hpj\nKoBbgcfGPV8KfAOIubl/N9122yasfYl4fC89PT0sXbqM6upqXekjInnP7at6eoE1OEk+3ZXAbcA/\nu7x/11x22eUAbN36MLt2vc4ll/yVzxGJiEyPq0s91tpha21f+nPGmCOBldba+4DpFbHPgf/d9jj/\n9d1bePmZp73apUjWOuJxHl63lidOO4mH132JjkTc75CkAPlxHf/NwNe83OF/b76N9zd9m1P37+e5\n2jDPfLuJD539F16GIDItIxVBQ0By5w62EOKTt2/xOywpMJ4mfmNMI2CAfzfGhIAlxpht1tqTM70n\nHK6guPjgypfZmP/gT3jv/v0AHNueoPm+/yByyZdn1eaIqqoyystLiESqsnpfttsXGvV/4v43xFrG\nVARtiLUU5FgVYp+y5ecYeJn4Q9baGHDEyBPGmF2TJX2ARKJ71jvuHxp78VDv4HDO6mR0dfXS09Of\nVXtBr1Wi/mfuf1vjUpJsd2b8QFvjsoIbq6Aff/CsVk/G19y+quc44CYgCgwYY84CPmutHflGEk8K\nsJWcu5ad11/DBzo7eSoSofaLF3ixW5GsrW66hS2EUt/6tVwVQcUVriZ+a+0OIOOM3lp7mJv7H3Hi\n2gux71/JD3c8yxEnruaDR/2JF7sVyVpNuE5r+uK6wBRpMx88HvPB43Pe7po1n8x5myIiblLJBhGR\ngFHiFxEJGCV+EZGAUeIXEQkYJX4RkYBR4hcRCZjAXM7phjffbOHWW2+ivb2d4eFhjj56JZde+lfM\nnz/f79ACpSMe54kN61M3PUVZ3XQLNeE6v8OSgJlL56ES/wwNDw/zD//w96xfv4EPfOAYADZtupEt\nW+5g3bqLfY4uWFTYTPLBXDoPA5P4t217jj/84R1OOCHKhz501Kzb2779dyxfvnw06QNccsnXKSrS\n6pnXapp3jylsVtO828doJKjm0nkYiCy1efPjXHDBIq677vOce24JP/nJU7Nus7l5N4cfbsY8V1JS\nQnFxYP4tzRsd0eho0ack0BFd7mM0ElRz6TwMRJb6z/8cYv9+J0m3tx/Lvff+J2efPbs2Q6EQw8ND\nOYhOZkuFzSQfzKXzMBCJPxRKTvr7TESjy7nvvh+NeW5gYICWlj9y2GHvmXX7Mn0qbCb5YC6dh4FY\n6jn33DKqq3cCSSKRp/jiF2f/Sfvxx5/Au+++y5NP/hZwPuzdvPlWtm17bIp3ioj4K5RMelISf8Za\nW7tyEuCzz77Ejh27OfFEw1FH5aYadDy+l+9853ri8b0UF8/n+ONP4IILvjLl+4L+RRTqv/of5P6D\nZ1/EkvE7zQOT+PNJ0E989V/9D3L/wf/EH4ilHhEROUCJX0QkYJT4RUQCxvXLOY0xRwMPAjdba79n\njFkG3AXMB/qBc621e9yOQ0REHK7O+I0xFcCtQPo1jtcB37fWnoTzD8LfuBmDiIiM5faMvxdYA3wj\n7bmLU88DtALHuhyD5Kl4vJ0NG7bR3FxNNNpBU9MphMO1foclOTa+auWZd92B8x/+2fPzHJrL56+r\nid9aOwz0GWPSn+sBMMYUAZcC33IzBslfGzZs46GHzgNC7NyZBO7h9tvP9DssybHxVSvvvXg+p952\nR07a9vMcmsvnry8lG1JJ/x7gcWvttsm2DYcrKC6e501gHopEqvwOwVeRSBWxWBjS6hnGYuHAjEtQ\n+gnQEGsZU7WycteunPXfz3Notvv28xzwq1bP3YC11l431YaJRLcH4Xgr6DewjPS/sTGOU8cwBCRp\nbEwEYlyCdvzbGpeSZHvqKMO+FSty1n8/z6HZ7NujG7gyvuZ54jfGnAP0WWuv9Xrfkl+amk4B7kmt\nkXbS1HSy3yGJC8ZXrTxz82YGclTY1s9zaC6fv66WbDDGHAfcBESBAeAtYCHOh7tdOP9cvmitvSxT\nGyrZUHjUf/U/yP0H/0s2uP3h7g5g7vwzKCISALpzV0QkYJT4RUQCRolfRCRglPhFRAJGiV9EJGCU\n+EVEAkaJX0QkYPwq2SDimni8ncsvf5Snny4C2li1qpJNmz41q8qJ4ytMrm66hZpwXe6CdoFb1SM7\n4nG2/vXXufepMDFW8r5VVWza9HHPKlPOxWORb5T4peBs2LCNRx65kJEaKlu3/gclJdtmVTlxfIXJ\nLYT45O1bchSxO9yqHvnEhvU8sRWe5G4gxO6tSUpKvKtMORePRb7RUo8UnObmahhTD7Iq9dzM1TTv\nHtNiTfPuWbXnhfHjMNsxGFHTvJvdrHSl7enuf64di3yjxC8FJxrtwCkDRepnF9Fo56za7IhGx7TY\nEV0+q/a8MH4cZjsGIzqiUZbzvCttT3f/c+1Y5Bst9UjBaWo6hf7+O3nqqSJgL6tWVdLU9MlZtTm+\nwuTqpptzE6yL3KoeubrpFvb1f423njqfGCs5alU1TU2n56Tt6e5/rh2LfONqdc5cUHXOwqP+q/9B\n7j/4X51TSz0iIgGjxC8iEjBK/CIiAaPELyISMEr8IiIBo8QvIhIwrl/Hb4w5GngQuNla+z1jzFLg\nHpx/dN4GzrPWDrgdh4iIOFyd8RtjKoBbgcfSnr4W+K619qPA68AFbsYgIiJjub3U0wuswZnZjzgJ\n+Fnq8c+AU12OQVwWj7ezbt0DnHba46xbdz+JRPvoa6+/3swxx3yXaPQBjjnmVnbtap52ux3xOA+v\nW8sTp53Ew+u+REci7kb40zJZH3Pd9o4dL0x7zKYb1/jtmt/YlTdj64d8Orf84OpSj7V2GOgzxqQ/\nvSBtaWcPsMTNGMR9k1WBPOusnxKLXQGE6OlJcuaZG3nrraOn1W4+VWF0q9LlRG0/+ujV9PZeS/qY\n7dz5tVnFNX672PbL+F0sP8bWD/l0bvnB71o9GW8pHhEOV1BcPM+LWDwViVT5HULOxGJh0is1xmLh\n0f61ty8d85rz+/T63xBrGVOFsSHW4tu4TdbHmUh/7/i2+/pWMH7MMu1runGN3669fZmvY+v3+Z8P\n55afY+BH4u8yxpRaa/uAQ4DYZBsnEt3eROWhQqtV0tgYx6mT6NS/b2xMjPavtraF7u4Dr9XWvgkw\nrf63NS4lyfbUO6GtcZlv4zZZH7M1/viPb7u09A16e8eOWaZ9TTeu8dvV1raQ7MaXsc2H89/vc8uj\nWj0ZX/Mj8T8GnAX8MPXzER9ikByarArkAw+cwZlnbiSRWEo4/CYPPPDpabebT1UY3ap0OVHbl176\nCdaund6YTTeu8dtddeVX2XLDO3kxtn7Ip3PLD65W5zTGHAfcBESBAeAt4Bzg34BSoBk431o7lKkN\nVecsPOq/+h/k/oP/1Tnd/nB3BzDRFOQ0N/crIiKZ6c5dEZGAUeIXEQkYJX4RkYBR4hcRCRi/b+AS\nEZFZSCaTdO7v55149+ifd+M9XHfxiRnfo8QvIjIH9PUPOUk90c07e9OSfKKbnr6MV8RPSIlfRCRP\nDA0Ps7ejN5XUe1KzdyfBJ7r6Dtq+eF6IReEKFkcrWFRXweKRP/UVk+5HiV+mJR5vZ/3lj/Li0500\n8jxfWLWXNZtuoyZc53doWemIx3liw/rUHZtRVjfdMus+ZGozHm9nw4ZtvPF6BWXxbVxa9z8Mv+cw\nzrzrDmD+6OvO3bQdNDWdQjhcO2V76dsnk04Bttdfn0c83kx9/ZEsPSTBmtDPqGpp4XvxD9NbfzKH\nHbZ/tP2JZIpFsjOd8yuZTNLe1ccrLe1jEvs78W72JHoYGj74ntW66lKOWh4em9zrKqivLqOoaMqS\nZwdR4pdp2bBhG7945EIgxG6SHLL1M1SWrJ9zFQ3dqMqYqc30iphwFr+IfYZ7//AA9148n1NvuyNj\nZc3ptDeyPZB67l7gCmKxEC+8kKSShwlxKNu5DVLPTVZR1M3qo0GSfuz2v/gS/1pSzxGX/j3vxnt4\nNzGS5Hvo6Rs86L3lpcUcuqgqldTLWVy/gEXhchbVVVA6P7eFKpX4ZVqam6tJr+64m5XUNP/Cz5Bm\npKZ595iqjDXNu11rc/yY7WIlIX5K5a5dE77u/D799ka2d56rZPzxOfDa+O0PlrltmUz/wBB72ntG\nE/sL5Udzxec/yFvhRtoXhAHY9tMXR7cvnhdiYbiCDxzRQHhByZgZfFXFfEKh7GfvM6HEL9MSjXak\nZoJOPcPlPE9HdLnPUWWvIxoluXPHaFXGXPQhU5vjx2wFz5ME9q1YMeHr0WhnVu052ydTz3WRXn1z\nOc8TArZzcPsTyRSLwODQMK1pyf3dRA/vxrvZk+gm3tnHmIWZxuMoGh5iYWcrx+7eQWdtOSd87tMs\nrq9gcbiCutTSjN/1ipT4ZVqamk5hoP9OXnzKWeNfvSo5JysaulGVMVObIxUx33ijgrK92/hEXTNb\n3vNZzty8mYGhzJU1p2rv4Eqc96TW+DdSX38ky5a2s5ok1S1/ZFf8stQaf/ekFUXdrD46FwwOOR+q\nOssxYxP83s5eJqplGa4qxRxay6K6ChaFK1hUV86CeQO8tPEq6prfcI7dhpvz8nMwV6tz5oKqcxYe\n9V/996P/B5J7D3tSiX1Pwknybe29DE+QC6sXlDjr7KnEvihcwcLU76UlM193L+jqnCIiXhpZltmT\nOPDn3YRztUxbx8TJvapiPoc1VqcSejkLw86a+8JwOeWlhZkiC7NXIlKwBgaHaG3vTSX2bt5t72FP\n3JnBZ1qWqa6Yz2GHVLOottxJ8KnEvrC2goqy4KXB4PVYRPJeb//g6Iy9tb1ndHlmT3sPifEfqKbU\nLCjh8ENqRpdjRpZkCnnmPlMaDRHxXDKZZF/PAHvae2gdWZZpT/1J9NC5v3/C94WrSjlyWe1BiX1h\nuJyyEqWz6dJIiYgrhpNJEp19TnJvP5DcE/v6iLXum7C+TFEoRH1NKX+yoo6FqWWZhbXlRFI/S3J8\nI1NQKfGLyIz1DwzR2tHrzNpTCX4kybd19DA4dPCiTElxEZHaciLLnMQeqXU+VI2Ey6mvLqN4nqrF\nu83zxG+MWQD8AAgDJcC11tpfeh2HiEwtmUzS1T3gJPORxJ5ad9/T3kP7vomXZBaUFbNsYaWT4GvL\nD8zewxUcvryevXv3edwTSefHjH8t8LK19h+MMUuAXwHv8yEOEQEGBofZ29k7OlsfmbG3tvfS2tFD\nX//BSzKhENRVlfG+aJhIbZmT3MMVRGrLWFhbTkXZ/Iz7m0lRMcktPxJ/G/D+1OM6oNWHGCRLXlVv\ndKN6ZiEYGZfQ629Mu+LmiGQySWf3AG1pib21PZXoOzJfJVM6f15aUi8fM3uvr5l6SWb8OXPVFcfy\nwsZv0RBroa1x6ZhjO5Pzazbnykzfm837/Kp4OhJj8qH7n7kgmfzQRNtMmfiNMR+31j6Sq6CstT8y\nxqw1xrwK1AJ/nqu2xT1eVW90o3pmIRgZly/w6QkrbvYNDKUS+4GE3paasbe299A/MHxQmyGgNnWV\njJPUy2hIJfZIbfmsi4aNP2di2y/jd7HUsWX7mGM7k/NrNufKTN+bzfv8qniaFuPxmbaZzoz/68aY\n24B/B+6y1jbPJihjzDlAs7V2jTFmJXAnkwQYDldQXFx4n+RHIlV+h5CVWCxMevXGWCw8qz5kem9D\nrGVMZcqGWMucG6vpyLZPdbE3ebd6IW/WnMiymj9SUd1NRe1+eiI1rL/tf2jfd/CXdABUlBVzSKTS\nKfFbVzH6c0nDAhaGy5nv4t+t8edMe/uyjMd2JufXbM6Vmb43m/dN1Se3zuv0GDOZMvFbaz9hjAkD\nZwKbjTEAdwP3W2uz+74vx4nAo6m2nzfGNBpjQtbaCWvyJBLdM9hFfpuLtVoaG+OkV39sbEzMuA+T\n9b+tcSlJto9WpmxrXDbnxmoqE/U/mUzSsb+ftvZe2jp6aO3oHV2aaevoZe9HvsEdq4uoA+rYmf5G\nSoqLOGp5mEhtOQ01ZaNLMpHachaUFWeYtSdpd/nv1vhzpra2hWQ3Ex7bmZxfszlXZvrebN43WZ/c\nzAHpMWYyrTV+a23CGHMv0A9cAvwtcI0x5svW2qezjOs14E+BB4wxUaArU9KX/OFV9UY3qmfmg2Qy\nyf7eQdo6engl1sXrLXHaOnpHE31bRy8DgwcvxwDUVJYQXVzF0GsvsGBPC79rq6dn3ntpjHTznes/\nSn192OPeTM/4c+aqK7/KlhveSa3xLxtzbGdyfs3mXJnpe7N5n18VT0diTD50//YLMmwzZXVOY8yf\nAecDJwP3A7dba18yxiwHHrDWHptNUKnLOe8CFgHzgKustb/OtL2qcxaeQu1/dyqxt3X0ppJ62uOO\nHnonuDoGnEsfG2rKaagtI5L62VDjrLnXV5cV3E1LhXr8szEXqnPeAHwfuMhaO7qQaK3dbYz5cbbB\nWGv3A3+R7ftE/DaS2PeOJvPeMb93T/B1egClJfOI1DjJvKGmjOVLaymfF6I+tSyjOjLitems8X9k\nktc25jYcEX8kk0m6+wbHJPW900zsJfOLaKgp5/ClNU4yTyX4+poyGmrKqCwfe3WMZrziN001JBBG\nioIdSOjOz72dqeTe2Tth7RhwEntkJLFXl40uxTRkSOwi+U6JXwrCcDJJx75+J6l39qSSet+BWXtn\n74TXsoOzFNNQ46ynR2rKR2fqmWbsInOdEr/MCQODwyS6UrP1zgOz9ZGf8c4+hoYnvg5gQVkxi8MV\n1KeS+4Gk7iT5zJc8ihQmJX7xXfr6eryzz0no6cm9s5fOff0TlhUA5ws4oourqK8uG03u6UleH56K\njKW/EeK6oeFh2rv6R5N439Db/PHtTuJpyT3TpY7zikKEq0oxh9ZOmNjrqktdvftUpBAp8cusJJNJ\nevoGnfX0zl4Snc7aerxzZAmml0RX/4Rfcg1QXnpgfb2upoyG6jLqqg8k95oFJarmKJJjSvwyqZG1\n9XhnH/EuJ6mPT+6ZZuuhkPNVeYc1Vh+YqVeXsmJZmGKS1FWVuf5F115XSOyIx/nlX1/K/qf+h3pg\nYNWJfGzTv2RdYXSqdtyuYppt+6qqOrco8QfYyJUw8ZHEnkrkibQkn+m7TwEqSp07TuurS6mrKaOu\nqjS1/OIk+dqqEuYVHVy618vr2L2ukPjEhvXUbv05F5Oq0LL152wpKc26wuhU7bhdxTTb9lVVdW5R\n4i9QyWSSrp6B0SQ+MmNPjCb4Ptr3Zb4SpnheEXVVpTQeWktdKpnXVZdSV+X8rK+eGx+aNjdXk14h\n0fndPTXNu5k/Zo/Oc7lup6Z596z3MdX+s2nf7Xgkt/L/b64cZOQqmPFJfWTWHu/qI9HVl7HoVwin\n8NfyxVXjEvqBpD7bWuz5IhrtSM30nQqJ0Winq/vriEYp2bkjrSYjdESX57ydjmiU5M4ds9rHVPvP\npn2345HcUuLPMyNJPd7Zl7a2fuBxIpXU+wYyV8SurphPY8MC6qpKD5qph6tKqa0sDcwXWntdIXF1\n0y082t/Pjam1+cFVH+GUGVQYnaodt6uYZtt+oVZVLVRTVuf0WyFV5xwpzUvxPF5vjpPoSiX11Cx9\nJMFnusMUoLJ8/mgiD1c5iXzkssZwdRnhylLmF+d3Ug96rRr1P9j9h7lRnVOmYXBomPZ9fbR39ZPY\n58zK27ucJZj2rr7Uc/0MDmVO6lUV81lcV+Ek9epSZ8Y+kuBTv+uadRGZLSX+KTjXqQ+R2Ock8sRo\nEh/7e9f+zHeWhoDqyhKWRhYQriplSaSS8vlFY5ZfwkrqIuKRQCf+4WHn6+7a9x1YOx/5M/rcvj76\nMlynDlBSXERtVSmN9bXUVpUSriwd/Rmudn7WVI69rFH/1RURPxVs4u8bGDpohj46S0/93rEv8x2l\n4KynL6wtH52Rjyb1tD8VpSrwJSJzy5xL/CPXp49J6p19B5ZiUj/39078pRng1H+prSzlsEOqnZl5\n6kqX9IReOwc+JBURmYm8T/z3/fp1Wtt7xizBDA5lnqWXlxZTV1XK8iXVB8/QU8m9smI+RZqli0hA\n+ZL4jTHnAH8HDABXW2u3Ztr25081A07dl9rKUpYtrEpbdimhrqrsQHKvLKW0RB+QiohMxvPEb4yp\nA64GjgWqgG8BGRP/t9edQFlJMdUL5k9Y90VERLLjx4z/VOC/rLXdQDdw0WQbL6lf4ElQc81Mq07m\nSxXFmcQxVZ/TX1+yuJXVAw/y8LMRYqzkfauquOaaVWzcuMOzSp2Z+hmJVGXVRjzezvrLH+XFpztp\n5Hm+sGovazbdNma8Jhub8eOyJvQzFrc080p8L4fW1TH4nsNzfh7kqiqq1+erV9Vc2/fu5eF1X8lp\nv8aP1VUPPRB+K5lMTLStH4l/ObDAGPMQUAt8y1r7Kx/imNNmWnUyX6ooziSOqfo85nWSvMBzvMm/\nACF2b03yv/+7kVjsiozvd8NE/Tz8wfuyamPDhm384pELgRC7SXLI1s9QWbJ+zHhNNjbjx6WSh/kx\nz5ME7o29xbl/eCHn50GuqqJ6fb56Vc116yWX5Lxf48fqFdgMfGGibf1I/CGgDvgMsALYBkQzbRwO\nV1BcgDc2ZTvrGy8WC5NeuzEWC0+rzYZYy5gqig2xllnHMhMziWOqPo9/Pc77x/ze3r500ve7YaJ+\nQnbHf3y/drOShtijk/Y9vW8TvR9+SgioxJ3zYKpjNd19eX2+zvTvVbYqd+3Keb/Gj9URcFimbf1I\n/O8CT1prk8AbxpguY0yDtbZtoo0TiW5vo/NALm7gamyMQ1rtxsbGxLTabGtcSpLto1UU2xqXeX4z\nWSRSNaM4purz+NfreIHutN9ra1vo7s5+zGZjon4CWe13fL+W8/xB4zXZ2Ez0flLP7MOd82CyeLI5\n/70+X2f69ypbXStWkNye236NH6tX4Y1M2/qR+H8J3G2MacKZ+S/IlPQls5lWncyXKooziWOqPqe/\n3rikjY/07+bhZ88nxkqOWlXNNd88gxtu8K5SJ+RmvJuaTmGg/05efMpZ41+9KnlQO5ONzfhxWU2S\nn7y5klf2Omv8W95zRM7Pg1xVRfX6fPWqmusnNm9mS99gTvs1fqx+8ND9F2/MsK0v1TmNMeuAL+P8\nw3SdtfbnmbYtpOqcI4JeskH9V/+D3H8IaHVOa+3twO1+7FtEJOh0YbyISMAo8YuIBIwSv4hIwCjx\ni4gEjBK/iEjAKPGLiASMEr+ISMDk/RexnHba455VUix06ZUHFy9uJRQa5O23l8x6fNOrAr615BAe\n4dPE3m7I2+PmRsXHbNrMtjJjvlRUlcKR94l/587PeFZJsdClVx50bpr+D2D245teFfDzO5fyC5xK\nkvl63Nyo+JhNm9lWZsyXiqpSOObIUk+I5uZqv4OY85wxTK/fVzX6eDbjW9O8e7RVp/LjgX3k43FL\njzeU+t3LNsdXZpxq/27EK8E2RxJ/kmi00+8g5rxotANnpk/qZ9fo49mMb0c0OtqqU/nxwD7y8bil\nx5sEOqLLPW2za8WKrPbvRrwSbHm/1HPMMQ96Vkmx0KVXHlyypA0Y4O23Zz++6VUBVy9Jsp87U2v8\n+Xnc3KglG2H6AAAIA0lEQVT4mE2b2VZmzJeKqlI4fKnOmQ1V5yw86r/6H+T+g//VOefIUo+IiOSK\nEr+ISMAo8YuIBIwSv4hIwCjxi4gEjBK/iEjA+Jb4jTFlxpjXjDFf9CsGEZEg8nPG/4/AXh/3LyIS\nSL7cuWuMMcB7gZ/7sf9CkV5tM18rYc4FE1W/HEoWTWtsg3IMVCG0sPhVsuEm4FJgrU/7Lwjp1Tbz\ntRLmXDBR9cuHOGNaYxuUY6AKoYXF88RvjDkPeNJa2+xM/Ml4WzFAOFxBcfE8T2LzUiRSNfVGU4jF\nwqRXwozFwjlp1wv5FGdDrGVM9cuGWAsxpje2Mz0G+dT/6ZhojGbTh7nWfzf4OQZ+zPj/HFhhjPkU\nsBToNca0WGt/NdHGiUS3p8F5IVd1Ohob4zj1Gp36+o2NiTlRAyXfarW0NS4lyfbRbyloa1xGI9Mb\n25kcg3zr/3RMNEYz7cNc7H+ueVSrJ+Nrnid+a+0XRh4bY64BdmVK+jK59Gqb+VoJcy6YqPrliRQx\nnbENyjFQhdDC4mt1zrTE/4NM26g6Z+FR/9X/IPcf/K/O6Ws9fmvtt/zcv4hIEOnOXRGRgFHiFxEJ\nGCV+EZGAUeIXEQkYJX4RkYBR4hcRCRglfhGRgPH1Ov4gisfbueyyh3nllfKCruaYD4JSOTMfqZpn\nflPi91hQqjnmA421f1TNM79pqcdjzc3VpFdzdH4XN2is/VPTvHtMNc+a5t0+RiPjKfF7LBrtwKlv\nCJAkGu30M5yCprH2T0c0mjby0BFd7mM0Mp6WejzW1HQKpaX3ptb4C7eaYz4ISuXMfKRqnvnN1+qc\n06HqnIVH/Vf/g9x/8L86p5Z6REQCRolfRCRglPhFRAJGiV9EJGCU+EVEAsaXyzmNMU3AR4B5wD9Z\nax/wIw4RkSDyfMZvjDkJOMpa+2FgDbDJ6xhERILMj6WeXwNnpx63AxXGmIzXm4qISG55vtRjrU0C\nPalfvwz8IvWcSM6oMqdIZr6VbDDGnAGcD5zmVwxSuFSZUyQzvz7cPR24AjjdWjvpfcvhcAXFxfO8\nCcxDkUiV3yH4yu3+x2Jh0itzxmLhvBrzfIrFD0HvP/g7Bp4nfmNMNdAEfMxa2zHV9olEt/tBeSzo\ntUq86H9jYxynLmQISNLYmMibMdfxD3b/wbNaPRlf82PG/xdAPfDj1Ie6SeCL1to3fYhFCpQqc4pk\n5seHu7cDt3u9XwmWcLhWa/oiGejOXRGRgFHiFxEJGCV+EZGAUeIXEQkYJX4RkYBR4hcRCRglfhGR\ngFHiFxEJGCV+EZGAUeIXEQkYJX4RkYBR4hcRCRglfhGRgFHiFxEJGCV+EZGAUeIXEQkYJX4RkYBR\n4hcRCRglfhGRgPHjy9YxxtwM/CkwDFxurf29H3GIiASR5zN+Y8yfAYdbaz8MfBm41esYRESCzI+l\nno8BDwJYa18Gao0xlT7EISISSH4k/sVAa9rvbannRETEA/nw4W7I7wBERILEjw93Y4yd4TcCb2fa\nOBKpKsh/GCKRKr9D8JX6r/4HnZ9j4MeM/5fA5wCMMccBb1lr9/sQh4hIIIWSyaTnOzXG3AB8FBgC\nLrXWvuB5ECIiAeVL4hcREf/kw4e7IiLiISV+EZGAUeIXEQkYX2r1BJExZgHwAyAMlADXWmt/6W9U\n3jHGhIDvA0cDfcBF1tpX/I3KG8aYo3HuVr/ZWvs9Y8xS4B6cidfbwHnW2gE/Y3TT+P6nnvs6cCNQ\na63t9jM+L0xwDiwD7gLmA/3AudbaPV7Foxm/d9YCL1trTwHOBv6Pv+F47gyg2lp7Ik6Nppt8jscT\nxpgKnHpUj6U9fS3wXWvtR4HXgQv8iM0LE/XfGHMesBB4y6+4vJThHLgO+L619iScfxD+xsuYlPi9\n0wbUpx7XMbZsRRAcATwDYK19A4im/hdQ6HqBNYy9SfEk4Gepxz8DTvU4Ji9N1P/7rbVX+RSPHyYa\ng4uB+1OPW3FygmeU+D1irf0RTrJ7Ffhv4G/9jchzLwCnG2OKjDEGWAE0+ByT66y1w9bavnFPL0hb\n2tkDLPE4LM9M1P+g3bCZYQx6rLVJY0wRcCnwQy9jUuL3iDHmHKDZWnsEToXSf/E5JE9Zax/BmfH/\nGvg68BKq0wQag8BKJf17gMettdu83LcSv3dOBB4FsNY+DzQGZKljlLX2amvtamvtpUCdlx9m5Zku\nY0xp6vEhOPWrgijod4/eDVhr7XVe71iJ3zuv4XzrGMaYKNBlrQ3MiW+MWWmMuTP1+OPAsz6H5KfH\ngLNSj88CHvExFi+Nn+gEauKTLrUC0GetvdaP/atkg0dSl3PeBSwC5gFXWWt/7W9U3kn97+ZO4E+A\nHuAca23BX9WRKkR4ExAFBnCuZDkH+DegFGgGzrfWDvkWpIsy9P+/gNOAE4DtwFPW2m/4FqTLMozB\nQpwPfbtw/ufzorX2Mq9iUuIXEQkYLfWIiASMEr+ISMAo8YuIBIwSv4hIwCjxi4gEjBK/iEjAKPGL\niASMEr+ISMAo8YtkyRjz18aY/5t6bIwxL6XuzBaZE5T4RbK3CTjSGPNhnCqr64JWaljmNiV+kSyl\niutdCPwYeN5a+1ufQxLJihK/yMzU4xTYOtTvQESypcQvkiVjTBmwGfgU0G+MOdfnkESyosQvkr1v\n4Xxv7GvA5cA3jTGNPsckMm0qyywiEjCa8YuIBIwSv4hIwCjxi4gEjBK/iEjAKPGLiASMEr+ISMAo\n8YuIBIwSv4hIwPx//qS8YITBGdIAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ax = d.sort_values(by=['x']).plot(x='x', y='yhat')\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

施肥処理Fをモデルに組み込む

\n", "\t

\n", "\t\tつぎに施肥処理fを使ってポアソン回帰をしてみましょう。xと異なりfには数字では無くT, Cのような文字列（因子型）が入っています。\n", "\t

\n", "\t

\n", "\t\t因子型で線形予測子を作る場合、因子の数より1個少ない数の$d_i$というダミー変数を使って、因子の影響を線形予測子として表します。\n", "\t\t$d_i$が、$f_i=T$の場合1， $f_i!=T$の場合0とします。因子型ではいずれか一つの因子に当てはまる（昔のカーラジオのボタンのように）\n", "\t\tとしているので、$f_i$がTでなければCとなるので、因子がTとCの2個の場合ダミー変数は、$d_i$1個で足ります。\n", "\t

\n", "

\n", "\t\tもし、肥料が複数の場合には、ダミー変数は、$d_i,A, d_i,B, d_i,C$のように複数作らなくてはなりません。\n", "\t

\n", "\t

\n", "\t\tglmは、とてもお利口なので人がダミー変数を指定する代わりに因子型の変数fを指定するだけで内部的にダミー変数を考慮した計算をしてくれます。\n", "\t

\n", "\t

\n", "\t\tモデル式をy ~ xからy ~ fに変えてポアソン回帰を実行してみます。\n", "\t

\n", "\t

\n", "\t\t回帰の結果は、summary2関数で表示します。結果は、\n", "\t\t

対数尤度(Log-Likelihood): -237.63
逸脱度(deviance): 89.475
切片$\\beta_1$(Intercept): 2.0516、標準誤差は0.0507
fの係数$\\beta_3$: 0.0128、標準誤差は0.0715

\n", "\t\t先の例と比べると、対数尤度は-237.63と小さく求まっています。\n", "\t

\n", "\t

\n", "\t\t久保本は、とても丁寧でこの結果をどう考えるかを解説してくれています。\n", "$$\n", "\t\t\\lambda_i = exp(2.05 + 0.0128) = exp(2.05) exp(0.0128)\n", "$$\t\n", "\t\t施肥を行った場合には、施肥を行わない場合（0.0128の項が０）の1.0128倍に若干増えると予測しています。\n", "$$\n", "\t\t\\lambda_i = exp(2.05 + 0.0128) = 1.0128 exp(2.05) \n", "$$\t\n", "\t

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

Model:	GLM	AIC:	479.2545
Link Function:	log	BIC:	-361.8317
Dependent Variable:	y	Log-Likelihood:	-237.63
Date:	2016-07-17 01:38	LL-Null:	-237.64
No. Observations:	100	Deviance:	89.475
Df Model:	1	Pearson chi2:	87.1
Df Residuals:	98	Scale:	1.0000
Method:	IRLS

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

	Coef.	Std.Err.	z	P>\|z\|	[0.025	0.975]
Intercept	2.0516	0.0507	40.4630	0.0000	1.9522	2.1509
f[T.T]	0.0128	0.0715	0.1787	0.8582	-0.1273	0.1529

" ], "text/plain": [ "\n", "\"\"\"\n", " Results: Generalized linear model\n", "==============================================================\n", "Model: GLM AIC: 479.2545 \n", "Link Function: log BIC: -361.8317\n", "Dependent Variable: y Log-Likelihood: -237.63 \n", "Date: 2016-07-17 01:38 LL-Null: -237.64 \n", "No. Observations: 100 Deviance: 89.475 \n", "Df Model: 1 Pearson chi2: 87.1 \n", "Df Residuals: 98 Scale: 1.0000 \n", "Method: IRLS \n", "---------------------------------------------------------------\n", " Coef. Std.Err. z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------\n", "Intercept 2.0516 0.0507 40.4630 0.0000 1.9522 2.1509\n", "f[T.T] 0.0128 0.0715 0.1787 0.8582 -0.1273 0.1529\n", "==============================================================\n", "\n", "\"\"\"" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 施肥処理Fをモデルに組み込む\n", "fit_f= smf.glm('y ~ f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "fit_f.summary2() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

個体のサイズと施肥を合わせたモデル

\n", "\t

\n", "\t\t最後に個体のサイズ$x_i$と施肥$f_i$を合わせたモデルを使ってポアソン回帰を行ってみます。\n", "\t

\n", "\t

\n", "\t\tモデル式は、y ~ x + fとなります。\n", "\t

\n", "\t

\n", "\t\t結果は以下の様になります。\n", "$$\n", "\t\t\\lambda_i = 1.26 + 0.08 x_i - 0.032\n", "$$\t\t\n", "\t

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

Model:	GLM	AIC:	476.5874
Link Function:	log	BIC:	-361.8936
Dependent Variable:	y	Log-Likelihood:	-235.29
Date:	2016-07-17 01:38	LL-Null:	-237.64
No. Observations:	100	Deviance:	84.808
Df Model:	2	Pearson chi2:	83.8
Df Residuals:	97	Scale:	1.0000
Method:	IRLS

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

	Coef.	Std.Err.	z	P>\|z\|	[0.025	0.975]
Intercept	1.2631	0.3696	3.4172	0.0006	0.5386	1.9876
f[T.T]	-0.0320	0.0744	-0.4302	0.6670	-0.1778	0.1138
x	0.0801	0.0370	2.1620	0.0306	0.0075	0.1527

" ], "text/plain": [ "\n", "\"\"\"\n", " Results: Generalized linear model\n", "==============================================================\n", "Model: GLM AIC: 476.5874 \n", "Link Function: log BIC: -361.8936\n", "Dependent Variable: y Log-Likelihood: -235.29 \n", "Date: 2016-07-17 01:38 LL-Null: -237.64 \n", "No. Observations: 100 Deviance: 84.808 \n", "Df Model: 2 Pearson chi2: 83.8 \n", "Df Residuals: 97 Scale: 1.0000 \n", "Method: IRLS \n", "---------------------------------------------------------------\n", " Coef. Std.Err. z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------\n", "Intercept 1.2631 0.3696 3.4172 0.0006 0.5386 1.9876\n", "f[T.T] -0.0320 0.0744 -0.4302 0.6670 -0.1778 0.1138\n", "x 0.0801 0.0370 2.1620 0.0306 0.0075 0.1527\n", "==============================================================\n", "\n", "\"\"\"" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 説明変数が数量型＋因子型のモデル\n", "fit_all = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "fit_all.summary2()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

Nullモデルとの比較

\n", "\t

\n", "\t\tついでなので、λが定数の場合（Nullモデル）のポアソン回帰の結果を以下に示します。\n", "\t

\n", "\t

\n", "\t\tここで注目するのは、Log-Likelihoodの値-237.64です。\n", "\t\tλが定数の場合の対数尤度と自分の考えたモデルの対数尤度を比べることで、\n", "\t\t自分のモデルがどの程度よいのか比較できることです。\n", "\t

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

Model:	GLM	AIC:	477.2864
Link Function:	log	BIC:	-366.4049
Dependent Variable:	y	Log-Likelihood:	-237.64
Date:	2016-07-17 01:38	LL-Null:	-237.64
No. Observations:	100	Deviance:	89.507
Df Model:	0	Pearson chi2:	87.1
Df Residuals:	99	Scale:	1.0000
Method:	IRLS

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

	Coef.	Std.Err.	z	P>\|z\|	[0.025	0.975]
Intercept	2.0580	0.0357	57.5862	0.0000	1.9879	2.1280

" ], "text/plain": [ "\n", "\"\"\"\n", " Results: Generalized linear model\n", "==============================================================\n", "Model: GLM AIC: 477.2864 \n", "Link Function: log BIC: -366.4049\n", "Dependent Variable: y Log-Likelihood: -237.64 \n", "Date: 2016-07-17 01:38 LL-Null: -237.64 \n", "No. Observations: 100 Deviance: 89.507 \n", "Df Model: 0 Pearson chi2: 87.1 \n", "Df Residuals: 99 Scale: 1.0000 \n", "Method: IRLS \n", "---------------------------------------------------------------\n", " Coef. Std.Err. z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------\n", "Intercept 2.0580 0.0357 57.5862 0.0000 1.9879 2.1280\n", "==============================================================\n", "\n", "\"\"\"" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Nullモデルの逸脱度とAIC\n", "fit_null = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()\n", "fit_null.summary2()" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "89.5069375696 -237.643221309 477.286442619\n" ] } ], "source": [ "# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットされている\n", "print fit_null.deviance, fit_null.llf, fit_null.aic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

3章注目のグラフ

\n", "\t

\n", "\t\t「Sageでグラフを再現してみよう」ではいろんなグラフをSageで再現してその本質を理解できるように挑戦しています。\n", "\t

\t\n", "\t

\n", "\t\t3章の注目のグラフは、図3.9（以下久保本から引用します）の（Ｂ）です。\n", "\t\t種子数λがサイズxに対して指数関数で増加するとそのポアソン分布は裾野が広く平らに分布しています。\n", "\t\t今回はこのグラフをSageを使って計算してみます。\n", "\t

\n", "\t

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

\n", "\t

\n", "\t\tGLMの結果が示されていないので、$x$が0.5, 1.1, 1.7のyの値をグラフから読み取り、\n", "\t\t0.1, 0.5, 3.0としてポアソン分布をプロットしたのが、以下の図です。図３．９（Ｂ）のポアソン分布\n", "\t\tとよく似た分布になっています。\n", "\t

\n", "\t\n", "" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD/CAYAAADsfV27AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFNRJREFUeJzt3X+QXeV93/H3IrnAajf6gVYWIgEhIX9xTO0WZ2xFNhQs\nuY4bSuoWpkzcFqdVbYcfFoiEJjRgjH+QQZUAxWaIoU7qTmymcaYYG6hxnNq1gcYWLhm1E30TJFYy\nkiyttCtF2jUCSds/9kqsVrt7r7Rnd+8jvV8zDPee5+w9H0mjzz167j3Paenv70eSVK4zJjuAJGls\nLHJJKpxFLkmFs8glqXAWuSQVziKXpMJNbWSniLgEeBxYk5kPDRlbBnwWOAg8nZmfqTylJGlEdc/I\nI6IVWAv8+Qi7PAh8CHgv8I8j4uLq4kmS6mlkauVV4IPA9qEDEXEhsDszt2VmP/AUsLTaiJKk0dQt\n8sw8nJkHRhieC3QNer4TOLeKYJKkxlT9YWdLxa8nSaqjoQ87R7GNY8/Az6ttG9HBg4f6p06dMsbD\nStJpZ8QT5RMt8mNeKDM3R0R7RJzPQIFfBfz6aC/Q09N3goesr6Ojna6ufZW/btXMWS1zVqeEjHB6\n5+zoaB9xrG6RR8SlwGrgAuD1iPgXwBPAy5n5deA3gceAfuCrmflSFaElSY2pW+SZ+WPgylHGfwAs\nqTKUJKlxXtkpSYWzyCWpcBa5JBXOIpekwlnkklS4sV4QJElN6dChQ3R2bqr0NefPX8CUKaNf0PgH\nf7CGv/mbv+bgwcOsWHEbF1/8i0fHXnvtNVat+hwvv7yJRx/9cmW5LHJJp6TOzk2sWPUErdPnVPJ6\nfXt38uBvX83ChYtG3OfFF3/MK6/8hMcee4x169Zz77338PDDXzo6/tBDD7JoUVT+BmORSzpltU6f\nQ9vM8ybseC+88CMuu+wKAC64YD779++jr6+P1tZWAD72sZvYu3cP3/7205Ue1zlySarI7t27mDFj\n5tHn06fPoLt799HnZ5999rgc1yKXpHHS398/IcexyCWpIrNndxxzBr5rVxfnnDN73I9rkUtSRd71\nrsV897vfASBzAx0dc46bTunv76fqE3U/7JR0yurbu3NCX+uSS95OxFu57rrrOHwYVq78Dzz99Ddp\na2vjssuu4M47f4edO3fwk59s5hOf+DhXX/0hli37wJiztUzUHM4RXV37Kj/g6bxG8XgwZ7VKyFlC\nRjixnJP1PXIYt/XIK7uxhCQVYcqUKaN+5/tU4hy5JBWuac/I6/2zqNF/4kjSqa5pi3y0y2sbuVRW\nkk4XTVvkMPGX10pSiZwjl6TCNfUZuSSdrGZcxvbaa6/mzW+eS0tLCy0tLdx112eYPXvsV35a5JJO\nSZ2dm7j9ibuY1tFeyev1du3jvqvvGdMyttDC6tVrOfPMsyrJdIRFLumUNa2jnfZ5MybsePWWsYXq\nL88H58glqTL1lrEFWLXqc9xww3L+8A+/UNlxLXJJGidDl0BZvvzj3HzzSj7/+S+yceNLfO97f1HJ\ncSxySapIvWVsP/CBf8KMGTM444wz+OVffg8bN75UyXEtckmqyGjL2Pb27mflyps5ePAgAC+++AIL\nFiys5Lh+2CnplNVb4QqEjbxWvWVslyx5Dx/96Ec466yzeMtbgiuuWFpJNotc0ilp/vwF3Hf1PZW/\nZj0f+9iNxyxju3DhRUfHrrnmOq655rpKM4FFLukU5TK2kqRiWOSSVDiLXJIKZ5FLUuEsckkqnN9a\nkXRKmoxlbA8ceJXPfvZT7N+/l97en3H99f+OJUvee3T8Rz/6S774xYeYOnUK7373Ej7ykeWV5Gqo\nyCNiDbAYOAzckpnrBo3dCHwYOAisy8yVlSSTpDHo7NzEc7d+gnOPrjw4Ntv7+uD+taN+pfEHP/g+\nF1/8i6xYcQPr1/8tt956wzFF/uCDq3nggS9wzjmzuemmj3Lllcu44IL5Y85Wt8gj4nLgosxcEhEX\nA18CltTG2oHfAhZkZn9EfCsi3pWZPxxzMkkao3NbWzm/rZr1yBuxdOn7jz7eseOnzJkz9+jzbdu2\nMn36dGbP7gBg8eL3sG7dDyemyIGlwOMAmbkhImZERFtm7gdeAw4APxcRvcDZQPeYU0lSwa677jq2\nb/8p9913/9Ft3d27mTHjjbXRZ86cybZtWys5XiMfds4FugY931XbRmYeAO4BNgEvA3+ZmdUs5yVJ\nhXrssce4997VfOpTd46yV3V3mDiZb620HHlQm1q5A7gIuBBYHBF/v6JsklSUzA3s3LkDgEWL3sKh\nQ4fYs2cPMLDE7e7dbyxx29XVdXSaZawamVrZRu0MvGYesL32+K3AxszsAYiI7wPvBNaP9GIzZ7Yy\nderoNy8F6OlpG3V81qw2Ogbdi6+jovvyjTdzVsuc1SkhIzSes6enjZcrPvbQ3hnqySf/H9u2beNt\nb7uDlpYDvP76ARYt+gVgIPdrr73K66/vY86cOfzwh8+xevXqSn7fGynyZ4C7gUci4lJga2b21sY6\ngbdGxJm1aZZfAp4c7cV6evoaCtbdvb/u+JHVxQavNNbMzFktc1anhIxwYjm7u/cPfNOkItv7+o7p\nneEsW3YVv//7n+bDH/4wvb19rFjx23z5y189uoztihW3c/PNK2hpgSuuWEZr66yGfz2jFX7dIs/M\n5yPihYh4FjgE3BgR1wN7MvPrEbEK+G5EvA48l5nPNpRKksbR/PkL4P61lb3ehdRfxvbMM8/kk5/8\nzIhvOO94xz/g4Ye/VFmmIxr6Hnlm3jFk0/pBY48Aj1QZSpLGymVsJUnFsMglqXAWuSQVziKXpMJZ\n5JJUOItckgpnkUtS4SxySSqcRS5JhbPIJalwFrkkFc4il6TCWeSSVDiLXJIKZ5FLUuEsckkqnEUu\nSYWzyCWpcBa5JBXOIpekwlnkklQ4i1ySCmeRS1LhLHJJKpxFLkmFs8glqXAWuSQVziKXpMJZ5JJU\nOItckgpnkUtS4SxySSqcRS5JhbPIJalwFrkkFW5qIztFxBpgMXAYuCUz1w0a+3ngq8CbgB9n5g3j\nEVSSNLy6Z+QRcTlwUWYuAZYDa4fsshpYlZmLgUO1YpckTZBGplaWAo8DZOYGYEZEtAFERAvwXuAb\ntfGbM/OVccoqSRpGI1Mrc4F1g57vqm17CegA9gMPRMSlwPcz847KUw7Rf/gwW7ZsPvq8p6eN7u79\nR5/Pn7+AKVOmjHcMSWoKDc2RD9Ey5PF5wP3AFuDJiPhgZj5dRbiR/GxfF59/8X8wbWv7cWO9Xfu4\n7+p7WLhw0XhGkKSm0UiRb2PgDPyIecD22uNdQGdmdgJExHeAtwEjFvnMma1MnVr/bLmnp23U8Wkd\n7bTPmzHs2KxZbXR0HF/yzaBZcw1lzmqVkLOEjGDO4TRS5M8AdwOP1KZPtmZmL0BmHoqITRGxMDM3\nAu8EvjLai/X09DUUbPBUyYnq7t5PV9e+k/758dLR0d6UuYYyZ7VKyFlCRji9c472xlC3yDPz+Yh4\nISKeBQ4BN0bE9cCezPw6cCvwx7UPPtdn5jcqyi1JakBDc+TDfIC5ftDYRuCyKkNJkhrnlZ2SVDiL\nXJIKZ5FLUuEsckkqnEUuSYWzyCWpcBa5JBXOIpekwlnkklQ4i1ySCmeRS1LhLHJJKpxFLkmFs8gl\nqXAWuSQVziKXpMJZ5JJUOItckgpnkUtS4SxySSqcRS5JhbPIJalwFrkkFc4il6TCWeSSVDiLXJIK\nZ5FLUuEsckkqnEUuSYWzyCWpcBa5JBXOIpekwlnkklQ4i1ySCmeRS1LhLHJJKtzURnaKiDXAYuAw\ncEtmrhtmn3uBxZl5ZbURJUmjqXtGHhGXAxdl5hJgObB2mH3eClwG9FeeUJI0qkamVpYCjwNk5gZg\nRkS0DdlnNXBHxdkkSQ1opMjnAl2Dnu+qbQMgIq4H/iewudpokqRGNDRHPkTLkQcRMRP4DQbO2n9h\n8JgkaWI0UuTbGHQGDswDttcevw+YDXwfOAtYEBGrM/O2kV5s5sxWpk6dUvegPT1DZ28aN2tWGx0d\n7Sf98+OpWXMNZc5qlZCzhIxgzuE0UuTPAHcDj0TEpcDWzOwFyMw/A/4MICIuAP5otBIH6OnpayhY\nd/f+hvYb6We7uvad9M+Pl46O9qbMNZQ5q1VCzhIywumdc7Q3hrpz5Jn5PPBCRDwLPADcGBHXR8Sv\nVRdRknSyGpojz8yh30hZP8w+mxmYapEkTSCv7JSkwlnkklQ4i1ySCmeRS1LhLHJJKpxFLkmFs8gl\nqXAWuSQVziKXpMJZ5JJUOItckgpnkUtS4SxySSqcRS5JhbPIJalwFrkkFe5kbr6sBhw6dIjOzk3H\nbOvpaTt6C7v58xcwZUr9e5dKUj0W+Tjp7NzE7U/cxbRh7rPX27WP+66+h4ULF01CMkmnGot8HE3r\naKd93ozJjiHpFOccuSQVziKXpMJZ5JJUOItckgpnkUtS4SxySSqcRS5JhbPIJalwFrkkFc4il6TC\nWeSSVDiLXJIKZ5FLUuEsckkqnEUuSYVzPXKNaLi7HA3WLHc5KiWnNF4sco2os3MTz936Cc5tbT1u\nbHtfH9y/tinuclRKTmm8NFTkEbEGWAwcBm7JzHWDxq4EPgccBDIzl49HUE2Oc1tbOb/t+NvVNZtS\nckrjoe4ceURcDlyUmUuA5cDaIbs8DPzzzLwM+LmI+JXqY0qSRtLIh51LgccBMnMDMCMi2gaNvzMz\nt9cedwHnVBtRkjSaRop8LgMFfcSu2jYAMnM/QEScC7wfeKrKgJKk0Z3M1w9bhm6IiDnAE8BvZmbP\nmFNJkhrWyIed2xh0Bg7MA45MpRAR7Qychf9uZn6n3ovNnNnK1Kn1vwrW09NWd5+RzJrVRkfH5H7w\nVS9/M2QcTUdHOz09bbw8yj7N8GsoKWezKyEjmHM4jRT5M8DdwCMRcSmwNTN7B42vAdZk5rcbOWBP\nT19Dwbq79ze030g/29W176R/vgr18jdDxpF0dLTT1bWv6X8NpeVsZiVkhNM752hvDHWLPDOfj4gX\nIuJZ4BBwY0RcD+xhoOT/FbAwIv490A98JTMfrSS5JKmuhr5Hnpl3DNm0ftDjs6uLI0k6Ua61IkmF\ns8glqXCutSJNEBf30nixyKUJ4uJeGi8WuTSBXNxL48E5ckkqnGfkk6D/8GG2bNk84rhzpZJOhEU+\nCfp27eeVr63mkHOlkipgkU8S50olVcU5ckkqnEUuSYWzyCWpcBa5JBXOIpekwlnkklQ4i1ySCmeR\nS1LhvCBI0jFcbrc8FrmkY7jcbnks8jEa6exltEWxpGbnEhJlscjHqLNzEytWPUHr9DnHbN/9yl9z\n3gcnKZSk04pFXoHW6XNom3neMdv69u4Adk1OIEmnFYv8NDfc1FBPTxvd3fudHpIKYZGf5jo7N3H7\nE3cxreP4+dCuDdu5jb83CakknQiLXEzraKd93ozjtvfu/Dvg9YkPJOmEeEGQJBXOIpekwlnkklS4\nU26O3DvUS6cHlxJ4wylX5N6hXjo9uJTAG065IgcvL5ZOF/5dH+AcuSQVziKXpMJZ5JJUuFNyjlyS\nmsVEfLumoSKPiDXAYuAwcEtmrhs0tgz4LHAQeDozPzOmRJJ0CpmIb9fULfKIuBy4KDOXRMTFwJeA\nJYN2eRB4P7Ad+F5EfC0zN4wplTSEqzSqZOP97ZpGzsiXAo8DZOaGiJgREW2ZuT8iLgR2Z+Y2gIh4\nqra/Ra5KuUqjNLJGinwusG7Q8121bS/V/t81aGwnsKCydKrEaHN0JZ3NukqjNLyT+bCz5STHTljf\n3p3Dbv/Zvm56u/YN/zPdvWzvOzjs2Pa+Pi6sLN2gYw6Ts5kydnZu4qN3PspZbbOOG9u7YxNvvqI5\ncgJs3Pi3w27fsmVz0/x+wsg5hzN4Cmh7X9+w+0x2ziMZAXOOollztvT394+6Q0R8EtiWmY/Unm8E\n3p6ZvRFxAfDVzFxSG7sL2JWZD1WQTZLUgEa+R/4McA1ARFwKbM3MXoDM3Ay0R8T5ETEVuKq2vyRp\ngtQ9IweIiM8B/wg4BNwIXArsycyvR8R7gfuAfuBrmXn/OOaVJA3RUJFLkpqXl+hLUuEsckkqnEUu\nSYWzyCWpcMUWeUS0RcRFtf+mTXaeExERx1+eOIki4rgLuSLi5ycjS6MiYvZkZ2hERLxvsjPUExFT\nI+KC2leIm1opf+4TrbhvrUTELwFrgRkMLBfQAswDtgI3Zub6SYzXkIj4i8yc9L/gEfEh4AGgFXgK\nuCkz99XGmiIjQET8KrAG+AlwC/AnDFyVPA24ITOfmsR4R0XEvxmyqQX4PeDTAJn55QkPNYyIeDAz\nV9QeLwP+M/BTYA7w8cz81mTmOyIiPgj8WmZ+vPaG+EfAPgb+3G/KzCcnNWBNRPwd8F+AT2fm8Jej\nj7OmfwcexgPAvx26wmLtYqUvAJdPSqohIuKGEYZagPMmMssofgf4h8AeYDnwTET8SmbupeLlFsbo\n9xhYYfN84JsM/OX+q4h4M/ANBt6EmsFdwG7gSd74/TsLxuVq8bF4+6DHdwFXZuamiJgL/HegKYoc\nuIeBiwwBPskbOc9h4Pe4KYoceAH4U+ArEbEF+GPgucwcft2IcVBikZ8x3DK5mfnjiBjb6uzVWgn8\nOQPL+w71pgnOMpJDmdlde/zFiNgBfCsirmLgAq9mcSAztwBbImJrZv4VQGbuiIhXJznbYJcAdwLv\nAFZm5ubaG+OnJjnXUIP/bLszcxNAZv40Ippp9bE3MXAGDgMnGy/XHnfTXCca/Zn5v4BltRmD5Qz8\nfdoH7MzMXx3vACUW+f+OiCcYWFr3yMqLcxlYRuB7k5bqeP+MgSmgFZl5YPBARFwxKYmO94OI+CZw\nbWb+rHal7qvAd4BzJjnbYDsi4rcy8z9l5nvg6Bz+bQxMtzSFzHwV+I8REcAXIuI5mvNzqEsi4r8x\nUIaLIuLazPzTiLiNgcJsFquA/xMR32agvB+v/Z6+D3h0UpMd6+ibSu2mO+sAIuJc4NyJCFBckWfm\nytrNLpYC765t3gbcnZnPT16yY2Xm/62d2Q53hnPbROcZTmbeXntTeXXQtm9FxPPAv5y0YMf7CPBP\nh2ybA2wGfnfC09SRmQlcFRH/mjfOIpvJtUOeH1nSbzvw6xOcZUSZ+ScR8TSwDJjPQGHuAH7jyD0Q\nmsR/HW5jZm5n+H+RV664DzslScdqxn/2SZJOgEUuSYWzyCWpcBa5JBXOIpekwv1/URQzolVGyt8A\nAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 3章のメインのグラフは、図3.9だと思う。\n", "# 0.5, 1.1, 1.7のyを0.1, 0.5, 3.0としてポアソン分布を表示\n", "# Fig. 2.6\n", "# ポアソン関数p(y, λ)を定義\n", "def p(y, lam):\n", " return float(lam^y*e^(-lam)/factorial(y) )\n", "# データフレームを作成\n", "df = pd.DataFrame(\n", " {\n", " '0.1': [p(y, 0.1) for y in range(8)],\n", " '0.5': [p(y, 0.5) for y in range(8)],\n", " '3.0': [p(y, 3.0) for y in range(8)]\n", " })\n", "# データフレームのプロット\n", "df.plot.bar()\n", "plt.show()" ] }, { "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 }