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

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

\n", "\t

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

\n", "\t

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

\n", "\t

\n", "\t\t久保本もようやくテーマの個体差が大きい場合の例題に入りました。今回のグラフは図7.10です。以下に久保本から引用します。\n", "\t

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

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

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

前準備

\n", "\t

\n", "\t\t最初に必要なライブラリーやパッケージをロードしておきます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 1, "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", "# statsmodelsを使ってglmを計算します\n", "import statsmodels.formula.api as smf\n", "import statsmodels.api as sm\n", "from scipy.stats.stats import pearsonr\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７章の例題は、個体$x_i$で観測された以下のデータを含みます。\n", "\t\t

\n", "\t

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

データの特徴と可視化

\n", "\t

\n", "\t\tいつもの通りデータを読み込み、その内容、統計情報をチェックします。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 2, "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", " \n", " \n", " \n", " \n", " \n", " \n", "
Nyxid
08021
18122
28223
38424
48125
\n", "
" ], "text/plain": [ " N y x id\n", "0 8 0 2 1\n", "1 8 1 2 2\n", "2 8 2 2 3\n", "3 8 4 2 4\n", "4 8 1 2 5" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 7章のデータを読み込む\n", "d = pd.read_csv('data/ch7_data.csv')\n", "d.head()" ] }, { "cell_type": "code", "execution_count": 3, "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", " \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", "
Nyxid
count100.0100.000000100.000000100.000000
mean8.03.8100004.00000050.500000
std0.03.0705341.42133829.011492
min8.00.0000002.0000001.000000
25%8.01.0000003.00000025.750000
50%8.03.0000004.00000050.500000
75%8.07.0000005.00000075.250000
max8.08.0000006.000000100.000000
\n", "
" ], "text/plain": [ " N y x id\n", "count 100.0 100.000000 100.000000 100.000000\n", "mean 8.0 3.810000 4.000000 50.500000\n", "std 0.0 3.070534 1.421338 29.011492\n", "min 8.0 0.000000 2.000000 1.000000\n", "25% 8.0 1.000000 3.000000 25.750000\n", "50% 8.0 3.000000 4.000000 50.500000\n", "75% 8.0 7.000000 5.000000 75.250000\n", "max 8.0 8.000000 6.000000 100.000000" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

同じ点に重なるデータの可視化

\n", "\t

\n", "\t\tプロットデータを見ると、とても１００個のデータがあるようには見えません。\n", "\t\tこれは、種子数・葉数が整数で同じ点に複数のデータが重なっているからです。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAFgCAYAAAAcmXr5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFs1JREFUeJzt3X9sXfd5mPHnUJecTJqGaJVKKTt1ugz7EsImQ1ld1FKc\nZkmxNKoLZ2u7FsoSJbMDGM4wt4ACTEXrKN2wBIlSNMXqFoubjAUqbMmQpAXCdsmgNVtlt9YAYtqw\n8rvGKTKL4URapWb+EEtR9+yPS3pSWpPiFa/eew6fz1889qXv+/ryPufwkLaKsiyRJN1ZPdEDSNJO\nZHwlKYDxlaQAxleSAhhfSQrQiB7gRqur18u5uaXoMbbF0FA/7tJ96rQL1GufOu0yPDxYbPaYrrry\nbTR2RY+wbdylO9VpF6jXPnXa5VZ0VXwlaacwvpIUwPhKUgDjK0kBjK8kBTC+khTA+EpSAOMrSQGM\nryQFML6SFMD4SlIA4ytJAYyvJAUwvpIUwPhKUgDjK0kBjK8kBTC+khTA+EpSAOMrSQGMryQFML6S\nFMD4SlIA4ytJAYyvJAVodPIfnlIaAH4LGAL6gF/KOX+tk88pSVXQ6SvfDwCTOed3AD8FfKbDzydJ\nldDRK1/gFeBvr318LzDb4efTNlptNhkbn2R6bomRoX6OHx2l0VPNO1Ur169z+swEM1eW2bdnNyeO\nHaJv167osdq2dO0aJ599gcXlVQZ2N/j4Uw/T39sbPVZb/vzqVT7ymRcogQL41NMPc+9dd0WP1XEd\nfSflnP8d8EBK6U+BPwBOdPL5tL3Gxic5PznD1Mwi5ydnGBufjB6pbafPTPDS1KvML63w0tSrnD4z\nET3SbTn57AvMX12lWcL81VVOPvtC9EhtWw8vQLl2vBN0+p7ve4Fv55zfnVI6CPwm8NBGnzM8PNjJ\nke6oqu8yPbdEURQAFEXB9NxSZXeaubLcuqwCKFrHVd0FYHF59S8dV3Wf8q84ruouW9Hp2w5HgP8A\nkHO+kFLan1Iqcs7f/e/7NbOz8x0e6c4YHh6s/C4jQ/1cvLRAURSUZcnIUH9ld9q3ZzfziyutAJet\n46ruAjCwu8H81dWbjqu6z9pLctNxVXdZdysnj07fwPsm8EMAKaUHgPmNwqvucvzoKA+N7uO+fQM8\nNLqP40dHo0dq24ljh3jzffcw2N/Hm++7hxPHDkWPdFs+/tTDDN7VoKeAwbta93yr6lNPP3zjNyV8\n6unq7rIVRVl2roVrv2r2OeANwC7gF3LO39jgU8qqn/HW1eHKd527dK867VOzXYrNHtPR2w4550Xg\npzv5HJJURdX8vSFJqjjjK0kBjK8kBTC+khTA+EpSAOMrSQGMryQFML6SFMD4SlIA4ytJAYyvJAUw\nvpIUwPhKUgDjK0kBjK8kBTC+khTA+EpSAOMrSQGMryQFML6SFMD4SlIA4ytJAYyvJAUwvpIUoBE9\nQN2sNpuMjU8yPbfEyFA/x4+O0ujxHBdt5fp1Tp+ZYObKMvv27ObEsUP07doVPVbblq5d4+SzL7C4\nvMrA7gYff+ph+nt7o8dqy059z9R/wztsbHyS85MzTM0scn5yhrHxyeiRBJw+M8FLU68yv7TCS1Ov\ncvrMRPRIt+Xksy8wf3WVZgnzV1c5+ewL0SO1bae+Z4zvNnt5ZmHDY8W49OdXNzyumsXl1Q2Pq2Sn\nvmeM7zZ74767NzxWjDfce9eGx1UzsLux4XGV7NT3TFGWZfQMNypnZ+ejZ7gtdbx/NTw8SNVfF+/5\ndq+avmeKzR5jfDukDsFa5y7dq0771GyXTeNb7dOLJFWU8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleS\nAhhfSQpgfCUpgPGVpADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleSAhhfSQpg\nfCUpgPGVpACNTj9BSum9wEeAa8AzOeff6/RzSlK36+iVb0rpXuAZ4DDwKPBYJ59Pkqqi01e+PwJ8\nPee8BCwBT3b4+cI1y5JzF6a5vLjC3oE+jhwcoacoosdqS512WW02GRufZHpuiZGhfo4fHaXRU927\nbnV6beq0y1Z0Or5vAgZSSr8D7AE+lnM+2+HnDHXuwjRnJ6bobfRwbbUJwCMP7g+eqj112mVsfJLz\nkzMURcHFSwsAPP7ogeCp2len16ZOu2xFp+NbAPcC7wG+H/hPwAMbfcLw8GCHR+qsy4sr9DZaV1S9\njR4uL65Udqc67TI9t0SxdjVVFAXTc0uV3QXq9drUaZet6HR8LwHP55xL4FsppfmU0vfknF95vU+Y\nnZ3v8EidtXegj2urzdfO4nsH+iq7U512GRnq5+KlBYqioCxLRob6K7sL1Ou1qdMu627l5NHp+H4N\n+HxK6ZO0roAHNgpvHRw5OAJw0/2rqqrTLsePjgLcdM+3yur02tRpl60oyrLs6BOklD4EPAGUwD/P\nOX91g4eXVT/jrRseHqz82Xudu3SvOu1Ts102/Ylhx3/PN+f8WeCznX4eSaqS6v6ujSRVmPGVpADG\nV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleSAhhfSQpgfCUpgPGVpADGV5ICGF9J\nCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleSAhhfSQrQiB6gbpplybkL01xeXGHvQB9HDo7Q\nUxTRY7VltdlkbHyS6bklRob6OX50lEZPNc/XdXpdoF771GmXrTC+2+zchWnOTkzR2+jh2moTgEce\n3B88VXvGxic5PzlDURRcvLQAwOOPHgieqj11el2gXvvUaZetqOZlTBe7OLu44XGVvDyzsOFxldTp\ndYF67VOnXbbC+G6z+4cHNjyukjfuu3vD4yqp0+sC9dqnTrtshbcdttmRgyMAN92/qqrjR0cBbrrn\nW1V1el2gXvvUaZetKMqyjJ7hRuXs7Hz0DNtieHgQd+k+ddoF6rVPzXbZ9CeG3naQpADGV5ICGF9J\nCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleSAhhfSQpgfCUpgPGVpADGV5ICGF9JCmB8JSmA\n8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleSAnQ8viml3Smlb6aU3t/p55KkqrgTV76/CFy+A88jSZXR\n0fimlBIwCny1k88jSVXT2OwBKaUfzTn/fpv//E8DHwY+0ObnV06zLDl3YZrLiyvsHejjyMEReooi\neqy2uEv3Wm02GRufZHpuiZGhfo4fHaXR449wqmTT+AL/NKX0r4DfBj6Xc/72rfyDU0rvA57POX+7\ndQFMdb/St+DchWnOTkzR2+jh2moTgEce3B88VXvcpXuNjU9yfnKGoii4eGkBgMcfPRA8lbZi0/jm\nnI+mlIaAvw/8+lpIPw98Ked8fYNP/THg+1NKPw7cDyynlF7OOZ/d6PmGhwdvefhudHlxhd5G6wqk\nt9HD5cWVyu7kLt1rem6JYu3KvSgKpueWKr3PujrscKtu5cqXnPNcSunfAivAU8AJ4KMppSdyzn/0\nOp/zM+sfp5Q+CvzZZuEFmJ2dv6XBu9XegT6urTZfu8LaO9BX2Z3cpXuNDPVz8dICRVFQliUjQ/2V\n3gda4a36Dutu5SRyK/d83wZ8EPi7wJeAx3POf5JSehPwZeDQ7Y1ZL0cOjgDcdG+xqtylex0/Ogpw\n0z1fVUtRluWGD0gp/SHwG8AXc85/8V1/72TO+ePbOE9ZpzOfu3SfOu0C9dqnZrts+jOuW7nn+9YN\n/t52hleSdgx/N0WSAhhfSQpgfCUpgPGVpADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDx\nlaQAxleSAhhfSQpgfCUpgPGVpADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKUAjeoC6WW02\nGRufZHpuiZGhfo4fHaXRU81zXLMsOXdhmsuLK+wd6OPIwRF6iiJ6rLbU6XWBeu1Tp122wvhus7Hx\nSc5PzlAUBRcvLQDw+KMHgqdqz7kL05ydmKK30cO11SYAjzy4P3iq9tTpdYF67VOnXbai/qeXO+zl\nmYUNj6vk4uzihsdVUqfXBeq1T5122Qrju83euO/uDY+r5P7hgQ2Pq6ROrwvUa5867bIV3nbYZseP\njgLcdP+qqo4cHAG46Z5vVdXpdYF67VOnXbaiKMsyeoYblbOz89EzbIvh4UHcpfvUaReo1z4122XT\nn0x720GSAhhfSQpgfCUpgPGVpADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleS\nAhhfSQpgfCUpgPGVpADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKUCj00+QUvok8FZgF/CJ\nnPOXO/2cktTtOnrlm1J6O3Ag53wYeDfwK518Pkmqik5f+X4D+OO1j68A/SmlIudcdvh5wzTLknMX\nprm8uMLegT6OHByhpyiix2rLyvXrnD4zwcyVZfbt2c2JY4fo27Ureqy2LK+ucuq5F7mysMKeu/s4\n9cQPsrvR8W/8Oma12WRsfJLpuSVGhvo5fnSURk817yLW6T2zFR396luL7NW1wyeA8TqHF+DchWnO\nTkzR2+jh2moTgEce3B88VXtOn5ngpalXoYD5xRVOn5ng59/3A9FjteXUcy8yc2UZgJkry5x67kU+\n8eTh4KnaNzY+yfnJGYqi4OKlBQAef/RA8FTtqdN7ZivuyKk/pfQY8EHg72322OHhwc4P1EGXF1fo\nbbSuQHobPVxeXKnsTjNXlmH9AqRoHVd1lysLK3/puKq7AEzPLVGsXR0WRcH03FJl96nTe2Yr7sQP\n3N4FnATelXOe3+zxs7ObPqSr7R3o49pq87Wz+N6BvsrutG/PbuYXV1oBLlvHVd1lz919r135rh9X\ndReAkaF+Ll5aoCgKyrJkZKi/svvU6T2z7lZOHkVZdu4uQErpHuC/AO/MOb9yC59SVv1fep3uX3nP\nt3t5z7e7DQ8PbrpAp+P7IeCjwP/itesn3p9zvvg6n1L5+K4bHh6s/Nl7nbt0rzrtU7NdNo1vp3/g\n9lngs518Dkmqomp+nyJJFWd8JSmA8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleSAhhfSQpgfCUpgPGV\npADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDxlaQAxleSAhhfSQpgfCUpgPGVpACN6AHq\nZnl1lVPPvciVhRX23N3HqSd+kN2Nav5rXm02GRufZHpuiZGhfo4fHaXRU83zdbMsOXdhmsuLK+wd\n6OPIwRF6iiJ6rLbVaZ86fZ1tRTWr0MVOPfciM1eWAZi5ssyp517kE08eDp6qPWPjk5yfnKEoCi5e\nWgDg8UcPBE/VnnMXpjk7MUVvo4drq00AHnlwf/BU7avTPnX6OtuK+p9e7rC5+ZUNj6vk5ZmFDY+r\n5OLs4obHVVOnfer0dbYVxnebDQ32bXhcJW/cd/eGx1Vy//DAhsdVU6d96vR1thVFWZbRM9yonJ2d\nj57htnjPtzvV6R4p1GufOn2drRseHtz0xTC+HTI8PIi7dJ867QL12qdmu2wa32qfXiSpooyvJAUw\nvpIUwPhKUgDjK0kBjK8kBTC+khTA+EpSAOMrSQGMryQFML6SFMD4SlIA4ytJAYyvJAUwvpIUwPhK\nUgDjK0kBjK8kBTC+khTA+EpSAOMrSQGMryQFML6SFKDR6SdIKf0y8ENAE/jZnPN/7fRzSlK36+iV\nb0rpbcDfyDkfBp4AfrWTzydJVdHpK993Al8ByDlPppT2pJTuzjkvdPh5w1xZXubEZ56nWUJPAaef\nPsye3bujx2pLnXZZuX6d02cmmLmyzL49uzlx7BB9u3ZFj9W2Zlly7sI0lxdX2DvQx5GDI/QURfRY\nbVltNhkbn2R6bomRoX6OHx2l0VP/O6Kd3vB7gdkbjl9Z+2u1tR4rgGbZOq6qOu1y+swEL029yvzS\nCi9NvcrpMxPRI92WcxemOTsxxf/81mXOTkxx7sJ09EhtGxuf5PzkDFMzi5yfnGFsfDJ6pDui4/d8\nv8ump+bh4cE7MUfHrMfqxuOq7lSnXWauLP//r76idVzVXQAuL67Q22hdO/U2eri8uFLZfabnlijW\nrtqLomB6bqmyu2xFp+P7HW6+0t0PbHiKnp2d7+hAndZT3BytnqK6O9Vpl317djO/uNIKcNk6ruou\nAHsH+ri22qS30cO11SZ7B/oqu8/IUD8XLy1QFAVlWTIy1F/ZXdbdysmj07cdvgb8JEBK6S3AVM55\nscPPGer004fpWbvCWr9PWlV12uXEsUO8+b57GOzv48333cOJY4eiR7otRw6O8I5D93Hgr+/lHYfu\n48jBkeiR2nb86CgPje7jvn0DPDS6j+NHR6NHuiOKsiw3f9RtSCn9S+CHgevAh3PO/32Dh5dVP+Ot\nGx4erPzZe527dK867VOzXTa9xdrxe74555/v9HNIUtXU//c5JKkLGV9JCmB8JSmA8ZWkAMZXkgIY\nX0kKYHwlKYDxlaQAxleSAhhfSQpgfCUpgPGVpADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwl\nKYDxlaQAxleSAhhfSQpgfCUpgPGVpADGV5ICGF9JClCUZRk9gyTtOF75SlIA4ytJAYyvJAUwvpIU\nwPhKUgDjK0kBjK8kBWhED7AupfS3gK8Av5xzfjZ6ntuRUvok8FZgF/CJnPOXg0dqS0rpLuDfAG8A\n/hrwL3LOXw0d6jallHYD/wP4pZzzb0XP046U0g8DX6S1RwFcyDk/HTtV+1JK7wU+AlwDnsk5/17w\nSG1LKf1j4H1ASeu1+Ts553v+qsd2RXxTSv3ArwL/MXqW25VSejtwIOd8OKV0LzABVDK+wI8D53PO\np1NK3wd8Hah0fIFfBC5HD7EN/iDn/A+jh7hda++RZ4BDwCDwMaCy8c05fw74HEBK6W3AT73eY7si\nvsAy8G7gn0UPsg2+Afzx2sdXgP6UUpFzrtx/Sphz/sINh98HvBw1y3ZIKSVglOqfQKB1VVUHPwJ8\nPee8BCwBTwbPs52eAY693t/sivjmnJvAX7TeG9W2Ftmra4dPAONVDO+NUkrngPuAR6NnuU2fBj4M\nfCB4ju1wIKX0FeBeWrdQqvpd45uAgZTS7wB7gI/lnM/GjnT7Uko/APzvnPPM6z3GH7h1SErpMeCD\nwD+JnuV25ZyPAI8Bvx09S7tSSu8Dns85f3vtL1X5yvFPgVM55/fQOpH8ZkqpKy6k2lDQOoG8h9b7\n5fOx42ybJ2j9vOR1Gd8OSCm9CzgJ/GjOeT56nnallN6SUrofIOf834BGSul7gsdq148Bj6WUXqD1\nxviFlNI7gmdqS875OznnL659/C3g/9D6zqSKLtE6KZZru8xX+GvsRm8Hnt/oAd14tqzyFQkppXuA\nTwLvzDn/3+h5btPbgAeAn0spvQEYyDm/EjxTW3LOP7P+cUrpo8CfVfXb25TSMWAk5/zplNL3AvuA\nqeCx2vU14PNrvyF0LxX+GluXUhoB5nPOqxs9rivim1J6C637cQ8A11JKPwH8g5zzldjJ2vLTwF7g\nCymlgtavnLw/53wxdqy2/Aatb2n/M7AbeCp4HrX8LnBm7dZWL/DkZm/0bpVz/k5K6d8Df0TrvVL5\n23TACPC693rX+f/zlaQA3vOVpADGV5ICGF9JCmB8JSmA8ZWkAMZXkgIYX0kKYHwlKYDx1Y6QUvq5\nlNK/Xvs4pZT+JKU0ED2Xdi7jq53iV4C/mVI6DPwa8KGc82LwTNrBjK92hLX/p/LjwBdo/bE7fxg8\nknY446udZC8wT+tP5ZBCGV/tCGt/cOav0/pz6VZSSv8oeCTtcMZXO8XHgC/lnL8J/CxwKqW0P3gm\n7WD+LyUlKYBXvpIUwPhKUgDjK0kBjK8kBTC+khTA+EpSAOMrSQH+H588tQBeuT4hAAAAAElFTkSu\nQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.lmplot('x', 'y', data=d, fit_reg=False)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "lmplotのx_jitter, y_jitterに揺れ（jitter)の量をセットすると、重なっている点がずれてプロットされます。" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAFgCAYAAAAcmXr5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHglJREFUeJzt3X1wHPd93/H3Hg4giMPRAuEVCQqiSKT1T/AkoORUlmua\nrivacWI5Y9Wd1pm4flDtznjiTJvMNJlmmthWptNmPHWm02mbTDyx48yk0zoZ20nHSmuFTFUGsmjV\nogk5Jn+QwscDAeoIkuLh8HA43PYP4K4HEE8EcPjuHj6vf4TF3XG/33347G9/CwhBFEWIiMj2SlkX\nICKyEyl8RUQMKHxFRAwofEVEDCh8RUQMpK0LqFcuz0W3bk1al7Eluro6UC/x00y9QHP100y9hGE2\nWOs9sRr5ptMt1iVsGfUST83UCzRXP83Uy3rEKnxFRHYKha+IiAGFr4iIAYWviIgBha+IiAGFr4iI\nAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGF\nr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIiBdCP/cedcBvhDoAtoA37Te/+d\nRq5TRCQJGhq+wCeB8977f+2c6wFOAv0NXqesUyWKGBwa5Wp+gqnpMrt3pXnw/k6ODvSQCgLr8hqi\nUok4dfYauXyR3jATy16r+2WlGutf7+/rZuBwV+x6qKpEEX919hrfO/86AG/v38e7YrjNLTQ6fG8A\nP7Hw9V4g3+D1yT0YHBrl5JkRJiZnKUyWyHa08erIGwAcO3LAuLp7t1ZoAZx46Qonz4wAMJy7Dcz3\nup7Pbpf6/fLij8YYvnqbp5/sr9UzODTKX3z/KjfvzPD82REO3t/Jr370baRT8ZtFHBwa5X+8cJnC\nZAmA6zenCLh7mz8QZsh27uL8xZsrbv847aOt0NDw9d7/d+fcJ51zrwL3AU82cn1yb3L5IgCl8lzd\nf1tr30+aamjB4mCtd2nszqLlaq/r+ex2yeWLtQsiwNCFcQaHRmv15PJFbt6ZYXK6DAFcGC3wtWfP\n86kPvtWk3tXk8kVK5TmiKKISwZ3JEi+eu04URXzv/OtcvzlFZ0crLw/nSbUEVOaiuy441dA9fe46\n129OkdmdNt9HW6HRc74fBS5773/GOTcA/D7w2GqfCcNsI0vaVnHvpb+vm4tjd9i9K01ptsLuXWla\n0yn6+7rvqj3uvQCMF0u0plOLlpfWfWj/Hn50Yby2XO11PZ/dLv193Xzv/HWChVHd7l3pRfX093Xz\n/NkRWBj0pQIYvTUZy33U39fN2b+5wdRMmUoUEUSQe32C8TemmZmdY3pmjpaWgHKlQmlmjkpl/nM/\nvHSToYu3eN/jD/Hc6cucemWU/O2p2vv3ZNpM99FWaPS0w1HgfwF474eccwecc4H3PlrpA/l8ocEl\nbY8wzMa+l4HDXRQK03fN+Q4c7lpUexJ6AejOtDFbrixaXlr38ccOUihM125dq72u57PbZeBwFz9+\naC9DF8ZpS7fQ3tayqJ6Bw10cvL+TC6MFUgGkgoCero5Y7qOBw13cefwgf/bCJSany2TaWymV55ia\nKdOWbiGKyrWvS8yPkAHSqRTnLozzSN9ezl0YZ7ZcIZ1K1d6/e1fadB+tZT0XhUaH72vAO4BvOuce\nAgqrBa9sr1QQJPq2bamjAz0Ai+YEl0qllu95PZ/dLqkg4Okn+++a36x//Vc/+ja+9ux5Rm9N0tPV\nwSc+8LBZvatJBQHvfuQBgiCoTetMTM6/ltk9Hz/79u7m7f37uPL6BC+du05buoXM7jS9YQaA3jDD\ncO42nR2ttfc/3r/PdB9thaB6pWmEhR81+wqwD2gBft17//wqH4nieiW7V0kZLa6HeomvpPSz9OEa\nUcTIjclFD866uzv51snhux6oJfFBWxhm1yyw0Q/cisBHGrkOEYm/9dxlrXRX0mx3aFXx+9kUEZEd\nQOErImJA4SsiYkDhKyJiQOErImJA4SsiYkDhKyJiQOErImJA4SsiYkDhKyJiQOErImJA4SsiYkDh\nKyJiQOErImJA4SsiYkDhKyJiQOErImJA4SsiYkDhKyJiQOErImJA4SsiYkDhKyJioKF/Ol4kKSpR\nxODQKLl8kd4ww9GBHlJBYF2WNDGF7wboRG0+g0OjnDwzAsBw7jYAx44csCxpXSqViFNnr3H19Qmm\nZsrsbk/zYNgZu2NS58zdFL7rVH/wTE7PcjU/QRAEiTpR5W7V/Xri+zmK02U6O1oByOWLpvVczU8w\nNV1m9640D96/cpieeOkKJ8+MMDE5S2GyRLajjVdzbwBwdKAnNoFXvbhFUcTLw3lOn7vO4/37dnQI\nK3zXqX5kdPPONG3pFvMTVTZvcGiUEy/nGL8zzdRMmelSme43tdMbZszquStMR+bDdLkL/KWxOwCU\nynN1/20lly/GajRfPUeKU2UKkyVK5TmK02XTmqzpgds61QdsW7qldrADZieqbF4uX6Q4VWa2XCEI\nAmbLldptu1U9sDRMV77AH9q/B5g/Juv/2xtm7vqM5SCheo5U+6nWuZMHLhr5rlNvmKmNHjo7Wul9\n8310tLfWbuckmXrDDC/+aAyAllRAtqONjvZWs1vh6nHWlm5hpjS3KEyXc/yxgxQK08vO+Q4OjdaO\n2dX+je1QPUdOn7vO9ZtTtbvGnTxwUfiuU/XgicP8mWydowM9DF+9zdCF8dpUUhxCark53+WkUsGK\nt+3LHbNWUsF8ncvNQ+9UQRRF1jXUi/L5gnUNWyIMs6iX+FmulyQ/iW/2fZNUYZhd8wDSyFd2vOqo\nTGQ76YGbiIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+I\niAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIgYb/9WLn\n3EeBXwFmgc957/+80esUEYm7hoavc24v8DngUSALPAPsqPCtRBGDQ6Pk8kV6wwxHB3pIBYF1WSJi\nrNEj3/cCz3nvJ4FJ4DMNXl/sDA6NcvLMCADDudsAHDtywLKkHaFcqfC1Z89z9fUJHry/k0984GHS\nKc2ySXw0OnwPARnn3J8C9wHPeO9PNnidsZLLF1ddlsb42rPneen86wCM3ZwE4FMffKtlSVum/m6q\nv6+bgcNduptKoEaHbwDsBZ4CDgN/CTy02gfCMNvgkrZPGGbp7+vm4tid2vf6+7oT2WMYZqlUIk68\ndIVLY3c4tH8Pxx87SCoVz5N+9NYkQV0gjd6arG33JG7/es+dvsypV0YBasfW+x5f9bSKvUol4rnT\nlxNxbG2VRofvdeAF730EXHDOFZxzb/be31jpA/l8ocElbY8wzJLPFxg43EWhMF2b8x043JW4Hqu9\nnDp7rTaFcnY4T6EwHdsplJ6uDnLXJxYt5/OFWi9Jdu7COLPlCgCt6RTnLozzSN9e46o259TZa5x6\nZZTZciX2x9Z6rOcC3+jw/Q7wVefcF5kfAWdWC95mlAqCRB9E9ZI0hfKJDzwMsGjOt1n0hpna84Pq\nctIl6djaKg0NX+/9NefcnwAvAhHwi41cnzRWkk76dCrVNHO8Sx0d6AFYNOebdL1hZtH0XJyPra0S\nRFFkXUO9KOm3hFXNcHtbVe2lGX5srpn2CzRPP5UoYujiLc5dGE/ssVUvDLNrFt/wX7KQ5tFMUygS\nL6kg4H2PP5T4uet7oR98FBExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BUR\nMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg\n8BURMaDwFRExoD8dLyLbrhJFDA6NkssX6Q0zHB3osS5p2yl8RZrUcgGXCgLrsgAYHBrl5JkRAIZz\ntwH48Hv3WJa07RS+2yyuJ0R9XQ+8uQOCgJGFGp964i3W5a1bXLevheUC7tiRA5Yl1eTyxVWXdwKF\n7zaL6wlRX9fLw3kAOjtaGc7dJptt55G+vZblrVtct6+FOAdcb5ip7Z/q8k6j8N1mcT0h6usolecW\nvmoF4NLYncSEb1y3r4U4B1x1jldzvrJt4npC1NfVlm5Z9Nqh/cmZi4vr9rUQ54BLBcGOvSOpUvhu\ns7ieEPV1LZ3zPf7YQcbHJ4wrXJ+4bl8LCrh4U/hus7ieEKvVlUol54FVXLevyFL6JQsREQMKXxER\nAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMK\nXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDDQ9f51y7c+4159zHG70uEZGk\n2I4/Hf8bwPg2rEdkx6tEEYNDo+TyRXrDDEcHekgFgXVZsoyGhq9zzgEPA99u5HqsLXfAi52dHECD\nQ6OcPDMCwHDuNgDHjhywLElWsGb4Oud+2nv/Pzf4738J+CzwyQ1+PhGWO+A//N49liVtmzgGXTME\n0Ea3ay5fXHVZ4mM9I99/7pz7T8AfAV/x3l9ezz/snPsY8IL3/vL8AJimHXrs5AM+jkHXDPtjo9u1\nN8zU3l9dlnhaM3y99x9wznUB/wD4nYUg/SrwDe/93CoffRI47Jz7WaAXmHbOXfXen1xtfWGYXXfx\ncdHf183FsTuLliGZvaxkpV7GiyVa06lFy9Z9L7c/6muyrm897mW71n//qSfeQjbbzqWxOxzav4fj\njx0klUrOuCcJ+2arrGvO13t/yzn334AS8AvAvwQ+75z7tPf+xRU+83PVr51znwcurhW8APl8YV2F\nx8nA4S4KhenaLeLA4S4gmb0sJwyzK/bSnWljtlxZtGzd93L7o1rTar3EyXq363L9PNK3l0f69gIw\nPj7R2EK3UFL2zXqs5yKynjnfdwNPA38f+AbwKe/9OefcIeCbwKObKzP5UkFgfqttpfpwMU4PG5th\nf8Rxu8rWWs/I998Cvwt8xns/U/2m9/6Sc+7r61mJ9/6ZDdYnMdcMQRdH2q7Nbz1zvu9a5bV/t7Xl\niIjsDPr1YhERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpf\nEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxER\nA2v+6XhpjEoUMTg0Si5fpDfMcHSgh1QQWJd1z5qlj7jYiduz2vN4sUR3pm1H9AwKXzODQ6OcPDMC\nwHDuNgDHjhywLGlDmqWPuFi6PSMggKYO42rPrekUs+UKsDOOIYXvJm10pJLLF1ddToq49rHcfgF4\n7vRlzl0Yj22QLd1+3zt3neJ0GWjei9tqx1Az3wkofDdpoyO/3jBTe391OYni2sdy+wXg1CujzJYr\nsQ2ypdtzqbhc3LbSasdQM99ZKXw3aaMjv+pIbOnILGni2sd69kscg2zp9oyiiL/8wbXa63G5uG2l\nas/1c75Vcb2z2goK303a6MgvFQRNcQWPax8r7ZeLY3fu+l6cLN2elSgiCILYXdy2UrXnMMySzxcW\nvRbXO6utoPDdpLiO/Ha6lfZLNtu+aM437uJ6cdsuzXx+BVEUWddQL1p65Uuq5a7iSaVe4quZ+mmy\nXtZ8KqhfshARMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDw\nFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BUR\nMZBu9Aqcc18E3gW0AL/lvf9mo9cpsppKFDE4NEouX6Q3zHB0oIdUsOZf+o6lSiXi1NlrTdHLTtPQ\n8HXOvQd4q/f+nc65vcAZYMeEb1Od5Mv0AiSyv8GhUU6eGQFgOHcbgGNHDgDJ22cnXrqyYi9JUd3m\n48US3Zm22G/zrdLoke/zwOmFr28DHc65wHsfNXi9sbDaSZ4E9SfF+M1JruYnCIKg1gsQm/7uJTRz\n+eKKy3HdZyv1d2nszqL3Le0tCarbvDWdYrZcAeKxzRutoeG7ELJTC4ufBp7dKcELq5/kSVB/Uly/\nOUlbuoXOjlZg+V6s+qtEEV/99jmGLozTlm5ZMzR7w8yiC0hvmKl9Hdd9ttJF4dD+PZwdztfeV99L\nUsR1mzdaw+d8AZxzHwKeBn5qrfeGYbbxBW2T/r5uLtaNTPr7uhPV33ixRGt6/pns7l1pZmbnasv9\nfd0AsejvudOX+eGlm5RmK5RmK7S0BIwXSyvW8tQTbyGbbefS2B0O7d/D8ccOkkrNj5Ljus/q90V1\nOQyzHO/uBFi2l6So3+at6VRstnmjbccDt/cDvwa833tfWOv9+fyab0mEMMwycLiLQmG6dqs4cLgr\nUf11Z9qYLVdoTafYvSvN337gTXS0t9Z6AWLR37kL46RTKaKoDMDUTJnuTNuytYRhlvHxCR7p28sj\nfXsBGB+fqL0e131W3Rf1y/l8gTDMrthLUlS3eXXONy7bfDPWc/Fo9AO3PcAXgePe+zcaua44SgVB\noueuqg/VVnsQEof+esMM/uotAErlOQb6umu136u47rNqP0sfeDaD6jYPw2ziQ/deNHrk+xGgG/i6\ncy4AIuDj3vtcg9crWyApJ8VywdRsT8vjelGQjWv0A7cvA19u5DpEFEySRPoNNxERAwpfEREDCl8R\nEQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfERED\nCl8REQMKXxERAwpfEREDCl8REQMKXxERAwpfEREDCl8REQMKXxERAw390/EiIlWVKGJwaJRcvkhv\nmOHoQA+pILAuy4zCdxssPej+7k/s57uvjOkglB1lcGiUk2dGABjO3Qbg2JEDliWZUvjeo41cvQeH\nRjnxco7iVJkXfzTGqaFrlMoVgiAwPQjX6qX6+nixRHemrWkuEhqBba/q9j7x/RzF6TKZ3WmCICCX\nL1qXZkrhe482cvXO5YsUp8oUJksAXLk+Qbajjc6O1trrFtbqpfp6azrFbLly1+tJ1QwjsCRdGKvb\nuzj9/8+Bzo5WesOMcWW2FL73aGlQric4e8MML/5orLbcmk5RKs8BrbXXLazVy0Z6TYJm6CtJF8bq\n9q0ONjLtaZ549AGODvRYlmVOP+1wj5YG5XqC8+hADwN93exqayHb0cab39TOQF83b+m9z/QgXKuX\njfSaBM3QV5IuIPXbt7OjleM/2cuxIwdiO1LfLhr53qNqUNbPF64lFQQ8/WR/7OYZ1+qlulx/a9sM\nNrIP46Y3zNSmTKrLcdUM27sRgiiKrGuoF+XzBesatkQYZlEv8dMsvSRpzne9mmXfAIRhds2doZGv\nSAKlgoBjRw40VWDtNJrzFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BUR\nMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg8BURMaDwFRExoPAVETGg\n8BURMdDwv17snPtt4B1ABfgl7/3/bfQ6RUTirqHh65x7N/C3vPfvdM49DHwFeGcj12mhEkUMDo2S\nyxd5IMxAFHFzcpbuTBtHB3pIBYF1iSISM40e+R4HvgXgvT/vnLvPOdfpvZ9o8Hq3XH3A9oaZRaE6\nODTKyTMjALw8/DrTpTnS6RTlcoXv/vUY3Xva2d2e5sGwM3ZhvFpfm3mvbJ+dsl+arc9Gh+9+oH6a\n4cbC915r8Hq3XH3ADuduA3DsyAEAcvli7X3F6TKl2TlSswFzcxGvjbzBhWt3yHa08WruDSKAKOJ7\n518H4O39+3iX4UG0tK8oiiAIOP3XY1wbn6QlFfCou5+PPPFjfPeVsRW3QZzc60ma9JN6tWMziVba\nH83WZ8PnfJdY84gOw+x21HHPxoslWtOpRcvVWvv7urk4dgeAIAhIpQKiCAggWvheuVKhNZ3iB6+N\nc+3GBG9MlAC48cY0e7LtvO/xh7a7pVof9X394MI41/JFbtyeohLNf+/5MzkAOna3MjVTZmZ2jl2t\nLdyYmN8GlUrEiZeucGnsDof27+H4YwdJpezC67nTlzn1yigAF8fukF2yfZceY9958RLffvEyM7Nz\nnP2bFjo7d/FT7zi0nSVvymrHZtKEYXbF/ddMfULjw/ca8yPdqgPA6GofyOcLDS1oo7ozbcyWK4uW\nq7UOHO6iUJgmly+y777dnLt8i6nSHKXZOdpbW5gtV0inUsyWK5Rmy0zNlOdHmMDUTJlzF8Z5pG9v\nLPoqleaYminXghcgiuDVK7d48P5ObhdmAJiaLnPz1iT5fIFTZ6/VRiRnh/MUCtOmI5JzF8YX9VS/\nfcMwe9cxduKlK4v6OvHSFR79se7tK3gTwjC76rGZJNV9s9L+S1Kf67koNDp8vwN8Afiyc+5twIj3\nvrj6R+Lp6EAPwKJboapUENTCpnrLdGOixM1bk7S3tTBdmqvN+UZRxPWbU8yU5gBoS7fQG2a2v6EF\nS/uKgOs3pwiCWRauDwQBPHh/J7t3pcl2tFEqz9GWbmH3rnTts/WWLm+33jBTuy2tLjez1Y7NJFpp\n/zVbnw0NX+/9d51z33fODQJzwGcbub5Gqg/Y9bxvuREWzIczsGjO1/IgWtpXtb6V5nxfHXkDaAXm\nAxniF3b3epK+vX8f129O1S4qb+/ftx1lbpn1HptJsdL+a7Y+g+rtb0xEcb2NuFcrhW8SVXtZ6UFI\nkh5YLbdfklT/Us14nDWDMMyueQBt9wM3SbCVRh5JH5EkvX5JJv16sYiIAYWviIgBha+IiAGFr4iI\nAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGF\nr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIAYWviIgBha+IiAGFr4iIgSCKIusaRER2\nHI18RUQMKHxFRAwofEVEDCh8RUQMKHxFRAwofEVEDCh8RUQMpK0LqHLO/TjwLeC3vff/xbqezXDO\nfRF4F9AC/Jb3/pvGJW2Ic2438AfAPmAX8G+89982LWqTnHPtwA+B3/Te/6F1PRvhnPt7wB8z30cA\nDHnv/4VtVRvnnPso8CvALPA57/2fG5e0Yc65fwp8DIiY3zc/6b3fs9x7YxG+zrkO4D8Cf2Fdy2Y5\n594DvNV7/07n3F7gDJDI8AV+FnjJe//vnXMHgeeARIcv8BvAuHURW+B/e+//sXURm7VwjnwOeBTI\nAs8AiQ1f7/1XgK8AOOfeDfyjld4bi/AFpoGfAf6VdSFb4Hng9MLXt4EO51zgvU/crxJ6779et3gQ\nuGpVy1ZwzjngYZJ/AYH5UVUzeC/wnPd+EpgEPmNcz1b6HPDzK70Yi/D13leAmflzI9kWQnZqYfHT\nwLNJDN56zrlB4AHgg9a1bNKXgM8CnzSuYyu81Tn3LWAv81MoSb1rPARknHN/CtwHPOO9P2lb0uY5\n5/4OcMV7//pK79EDtwZxzn0IeBr4RetaNst7fxT4EPBH1rVslHPuY8AL3vvLC99K8sjxVeAL3vun\nmL+Q/L5zLhYDqQ0ImL+APMX8+fJV23K2zKeZf16yIoVvAzjn3g/8GvDT3vuCdT0b5Zx7m3OuF8B7\nfxZIO+febFzWRj0JfMg5913mT4xfd849YVzThnjvr3nv/3jh6wvAGPN3Jkl0nfmLYrTQSyHBx1i9\n9wAvrPaGOF4tkzwiwTm3B/gicNx7/4Z1PZv0buAh4Jedc/uAjPf+hnFNG+K9/7nq1865zwMXk3p7\n65z7eaDHe/8l59x+4H5gxLisjfoO8NWFnxDaS4KPsSrnXA9Q8N6XV3tfLMLXOfc25ufjHgJmnXP/\nEPiw9/62bWUb8hGgG/i6cy5g/kdOPu69z9mWtSG/y/wt7f8B2oFfMK5H5v0Z8F8XprZagc+sdaLH\nlff+mnPuT4AXmT9XEj9NB/QAK871Vun/5ysiYkBzviIiBhS+IiIGFL4iIgYUviIiBhS+IiIGFL4i\nIgYUviIiBhS+IiIGFL6yIzjnftk593sLXzvn3DnnXMa6Ltm5FL6yU/wH4C3OuXcC/xn4Z977onFN\nsoMpfGVHWPh/Kn8K+Drzf3bnr4xLkh1O4Ss7STdQYP6vcoiYUvjKjrDwhzN/h/m/S1dyzv0T45Jk\nh1P4yk7xDPAN7/1rwC8BX3DOHTCuSXYw/S8lRUQMaOQrImJA4SsiYkDhKyJiQOErImJA4SsiYkDh\nKyJiQOErImLg/wGi/6LDPmbzZAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# jitterを使って同じポイントに重ならないようにする\n", "sns.lmplot('x', 'y', data=d, fit_reg=False, x_jitter=0.2, y_jitter=0.05 )\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

GLMを使った分析

\n", "\t

\n", "\t\tこのデータを６章と同じようにGLMで分析してみます。\n", "\t\t予測値は直線となっており、６章のようなロジスティック曲線にはなっていません。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Rにdを渡す\n", "PandaDf2RDf(d, \"d\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "\n", "Call:\n", "glm(formula = cbind(y, N - y) ~ x, family = binomial, data = d)\n", "\n", "Deviance Residuals: \n", " Min 1Q Median 3Q Max \n", "-4.4736 -2.1182 -0.5505 2.3097 4.0966 \n", "\n", "Coefficients:\n", " Estimate Std. Error z value Pr(>|z|) \n", "(Intercept) -2.1487 0.2372 -9.057 <2e-16 ***\n", "x 0.5104 0.0556 9.179 <2e-16 ***\n", "---\n", "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n", "\n", "(Dispersion parameter for binomial family taken to be 1)\n", "\n", " Null deviance: 607.42 on 99 degrees of freedom\n", "Residual deviance: 513.84 on 98 degrees of freedom\n", "AIC: 649.61\n", "\n", "Number of Fisher Scoring iterations: 4\n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 通常のGLMを使って解析\n", "r('fit <- glm(cbind(y, N - y) ~ x , data=d, family=binomial)')\n", "r('summary(fit)')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAAAAACDQ6CjAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC\nAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAJiS0dE\nAACqjSMyAAAACXBIWXMAAABIAAAASABGyWs+AAAfqElEQVR42u2dT2jcaJrG9ZXKViUZVSZdmWQz\njZxMt2dbNvYExtVkN9M9nt4VhAxuL72HYqkmZKEZesUwdMOaXWJYUtCQw0CaYZtcdCv21rccBM0e\ncpxDyCEMOiw51M2Epmh8+fAhMNTqk/8ktlWq95VLn+vP88zkizuup9569Cup9KkkvUYPmkoZp/0C\noNMRwE+pAH5KBfBTqgHgX2Zo+/uXbP3At0i+5Xv+K/t+W0uYH/J4+JZMNJICPsqQ7ERsdfmWHt/S\n4b+yjtQSppvHw7dkonkJ8G9apJYwAE8XwPMtAE+3SC1hAJ4ugOdbAJ5ukVrCADxdAM+3ADzdIrWE\nAXi6AJ5vAXi6RWoJA/B0ATzfAvB0i9QSBuDpAni+BeDpFqklDMDTBfB8C8DTLVJLGICnC+D5ltzg\nn20+GwXwgWM7we5raNRsN6T6mODnTXFuruoERYbxa5WSmPnHparTqtu1RkFlgnOitBi1zgpDlMpm\nqfSjtEIZ4NuL7cVnpw/eF17gmd+pH12rEbgW1cgDXxHve8J47Jmt4sJ41ntGufyWUXvaEBdbLdst\npExLWHcWhS1E2TQMQ5TPlCophTLA+81e0z998JYi0bgcD6Gp1nanTjSywH9pBJHtz/w0XimLC2OG\n5uXInb0kup6t3r+cNxm9zJkz8eAbJbvanInX+Zb10Tnz+HYsA/zW4sriVq83Pz8vM/TXHcnWK86D\nKwfjvWV5MBK0w3llH56Vsix/HY+isDD3LklDym+EFK+WP1OBlu8VUUZ8pkZj9m9+vbxumBeWl9fL\nKYW2+4OPV3c/XuMfPXrUydDOyw5b25wHiyfx0K7EwwM1dG4tEI0vOa/sn0WnU35y/mznSaWwMHEI\n8e+dZuk/jO0b83G9ztwXRZQx5+PhiVGuLd+eM8TsrbmPy/bxQt/3B7/a7rWbp7+pd50wCq1V9aNd\njz/BTOruF+8zXlyO6hWjFdbqxYWpuZfEl+ayUekG4or6yC+kzLL4NF5gQpTMT+KPeGPevFBNKZSx\nqQ9Xmovh6YOPXNMWbvIaQtuyBXlfmAc+ECLeC64Kzg4XN0xYs+Iahimq5hUzzsKZQTDKXIlrlFpX\njH2V0gplTee22qMxjw+Dg3l8QJ7M8efxQRA9eVxwmDCMvgqj7rdh8mNhZb5Uzx0+/p8gVP9LLYQD\nOIcsUkuYUT+AA/AFhQF4ugCebwF4ukVqCQPwdAE83wLwdIvUEgbg6QJ4vgXg6RapJQzA0wXwfAvA\n0y1SSxiApwvg+RaAp1ukljAATxfA8y0AT7dILWEAni6A51sAnm6RWsIAPF0Az7cAPN0itYQBeLoA\nnm8BeLpFagkD8HQBPN8C8HSL1BIG4OkCeL4F4OkWqSUMwNMF8HwLwNMtUksYgKcL4PkWgKdbpJYw\nAE8XwPMtAE+3SC1hAJ4ugOdbAJ5ukVrCADxdAM+3ADzdIrWEAXi6AJ5vAXi6RWoJA/B0ATzfAvB0\ni9QSBuDpAni+BeDpFqklDMDTBfB8C8DTLVJLGICnC+D5FoCnW6SWMABPF8DzLQWBb7luPYrC5Tl6\ne6AThCga/LXK7DUFvu66fsFhAnfp5m9dp/GpWSrdKaxMeK1iLWT8Pjf4hun5rtUyP/jCJbcOyL2s\nigYfmuLdd4UZdio1Pw5WaBjf9DbOiN/6ZcN82zQWCyoTlsS7rpjp/4Dc4BPabtmPtydegf179lUs\n+LoIVKuC+o1L0X7Po8LCmK3o64rrhkK0oqhiFFTGiRPJL8x63weQwHePa+OqGr8W3Vfb3W6ly9MO\n8/GxenzL9jb5oWd+rMYfn1n+L/X31Y0iw1S73Q8/fFrdEL9YihejoNZililf6sZozGrfB/yQc40P\nkjYnDSOM1/jQZL7pR26NP58kMM8vfKb+tlit59hrfBStLbXsllAddjyDWotZZnZW7dyVnL4PyL2p\nt714o2hddGLwRfbv2Vex4H2xrJr4+A/K8VbeswsN47jRU7PSiITxVfzpUiqoTL30q0j+tNR/RzU3\n+NC2nHg3yCnPmUX279lXwXv1N0WpJH4VddZNx7IKbBYTJR2JlkxRq5VUXyJB3i/mlnlnN1FfnWAe\n3/LVEnr8Bbsb6wiCj8KPPgrVdC70C2wxuqdg42nkx+vinbfpszl+mfCjf3mS8WscwDlkkVrCjPcB\nHMqzDy8EwPMtAE+3SC1hAJ4ugOdbAJ5ukVrCADxdAM+3ADzdIrWEAXi6AJ5vAXi6RWoJA/B0ATzf\nAvB0i9QSBuDpAni+BeDpFqklDMDTBfB8C8DTLVJLGICnC+D5FoCnW6SWMABPF8DzLQBPt0gtYQCe\nLoDnWwCebpFawgA8XQDPtwA83SK1hAF4ugCebwF4ukVqCQPwdAE83wLwdIvUEgbg6QJ4vgXg6Rap\nJQzA0wXwfAvA0y1SSxiApwvg+RaAp1ukljAATxfA8y0AT7dILWEAni6A51sAnm6RWsIAPF0Az7cA\nPN0itYQBeLoAnm8BeLpFagkD8HQBPN8C8HSL1BIG4OkCeL4F4OkWqSXMqIPf2gwHgA9cd10D+LjM\nJ/wqfPDhjeV6qKpx+hHlBh/X8Qotswu+T5kM8FsrD1fameAbZr1RncsRnPdwVeYyt+9NDvCBuXDP\nNW/G1WyH7soLXqWq0VPlBd+vTAZ4f7O39SwTvGpIJC96/OC8h6syvRq7DBt8rRFv6peFugu/3Sgq\nTGLp7qWK6Knygjdb6WUywDdXVla3er3bt2/36WeTNCR69Z9LzG483MY6G5fioff1JW4VRjOiXYnu\n9qvuhlA/3qWHytFZaWfnoJsTOVWOMqpPVLLw0spkNCNqbvYeNrPAJ0/36l+LBv+1aqXUu3uVWyUf\n+N8l4NcKB5+k6tJT5QTft0wGeP9hAj5jU2/FG8TnFWZnRv5mS5V5YbHLsDf1rhtv6h3VDi406Q1K\n8m7qVaqQnirvpr5fmaydu8XmYvZnfMt03PKNHMF5D29ZcZk6uwobfFizl62aCsVpNpkXvEqV0RBw\nCGUS8P3KZM7j2wPn8b7/WMc83ve/41fJMY9/cE+t7r7P6a+Ufx7PqpN/Hp9eBgdwDlmkljCjfgAH\n4AsKA/B0ATzfAvB0i9QSBuDpAni+BeDpFqklDMDTBfB8C8DTLVJLGICnC+D5FoCnW6SWMABPF8Dz\nLY8BnmyRWsJoAe+ZZYAnW6SWMMWDD11TuM8BnmyRWsIUDT50hXBDfMYzLFJLmGLBB64wPXXiIMDT\nLVJLmCLBB46xfxoJwNMtUkuY4sC3HOP1CcIAT7dILWGKAu/bhv3G2XUAT7dILWGKAd+wDefQSZUA\nT7dILWGKAO+ZhnPk3DqAp1ukljDDB++Zwjl2SiXA0y1SS5ghg9+bth8TwNMtUkuYoYKPsZtp2AGe\nY5FawgwRfOCIvld/ADzdIrWEGRr4N47WpAjg6RapJcyQwPu1Q9P2YwJ4ukVqCTMU8EeO1qQI4OkW\nqSXMEMB7tuEMuqIX4OkWqSXMicGnHK1JEcDTLVJLmJOBD+umcCnX2QI83SK1hDkJ+H5Ha1IE8HSL\n1BImP3h1kkWdhh3gORapJUxe8NnT9mMCeLpFagmTD3yM3fY4FoCnW6SWMHnAbwycth8TwNMtUksY\nPnjfNmrs234BPN0itYThglfT9m/5ZQCebpFawvDA7x6twbVzZE0E+OSSmDBfGYCnW6SWMGTwbxyt\nAXiyxh78oZMsAJ6sMQd/5GgNwJM11uD9Ny+JyVkG4OkWqSXMQPApJ1kAPFljC76RdpLFyID39lrc\nFA8+9NzV+nX6DeT3xQffWviAc99qfpjG7kLLBN/nJAveMvNdr5jr40PL8Rxr4EXYQ1hWUWC6l43y\nKudrqV3xGxVYNz4w6c1o2GEOFlp/8KHX7yQL1jKzbc8VrSLAO3W1oFSvpsLB241AhO710OSujFzw\nvhVv6gOT+m03P4zjxkPd6Q8+6yQLzjJTNeI4JwffOSaRjOX4z87LDlvbnAdXOjdudNrnOzcWmFVe\nMl/ZQrPzcqcz90VhYcrJGC+67VTPkwVRXnhy8jKdjt1WY+UvWQvg+3zgDzIUD77cubWQgL/BrMIG\nf0sL+HI6+PacUb41lDIx+AdqHAL4PlstVw3Fb+q9ePNb07KpDzuyVeCmPllebuqmPnAMyxtSmXj/\nUDUY9OxCdu7sWr1m69m5s5y3xIyOnbu6eWPZZH/rzdi5219oR8EPuiSGVyaWY9cdMyhmOtfYa8Cq\nYR7fcO94q/x5Fn86F9xY567vvDC+m0waDoMffEkMt0w8MXUbuN0ZxyK1hHkTvGdRro0YoQM4pGcf\nXogJBU+7JCZnGYCnW6SWMHvgyZfE5CwD8HSL1BImAc+4JCZnGYCnW6SWMDH4g/uOFlgG4OkWqSVM\n91veJTE5ywA83SJ1hAmuMi+JyVcG4BkWWXyYwDGqG8WXiQCeY5FFhwldw/LRk4auyQCv7knnoRkR\nR5MAfg87wHM0AeDrQnh7FoAna+zBe+brwzUAT9eYg38TO8BzNNbgfdM4dHAW4OkaY/C+ffQbOICn\na2zBH8cO8ByNKfjANuzjX7wCPF1jCV7dpyrtpCqAp2sMwffDDvAcjR340DX6nrEL8HSNGfiDo7Pp\nFoAna6zAZ2MHeI7GCXxdiHq2BeDJGh/wh4/OplsAnqxxAU/ADvAcjQd43zIoZ0wDPF3jAD7t6Gy6\nBeDJGn3wZOwAz9Gogw9qaQfl+1kAnqzRBt//6Gy6BeDJGmXwTOwAz9Hogn+acVC+X36AJ2tUwQ86\nOpueH+DJGk3wCvtajvwAT9ZIgveEcE+1f3y2AJ5ukYwH7x2dBfhCQ4wc+IOD8gBfaIgRA+9bhhPm\nDgPwdI0U+ENHZwG+0BAjBP7IQXmALzTEyIBXB+UPNU0A+EJDjAj4lKOzAF9oiJEAn3pQHuALDTEC\n4NWZ8imNSyYS/MPB4B83G+w7PucGHzYa9LtYc8E3Gq0M8H0PynPD3Jmf/3YUwb94g3tzZSB4Tyy4\n7Hu85wXvm65L/0qEB75lOq61IPv8NuO7GGYYS1TK4iP+Aigc/DuV1X324WDwgfm8E7G7OuQEH5ot\nNVDXeR74pEHExfX0X6rDdMMJ44h4Wf1SFNvlbE/MTX1wvXz5E/XD1mpPgb99+3a3r5bWXm3H40aX\npR3ew5V68Z+Nq+qntSWiZXub8fxPK2r8UzXtd2tlsfR0SGHMnyuLuMleADmW2ausBfBDymf85xWj\nfL3XW/Pbi+FIgU+QFwQ+QZ4GPhs7G/yKspRGEPz9d8TlT168iFf6zWbzZ5sDNvUt1dyM3Sco96Ze\nfaRY1D0K5qZePe3FD47+szpMl/1BlmNT/8kobuqvf5J8xN//To2Dd+7cyq36/u3bigyRvAbPrHuW\nS7Vwd+5cz7bl4X+knDLNvKWpKS6fM27wF8DIzeO/Waiz3765p3NB3aX3mGVO58K62zg8naOdKc8N\n88tK9b9HcTqHAzh7UofpKO+yiTyAM73g6adMA3yhIfSC55wpD/CFhtAJnnfKNMAXGkIfeO6Z8gBf\naAhd4PkXSAB8oSE0gV+n3NDixGEAni4t4L0yGzvAFxxCA3jPFAvPtYQBeLoKB59cIHEK3aQLLAPw\ng+XtthIA+OkCf9DMG+CnCby6HGrvuxiAnx7wh76CA/hpAX/km1eAnw7wx75wB/hpAJ9yngXATz74\n1NNrAH7Swfup3aEAfsLBK+yp51kA/CSDzzi9BuAnF3zmWVUAP6ngB5xMB/CTCX7gOZQAP4ngCafO\nAvzkgSedMQ3wkwY+xk45hxLgJwt8SMMO8JMFnnHGNMBPDnjWifIAPyngmddHAPxkgGdfFgPwkwD+\nyQK7XQzAjz/4eG1n36IF4McevMK+8IRdBeDHHHxdqMtiTqd/PMUC8GRxwO+1iwH46QJ/0CUI4KcJ\nvPf6AneAnx7w3pv3NQD4aQHvHb6dBcBPB/i9K51fC+CnAfzBlc6vBfCTD/6NK51fC+AnHXyfmw0D\nfEHgQ9+PfP95seBVEaX+4PveY5oL/tP3Py4afJwmjJfZ06LBqzJRdO+LjFt45Qfvme6PjLfc8q0i\nQzRMxxHqHtL9wGfcWpwHPjDFrBDfFBlmN03JcctrhZZRfZVMr2Vem8v4djI3eN8KQ/G+FT0usttC\noHpfJL2O0sFn3lGeB362FEbhTImdhUFEpQlLcZ2nhXao2O3ZJERLdoL+vUNyg3cbke9EtZa8QW4Z\nwQ9Rr++VSgU/oJEAC3yr5MVjc4bbSY1DRKXxnThN92aBy0xhifV2RX3G1/sWIoFPa2aztKGaA13d\nePUxtTlQjsY6S2vJeHe3GdFhbVSN6tdZZlYzoo2ZuEj3TyazrRInTJJmYylOs/NPBS6zvZ5N751R\nzYj6d276Ie8a79WiQHxlhpJ8Y/cc796WFe11mjv6Glr24MtiWJv60tl4eKvEvqMpY1VUaQIRp+lW\ni9ywBKrfUXS2FMZrfP/OTfl37mq1xmXDaVTn2BkYIVy74Vn16Ch42mUxLPANMbNkiXWOhRsmcq2G\nZwqvUb1aaBnP9Bq2U7fWm7bT90EnmM55rrvuuveKnc75rpsQfvM1EPtHMKdzwYVK9YEsNEySJv5z\nt+DpXFwi3i3ylxcyltL4HcAhtw3BAZxJOnLH6BYD8JMDntMkCOAnBjwLO8BPCnjilc6vBfCTAP6F\ny8QO8JMAPnDYV0MB/PiDV/cj/IRfBeDHGnyoTqryT6GpcBFh9i0AP0iBK4SrvoEDeL5lfMH7akd+\n92sTgOdbxhW89+YXcADPt4wl+NA1xZtnWQA83zKG4NX07fDVEQDPt4wdeDV98478G8DzLeMFfm/6\ndlQAz7eME/iD6dtRATzfMj7g35i+HRXA8y3jAt7LOn8S4PmWsQB/dPp2VADPt4wB+OPTt6MCeL5l\n5MGnTd+OCuD5ltEG32f6dlQAz7eMMvi+07ejAni+ZXTBZ0zfjgrg+ZZRBe8NvvzttQCebxlJ8IOm\nb0cF8HzLCIIfPH07KoDnW0YOvF8dPH07KoDnW0YLvJq+XeVfHg7wfMsogd+dvp12//j+FskvA/CD\ntT99A3gtZUYF/OvpG8BrKTMS4A9N3wBeS5kRAH9k+gbwWsqcOvhj374BvJYypws+7ds3gNdS5jTB\np3/7BvBaypwe+H7fvgG8ljKnBb7/t28Ar6XMqYDP/PYN4LWUOQXwA759A3gtZbSDH3jyJMBrKaMX\nPOXkSYDXUkYneNrJkwCvpUxu8M82w2zwfnD42aknT/JChO/fOQAf+vTTtXKB9++0Av8rzvkCnFtx\nq54hfqgVfNAnSwb49mJ75WEGeM90rFr4Gjz95ElWiAuGKUp/2H0NrlmLSxKNOcDfMoUQhjDKjLvq\nkcM04hdvOaZjuvrAh3ZcMTVLBvh4fW83+4NvWDEBt7YHnnXyJCfE34o/RJEtktfQsOOSDrWfCx+8\nb5ottypEYN22yCs9NUxLdSCqqm5HtZvawKud7KQ50TFlfsZvrcZr/KNHjzpputhUY/kvL+OxPSfK\nC086VG2TH9nplBbi4Ym4r36222oUROfLl4wyieauLnfKbTF364HdnBt2mIVb6rW/9SBOU+YsgBzL\nbE87Lx9U1N+pWb7PAL+50o7H+fl5mabq42T8vx15r2pU1yVDrxiPLX2mRnE/KfZcjRWic2eH85qS\n53d+LSv3yleXn1cfV4cdZvlePJSTUXAWQI5ltqe/7ty7pv7+Ji3Ldn/wD9cy9+pdtckNxPNbpGvf\ncm+27HI8fCq+Uz879WivQRFF/E19/aeVyLxjzLbqrlsfdhivFg8XVLejxiVtm/qk71xqloxN/dri\nyspmf/Ch5bYaZuMe7dq33CECMbvuGFeS1xCabsszqe+yHDt3lbOVtw0hLpTetchn/ZPD2E6rMVOq\nt+rmt9rAe2aj5aauKSfpLVuPoxR/E+OwWppd3ZvOha5Nf5vlmc59MGtWLlQq5+pFhKnXHD+IE4Qa\np3O+Y6dnOfXLpInCARy+5dRPvRpGCIDnWwCebpFawgA8XQDPtwA83SK1hAF4ugCebwF4ukVqCQPw\ndAE83wLwdIvUEgbg6QJ4vgXg6RapJQzA0wXwfAvA0y1SSxiApwvg+RaAp1ukljAATxfA8y0AT7dI\nLWEAni6A51sAnm6RWsIAPF0Az7cAPN0itYQBeLoAnm8BeLpFagkD8HQBPN8C8HSL1BIG4OkCeL4F\n4OkWqSUMwNMF8HwLwNMtUksYgKcL4PkWgKdbpJYwAE8XwPMtAE+3SC1hAJ4ugOdbAJ5ukVrCADxd\nAM+3ADzdIrWEAXi6AJ5vAXi6RWoJA/B0ATzfAvB0i9QSBuDpAni+BeDpFqklDMDTBfB8C8DTLVJL\nmLEHH4bq2cO9n4oLEWoBH+oEz1tg3DK7z67A9y1zAvAty7Ls9ZIQYkbYWY0lT7isfNM2a//rmrZo\nsHws8F5cxNUG/i5zgfHKJGBUg7BbSah05Qcfqo4BvzTMB9G7xtkocmoFLauW6qXjlp0oCugdgpQ4\n4BtWoNocaQK/UQkYTZW4ZULhqzdyGH1WyejdlB98XT2jW7rWiZzZmZiNSX8Hs5aV6yVhVH+CBuPN\nxQOf9E0LhSbwl+6qkbHAWGW8BLXTiKoPVCgz/VEk8N00Ld1Vw8zKdrd67fJGt1v9tkvVDvmRsa7G\nz919aiZjlWPc3qY/dveJq9uvOAVyhNmr9pdkfFpMmaW1vbGaLIBK+qN+yL3GJ6ufJy52oneFGURB\nn3fWiVeSZMsSGl/t/0gWZ43fbXNkalrjl25GqtlOQWUathqtVnTtVhIq/VEn2Lmz3aBVK5Xm7ljG\nubBheQUtq9CK61gXrUZY520cOeADsx42zIYm8E/LcTWLt6vKKGM7MZh4n+hxOQmV/qCTTOfqtu2F\n75nizG9du1ZQO9Yo6UBUa/T8ms3cD2bt1YeOXWtp26t/ylxgzDIKTPyXbCeh0oUDOIcsUkuYsT+A\nM/DZhxcC4PkWgKdbpJYwAE8XwPMtAE+3SC1hAJ4ugOdbAJ5ukVrCADxdAM+3ADzdIrWEAXi6AJ5v\nAXi6RWoJA/B0ATzfAvB0i9QSBuDpAni+BeDpFqklDMDTBfB8C8DTLVJLGICnC+D5FoCnW6SWMABP\nF8DzLQBPt0gtYQCeLoDnWwCebpFawgA8XQDPtwA83SK1hAF4ugCebwF4ukVqCQPwdAE83wLwdIvU\nEgbg6QJ4vgXg6RapJQzA0wXwfAvA0y1SSxiApwvg+RaAp1ukljAATxfA8y0AT7dILWEAni6A51sA\nnm6RWsIAPF0Az7cAPN0itYQBeLoAnm8BeLpFagkD8HQBPN8C8HSL1BIG4OkCeL4lN/hnm8/6gQ/d\nkmEYwg6iD4SIf4j/S7xTQIi6ec5ICol6FNimMIRxmeplgG/FT+3sd6FqnYvLnGtRvVwidWGX7zI9\nOcrEkk68uKp97vWeAf7ZyubKsz7gaxes2bfNM5fN5ao4F8wI4050ZZbcP4IcwrPCklGeN43zL+xl\ns1EV79lfinmimQ4+6avk2rvgA9OqR+5ZclMEJpG6HUbfVnhNCnKUUfq5+DQKShfSf5kBvvmwt+mn\ngw8sa9WJvHdNt/SxiKKfiPK5KDKHv6zMoGGsBiK4YvRC4UZGENn+l9SWHnTwbvKWtXYbFbjXVGsP\n60p92GF2JeKF1N2weSZ+maTUh/EQGOlUssC3e+1mr1epVI796v47lc+v94Lzlc9Lvznb650XN+Lx\ncuVFb8gSvc+NF/GrU2Pp8/g/e9fvq3HIip/0YOxd/831g7EAVd4YC5YIkvF+6i9lMmaB76Ws8aF5\n4SM7XjesmvkPhlrjS2/F/0ZuR0R+91qtlnHNL/lnjF7LdCKxHtmt96ll6Gt83VGjGSRrfP2KqQL+\nxBt2mF2pLnrdDVYftTxllEpqr+tTwV7j/U31/1TwkfsjcbZSmv2J7VdmzL+Pd+9uBmdL9aGH8E2/\nYoi/M4y371u+5f5CVGrvG+tEM+Mz3qqHQc3d/YwPzXO1Ly9UrKGH2VXD8qO7grzrmLeM0r8Zb0fr\nxpX0X2aA31pprmz1AR95lZIQlhtG9+x4H3j2xyUxS99hoYfwa7aaMMzY533VlKgshEnlztmrDx3b\nrkd7e/WhY4oyvfERl4hv25e+ZnpylIkl75ZE6f0+v8ycx7f3/s58dszj+ZYRn8cfKPezDy8EwPMt\nAE+3SC1hAJ4ugOdbAJ5ukVrCADxdAM+3ADzdIrWEAXi6AJ5vAXi6RWoJA/B0ATzfAvB0i9QSBuDp\nAni+BeDpFqklDMDTBfB8C8DTLVJLGICnC+D5FoCnW6SWMGMDvpOh3z3usNXlW27xLVtbbMvj32kJ\n88c/aimTieYHCvgszT/K72VIy2mpvUfzWsrcvq2lDAENwO8K4On6/Z+1hNCzqP78ey1lHulZWQho\nTgAeGmflB/9sM9TxAsP96zeL1kMNNZrNpo40WwQ0ucG3F9srGpZVe7W9qIV8c6X4Gs9WtIDfWnm4\n0h70oNzg4zfV3mVWhSou0xwYYggKdYBvN9tbGrL4m72tge+vE3zGb63q2Dr21hY1LKyt1Z4G8P5i\nc1HD52NzZWV14DLLD35z8OZkSDn84mus+W0dSOI9CQ1byeYmoUxu8A/Xik+gQjzTsqw2m82fbRZe\n5aGeMP7DIsGvLa6sFL+oeg/jraOe3Xodn/ErzRUdO3eUZTby8/gtTR8oExWGUGbkwUPFCOCnVAA/\npQL4KRXAT6kAXunz1d53l0/7RegVwCeqBOfvn/Zr0CuATxSI66f9EjQL4BN9B/DTqfP3K8Fpvwa9\nAnil69d7gZ6TOkdGAD+lAvgpFcBPqQB+SgXwUyqAn1IB/JQK4KdU/w+ZhbrQuOPkEAAAACV0RVh0\nZGF0ZTpjcmVhdGUAMjAxNi0wNy0xN1QxMTowNzo0MCswMDowMB7mkEIAAAAldEVYdGRhdGU6bW9k\naWZ5ADIwMTYtMDctMTdUMTE6MDc6NDArMDA6MDBvuyj+AAAAIHRFWHRwZGY6SGlSZXNCb3VuZGlu\nZ0JveAA1MDR4NTA0KzArMKV3vKMAAAAUdEVYdHBkZjpWZXJzaW9uAFBERi0xLjQgHEc6eAAAAEp0\nRVh0c2lnbmF0dXJlAGRjODgxNWRmNGQ4ZTA3YzMxNzEzOWVkNmJiNDk3ODM1ODgxNTc4OWExYTZk\nYzA4NWU4NTc5OWJhOGRlMTE2ZTOELQNJAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "graph = preGraph(\"images/fig7.3-R.pdf\")\n", "r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=position_jitter(width=.2, height=0), shape=21, size=2.5) + geom_line(aes(x=x, y=8*fit$fitted.values))')\n", "r('plot(p)')\n", "postGraph(graph)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t久保本に従い、$x_i = 4$のデータを抽出し、そのヒストグラムとGLMから推定された\n", "\t\t二項分布を表示してみましょう。\n", "\t

\n", "\t

\n", "\t\t以下の様にPandasの抽出機能を使ってx=4のデータを取り出し、頻度を求めます。\n", "\t\tこれにロジスティック_logisticと二項分布_pの２つの関数を定義します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "y\n", "0 3\n", "1 1\n", "2 4\n", "3 2\n", "4 1\n", "5 1\n", "6 2\n", "7 3\n", "8 3\n", "dtype: int64\n" ] } ], "source": [ "tbl4 = d[d.x==4].groupby('y').size()\n", "tbl4_lst = zip(tbl4.index, tbl4.tolist())\n", "print tbl4" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# ロジスティック関数の定義\n", "def _logistic(z):\n", " return 1/(1 + exp(-z))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 二項分布を定義\n", "def _p(q, y, N):\n", " return binomial(N, y)*q^y*(1-q)^(N-y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\tGLMで求まったモデルを使って、x=4でのqの値をロジスティック関数から求めます。\n", "\t

\n", "\t

\n", "\t\tx=4での種子数分布（青点）と推定された二項分布（_p関数の値）にデータ数（tbl4.sum()=20）を掛けた値が、\n", "\t\t推定された種子数になります。\n", "\t

\n", "\t

\n", "\t\tこのグラフからは、推定された種子数は、観測された種子数を説明しているように思えません。\n", "\t\t久保本の説明では、$x_i = 4$の生存種子数の平均と分散を求めてみると、平均4.05に対して、\n", "\t\t分散は8.37と大きく、二項分布の分散$N p (1-p)$から期待される値$8\\times0.5\\times(1-0.5)=2$と\n", "\t\t比べて４倍ほど大きな値となっています。久保本ではこの状態を「ばらつきが大きすぎる」過分散と呼んでいます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.472527695655406\n" ] } ], "source": [ "# x=4での生存確率qを求める\n", "q = _logistic(-2.15+0.51*4)\n", "print q" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAECCAYAAAARlssoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xd8U+Uex/HvOdlJkw5ogSJbQOWCAwcgKktRFBBklg2V\njWwKCMgupewhyBJliAiooF4ciMJlKIIge2+EAh1p9nruH4VCKaMjyXOS/N6vl6/ehpJ87iHNL09y\nzonAGGMghBBCAIi8AwghhEgHDQVCCCHZaCgQQgjJRkOBEEJINhoKhBBCstFQIIQQko2GAiGEkGw0\nFAghhGTz21BgjMFoNIKOlSOEEOmS++JKr1/PzHWZ2ZyJcuVK4uzZy9Dp9L64Wa8SRQFRUTqkpprh\n8Uh7kFGrb1Crb1CrbzyqNTo6b4+7flspCIKQ46vUiaIAQRAgitLvpVbfoFbfoFbf8FYrvadACCEk\nGw0FQu5DSEuFekYysH8/7xRC/IqGAiF3EUyZ0E5PQtTz1aCZOA5o2hRCehrvLEL8hoYCIQBgtUKz\nYB6iXqgG7axpsMV1gHHLNiAjA9qBHwC01xwJET7Z+4iQgOF0Qr16BbQzpkJMuQZbXEdYBg+DJ7Yk\n5HIRWLwYylatoK5TH7Z2HXnXEuJztFIgocnthuqrNYiqVR1hwwbCWfNlpO74C6bps+GJLXnn51q2\nhL1DZ4R9OAyykyf49RLiJzQUSGhhDMrvNyGybi0Y+nSH68kqSNu6E5kLl8JTvsJ9/4plchLcJR+D\nvkdXwG73czAh/kVDgYQGxqDYugURb9ZFeJd28MQUR9p/t8D4+RdwP1Xl4X9Xp4Nx4TLITxyDbuJH\n/uklhBMaCiToyf/YjfBmbyOidTNAJkf6hu+Qse5buKq/kOfrcFetBvPocdB+8jGUW37yYS0hfNFQ\nIEFLfvAADHEtENn4DYgZGchY+SXSv/8ZztqvFuj6rN17w17/dej79YRw7ZqXawmRBhoKJOjITp2E\n/v3OiKz/CmRnTsO46FOkbdkOxxtvAYU5zYogIHPOQkAQYejXA/B4vBdNiETQUCBBQ7x4AWH9eyOy\n9gtQ/PUnMmfNR9r/9sD+7nuA6J27OouOhnHeJ1D+9is0C+d75ToJkZJ8/6aIogiNRgOtVpv9tX//\n/r5oIyRPhGvXEDZiCKJqPAvVzz/CNHEKUnf/DVtcB0Du/UNxnHXrw9L7A+gmjYX8wN9ev35CeMr3\nb4wgCDhx4gRKlSrlix5C8kxIS4V23mxoliwEU6pgHjYS1viegE7n89s2jxwDxY7t0PfoirRftgNh\nYT6/TUL8Id8rBcYYfVAO4UowZUI7Y2rW+YmWLoKlZx+k/vUPrP0H+2UgAACUSmR+shSyq1ehHznU\nP7dJiB8U6IXWhIQElClTBlFRUejRowfMZrO3uwjJzWaDZuGt8xPNmApbXHvc3PMPLCPGgIVH+D3H\nXf5xZE6ZBvWaVVB9vc7vt0+IL+R7KNSsWRNvvPEGTp06hV27dmH37t3o06ePL9oIyeJ0Qv35p4h6\n6Rnoxo2GvVFjpP6xH+YJU8Cio7mm2VvHwda8BcKGDIB4/hzXFkK8QWCFfC1o8+bNaNKkCcxmMxQK\nBQDg5k1Trk//MZtNKFWqOC5evAqdTvqvv8pkIgwGDYxGK9xuae96GLStbjeUG9ZBPWUSxHNn4Xiv\nJWwJI+Gp8Li0Wo0ZMLxaCyymGDK//xG49XvgT0F7H+AsmFojI/P20mqhh8KxY8dQpUoVXLhwASVL\nZp1IjDGW62M3jUYjwsPDkZGRAYPBUJibJMGOMeDbb4HRo4FDh4AmTYAJE4Bq1XiXPdju3UDt2sDw\n4cDEibxrCCmwfO19tH//fqxcuRLTpk3LvuzIkSNQqVSIjY3Nviw11XyflYIVAG5NMVlhmv0imJ4h\nSMlDWxmD/Let0EwaD/m+v+B8rQ6sP22F+/lbp6NI8+97V/narpWrQj38Q6gnT4DppZfhKuBR0wUV\nNPcBiQmm1ryuFPI1FGJiYrBo0SLExMRgwIABOHfuHMaMGYMePXrkWBl4PAweT84FyO1It9sDl0va\nG/dugdQbyK3yP/+ALnE8lDu2w1n9BaSv3wTnK69l/SHn/0953a6mvgMh+20rtD3ikbZ1B1hUET/U\n5RTI9wEpC6XWfL3RHBsbix9++AHffvstihYtitq1a6NRo0ZISkoqcAAJbbKD/8DQriUi33kdYlpa\n1vmJfvjlzkAIJDIZMj9eDMFmhX5AX/q0NhKQ8n3wWu3atbFjxw5ftJAQIp48AX3iRKi/2QBX+Qow\nfrIM9qbNvXY6Cl48JWKROetjhHdqC/XypbB1ieedREi+BPZvIAk4wqWLQLduMNR8Hoo9fyJz5rys\n8xM1axHwA+E2x1tvw9q5G8I+GgnZ0SO8cwjJF/qMZuI3wo0bMLxSE1ApYZ00BeZ2nQG1mneWT5jG\nTYZi904YenZF2uatgEbDO4mQPAmOp2YkIGjnzoTgdgOHDsHeo3fQDgQAgEYD48JlkJ05jbBxo3jX\nEJJnNBSIX4hX/4Xm08Ww9e4LcD4K2V/cT1WBadxkaJYthvK/3/POISRPaCgQv9DOmAqm0cDWux/v\nFL+ydYmH/c23oR/QG+K/V3jnEPJINBSIz4nnz0G98jNY+g4EQu1odkFA5sx5YCo19H26A2437yJC\nHoqGAvE53bQp8EQVgbVbd94pXLAiRZD58WIodmyHZt4s3jmEPBQNBeJTshPHofpqDSwDhwBaLe8c\nbpy1X4X1g0HQTZkI+V9/8s4h5IFoKBCf0k6dDE9sSdjad+adwp152Ei4nnkWhp7xEIwZvHMIuS8a\nCsRn5AcPQL3xa1iGDAdUKt45/CkUMC5YCiH1JsKGDaLTYBBJoqFAfEY7ZSJcFR6HrVVb3imS4Slb\nDqbkmVBv+AqqtV/wziEkFxoKxCfke/6A6ucfYRk2EpDTgfN3s7/XCrZWbRE2fAjEM6d55xCSAw0F\n4hO6xAlwPfWfrJPckVxMU6bBExMDQ8+ugMPBO4eQbDQUiNcptv0G5f+2wTx8VNCc5M7bWJgemZ8s\ng/zQQegSJ/DOISQb/cYS72IMusTxcD5XHY6Gb/GukTTXM8/BPPIjaOfPhuK3X3nnEAKAhgLxMuVP\nm6HY+xfMI8YA93xON8nN2rsfHK/Vhb5vDwg3bvDOIYSGAvEijwe6xAlwvPwKnK/W4V0TGEQRmfM+\ngeB2Qd+/F+2mSrijoUC8RrXxa8iPHKJVQj55ihVH5pwFUP38IzRLFvLOISGOhgLxDpcL2qRJsDd4\nA64XX+JdE3Acr78JS/de0I0bDdnBf3jnkBBGQ4F4heqrNZCfPgXLiNG8UwKWefR4uCtWztpN1Wzm\nnUNCFA0FUnh2O3TTpsDe+F24qj7NuyZwqVQwfrIMsksXETZmBO8aEqJoKJBCU6/8DOLlSzAnfMg7\nJeC5K1WGaWISNCuWQ7npG945JATRUCCFY7FAOzMZ9hat4a5UmXdNULC17wT7O02hH/QBxEsXeeeQ\nEENDgRSKZtliiKk3YR4ynHdK8BAEZM6YAxYWBkOveMDl4l1EQggNBVJgQqYR2rkzYGvXCZ6y5Xjn\nBBUWEQnjgqWQ7/kD2pnJvHNICKGhQApMs3A+BKsVlkFDeacEJVeNmrAMGgbt9CTId+/inUNCBA0F\nUiBC6k1oFsyDtXM8PCVieecELcugYXA9/yIMveMhpKfxziEhgIYCKRDtvNkAY7B8MIh3SnCTy2Fc\nsASC0YiwIQPoNBjE5wo8FAYOHAiRToscksRrV6FZ+gmsPXqBFS3KOyfoeUqVRuaMOVBv/BrqVZ/z\nziFBrkCP6vv378eKFSsg0PltQpJ21jQwpQrWXv14p4QMR5NmsLbvhLBRCZCdPME7hwSxfA8Fxhh6\n9eqFwYMH+6KHSJx48QLUn38KS9/+YOERvHNCimnCFLhLPgZ9j66A3c47hwSpfA+FhQsXQqPRIC4u\nzhc9JJ/27xcxbZoCGzb45/a005PAwiNgje/pnxskd+h0MC5cBvmJY9BN/Ih3DQlS+fpE9WvXrmHs\n2LHYtm2br3pIPuzdK6JpUy0cjqyX8caNk6NXL9993q/s9Emov1wN87hJgE7ns9shD+auWg3m0eMQ\nNnoEnK/VhaNBQ95JJMjkaygMHjwY3bp1Q+XKlXH+/PkH/pwoChDFnO83yGRi9le5XPpvUN/dK1U/\n/aTIHggAsHGjAv36+e7oV11yIlix4nB2jS/wv2EgbNfbpNrq7N0Xzm2/Qf9BLxi37QYrXlyyrfdD\nrb7hrdY8D4UtW7Zg586dWLx4MYCs9xYeJCpKl+tNaJnMDQAwGDQwGALnWabBoOGd8EBVquT8vlIl\nEZGRPtq2//wDbFgHLFqEyBJFCn11Ut6u95Jk68rPgWrVENG/F7B5M3BrT0BJtj4AtfpGYVvzPBRW\nrVqFlJQUlC5dGgDg8XjAGENMTAzmzZuHVq1aZf9saqo510rBbLYCAIxGK9xuWaGi/UEmE2EwaG71\nenjn3FezZsDBgwps3ixHlSoiJk+2Ii3NN6264SMhK1cexqYtgbSCn+s/ELbrbZJuVeggn78I+hZN\nYZk4Ba4BA6Xbeg9Jb9d7BFNrXp8wCuxhT/nvkpGRAfNdH/xx8eJF1KxZE5cvX0ZkZCTUanX2n12/\nnpnr71ssJpQtG4tz565Aqw3LUxxPcnnWs+60NDNcLmnfGXzdKt/3FyLfrAfjx4thb9G6cNdF29Wr\ndGNHQbN4ATJ//BWGurUl3XpbIGzX24KpNTpan7fryesNhoeHIzw8PPt7p9MJQRBQokSJvF4FCVC6\nxAlwPfEk7M1a8E4h9zCPHAPFju3QxXcG9v8NOkkBKawC34PKlCkDt9vtzRYiQYqd/4Py960wJ4wC\nZNJ/2S/kKJXI/GQpxGtXgQ/pQ45I4dHTCvJgjEE3eTycTz8LR6N3eNeQB3CXfxy2QUOBhQshXrzA\nO4cEOBoK5IGUv/4MxZ+7YR4xCqBTmkiarXsvIDIS6uQk3ikkwNFQIPfHGLSJE+F8qSacdRvwriGP\notMBI0dC+cVKyM6c4l1DAhgNBXJfyu82QvHPfphHjqFVQqDo0QOsWHFop07mXUICGA0FkpvbDV3S\nRDjq1IOz5su8a0heqdWwDkmA6uv1kB05zLuGBCgaCiQX1fq1kJ84DvOI0bxTSD452nWAp3QZ6JIm\n8U4hAYqGAsnJ6YQuORH2t96B69nqvGtIfikUMA8dAdV/v4N8/z7eNSQA0VAgOahXr4B44TzMw0fx\nTiEFZH+vFVyVKkOXOIF3CglANBTIHTYbtDOmwt6sBdxPPsW7hhSUTAZzwodQbt0Cxe6dvGtIgKGh\nQLJpli+BmHIN5mEjeaeQQnK83QTOqk9DO3k8kLfTmxECgIYCuc1kgnbODNjatoenfAXeNaSwRBGW\nEaOg3L0Tit9+5V1DAggNBQIA0C5eAMFohGXQMN4pxEsc9d+A8/kXoUuk1QLJOxoKBEJ6GjTz58Da\nqSs8j5XinUO8RRCyzqK6/28oN//Au4YECBoKBJqP50JwOWHpP4R3CvEyZ+1X4XjlNeimTAQ80v48\nACINNBRCnHD9OrSLFsAa3xMsJoZ3DvEB84jRkB89DNW3G3inkABAQyHEaedMB5PJYOnzAe8U4iOu\n51+E/Y03s86J5HLxziESR0MhhIlXLkOzfCmsvfuBRUbxziE+ZE4YBfnpU1Cv/YJ3CpE4GgohTDt9\nKlhYGKw9evNOIT7mrloNtibNoJ02BbDbeecQCaOhEKLEs2eg/mIFLP0GgYXl7QO9SWCzDBsJ8cpl\nqFd+xjuFSBgNhRClS06Ep2g0rF3ieacQP3FXqgx7i9bQzkwGLBbeOUSiaCiEINmxo1CtXwvLwKGA\nRsM7h/iRechwiKk3ofl0Ce8UIlE0FEKQLmkSPKVKw9auI+8U4meesuVgi+sI7dwZEDKNvHOIBNFQ\nCDHyA39D9f1GmIcMB5RK3jmEA8ugoRDMZmgWLeCdQiSIhkKI0SVOgKtiJdhbtuGdQjjxxJaEtXN8\n1pHsaam8c4jE0FAIIfLdu6D89ReYEz4EZDLeOYQjyweDILjd0M6fwzuFSAwNhVDBGHSJ4+H8TzU4\n3mnKu4ZwxqKjYeneC5olCyGkpPDOIRJCQyFEKH77FcpdO2AZMQoQ6Z+dIOtIdrkC2jnTeacQCaFH\nh1Bwe5Xw/ItwNGjIu4ZIBIuIhLV3P2iWL4V4+RLvHCIR+R4KBw4cQIMGDRAREYESJUqgTZs2uHbt\nmi/aiJco//s9FPv/hnnkGEAQeOcQCbF27wWm10M7I5l3CpGIfA0Fh8OBhg0bol69erh+/ToOHTqE\na9euoXdvOneOZLnd0CVNhOOVOnDWfpV3DZEYFqaHpd8gqL9YAfHsGd45RALyNRQsFgsmT56M4cOH\nQ6FQoEiRImjevDkOHTrkqz5SSKpv1kN+9AjMI0fzTiESZe0SD0+RotBNm8I7hUhAvoZCREQEunbt\nCvHWG5XHjx/H8uXL0abNo/d5X7dODgDYtIl2hfQbpxPaqZNhb/gWXNVf4F1DpEqjgWXgUKjWfQnZ\n8WO8ayQjJUXA2LEKDB4MXLok7ZddMzKASZMUGDAAOHmycK0FeqP5woULUKlUqFKlCl566SWMHTv2\noT+/dKkCQ4aoAAD9+qnx1VfygtwsySf1l6shP3sG5oRRvFOIxNnad4KnVGnokibxTpEEpxNo3lyD\nOXOUmDEDePttNcxm3lUP1qaNFtOnKzF7NtCokQY3bhR8MBTo0bl06dKw2+04ffo0unfvjvbt22PV\nqlXZfy6KAkTxTtT27XIAjhzft20r7c+LlcnEHF+l7L6tdjt005PgaPYehGeeLtg/tA8E/HaVqEK3\nytWwDRsBXd+eUB3+B+6nn/FiXU6BsF0vXxZw4sSdVzUuXhRx/rwc1apJ73HLZAL27r3TevOmgCNH\n5KhXz12g6xMYY6wwQbt370atWrVw/fp1FClSBADAGINw114uY8cC48YZAYQDyMCMGQYMHFiYWyWP\nNGcOMHAgcOQIULky7xoSCFwu4D//ASpUAL7/nncNV1YrUL48cPVq1vcREcDp00CURD+gsGJF4NSp\nrP+t0QBHjwJlyhTsuvL1BHLr1q3o1asXjh2787qjIAgQBAHKu06ulppqzrFS6N0buHHDifnzgcGD\nnejY0Yy0tIIF+4tMJsJg0MBotMLtlt6zg7vlajWbET5xIpxt2sES8xiQJp11b0BvVwnzVqti6AiE\nxXeGcfMWuF+q4cXCOwJlu65bJyAxUQVAhsGDbRAEt2Qft9auFTBunAo2mwy9e9thMLhytUZG6vJ0\nXflaKRiNRjzxxBPo0KEDxo4dC5PJhE6dOsFqtWLr1q3ZP3f9emauv2uxmFC2bCzOnbsCrTYsrzfJ\njVwuIjJSh7Q0M1wu6d5xgdytmjkzoEuahNRd++ApXcCnCz4SyNtVyrzW6vEgsl5teKKikLHhO+8F\n3iUkt6sfPKo1Ojpvn7CYrxf1DAYDfv75Z/z555+Ijo5G1apVERkZidWrV+fnaogPCRnp0M6blfXG\nocQGAgkAogjz8FFQ/m8bFNt+411DOMj3+49VqlTJsSog0qJZMA+CzZb1qWqEFICj4VtwPlcdusQJ\nSH/lNToKPsRI9+1/km/CjevQfPIxrF27w1O8BO8cEqgEAebho6HYuwfKnzfzriF+RkMhiKhnzwQE\nAZZ+tGsXKRzna3XhqFUbusSJgEfar6UT76KhECyuXIFq6SJYe/QGu7VrMCEFdmu1ID98EMrvvuVd\nQ/yIhkKwmDgRTK2BtVdf3iUkSLhq1ISjXoOso5zdBTsQigQeGgpBQDx/Dli8GLb+A8EM4bxzSBAx\njxgN+ckTUK37kncK8RMaCkFAPTURKFIE9vgevFNIkHE9/SzsjRpDlzwFcDge/RdIwKOhEOBkhw5C\n+eUXwIcfArq8HbFISH6Yh4+CePE81KtX8E4hfkBDIZB5PNAnDIKnYiWgZ0/eNSRIuZ94EvbmLaGd\nMTXrpEAkqNFQCGCqtV9AsecPWJJnAAoF7xwSxMxDR0C8ngLNZ0t5pxAfo6EQoIT0NISNHw1b85Zw\n0cdsEh/zlK8AW9v20M6ZkXWuZhK0aCgEKN3k8YDdAfM4+lAU4h+WQcMgGI3QLlnIO4X4EA2FACQ/\n8DfUny2DJWEkPMWK884hIcLzWClYO3aBZv4cCBnpvHOIj9BQCDQeD8ISBsH9ZBVYu3bnXUNCjKX/\nEAgOOzQL5vJOIT5CQyHAqFd9DsW+vcicMh2QS+VDNkmoYMWKwdqtBzSfLIBw4wbvHOIDNBQCiHDz\nJnQTP4KtdRxcNWryziEhytK3PyCKWW86k6BDQyGA6CaPA9wemMZM4J1CQhiLKgJrzz7QLF8C8d8r\nvHOIl9FQCBDyvXugXvkZzCNGg0VH884hIc7asw+YRgPtzGTeKcTLaCgEArcbYQmD4fpPNdg6d+Nd\nQwiY3gBL34FQr/ws64SMJGjQUAgA6s+WQfHPfpiSpgMyGe8cQgAA1q7vg0VGQTc9iXcK8SIaChIn\nXL8OXeIEWNt3guv5F3nnEHKHTgfzwCFQrf0CspMneNcQL6GhIHFhE8YAogDzh2N5pxCSi61DF3hK\nxEKbPJl3CvESGgoSJv9jN9RrVsH84Vj6iE0iTSoVLIMToP5mA2SHDvKuIV5AQ0GqXC7oEwbB+Vx1\n2Np34l1DyAPZWsfBVa48dEkTeacQL6ChIFGaZYsgO3oYpinTAZH+mYiEKRSwDBsJ1Y//hXzvHt41\npJDo0UaCxGtXoZ0yCbZOXeF65jneOYQ8kv3d9+B64knoEmm1EOhoKEiQbuwoQKWEecRo3imE5I1M\nBnPCKCi3bYVix3beNaQQaChIjGLHdqjXr4VpzASwyCjeOYTkmaPRO3A+/Sx0iRMAxnjnkAKioSAl\nTifChg+G84WXYG8dx7uGkPwRBJhHjILiz91Q/voz7xpSQPkeChcuXEDz5s1RtGhRlChRAl26dIHR\naPRFW8jRLFoA2ckTWafFpjeXSQBy1m0A50s1oU2cSKuFAJXvR57GjRsjKioKFy9exN69e3H48GEM\nGTLEF20hRbxyGbrkRFi7dYe7ajXeOYQUjCDAPGI0FP/sh/L7TbxrSAHkayhkZGTghRdeQGJiIjQa\nDWJjY9GpUyds27bNV30hQ/fRh2A6HSwJH/JOIaRQnLVqw/Fa3azjFtxu3jkkn/I1FMLDw7FkyRJE\n33Xq5gsXLqBkyZJeDwslit+3Qv3tBpjGTgQzhPPO8ZlNm+To2FGFQYMAk4l3DfEl84jRkB8/BtXX\n63JcfuiQiPh4FTp2BM6eFTjVkYcp1Oc5/vXXX5g3bx6+++67HJeLogBRzPkPLpOJ2V/lcum/Xn53\nr0/Z7dCPGAJnrZfhbtMWciH/vyh+ay2EXbtExMerwVjW/7/z51VYssTOuerhAmG73ia51hdfhOOt\nt6FLToT7vRaAQoH0dKBFCy1SU7PuA7//rsaePVYoFJxbH0Jy2/UhvNVa4KGwY8cONGnSBFOnTkXd\nunVz/FlUlA6CcO9QyFpGGgwaGAy6gt6s3xkMGt/eQOIc4OwZyL7egMiosEJdlc9bC+Ho0ZzvO+7b\nJ0dkZGB8xrSUt+u9JNU6ZTLw9NOI3LgOiI/HqVNAauqdP75wQYTdrkNMDL/EvJLUdn2EwrYW6Ldy\n06ZN6NChA+bPn4927drl+vPUVHOulYLZbAUAGI1WuN3S/0wAmUyEwaC51evxyW0Ily4ifMIE2Hv2\ngbVkOSDNXKDr8UdrYVWtKkIU1fB4su4XNWq4kJYm/ZWC1LfrbZJsLVUBumbvQT52HDLeaY7oaBVi\nYjRIScl6JluhggcqlRVpaXwzH0aS2/UBHtUaGZm3J+P5Hgo7d+5E586dsX79etSvX/++P+PxMHg8\nOXdHux3pdnvgckl7497Nl72G4cPgCY+AaXACmBduQ8rb9tlnPVi1yopvvlGgUiUFevWyS7b1XlLe\nrveSWqtp6EhE1n4BimVLoHu/F77+2ooFC5TQ6xXo1csGQfDA5eJd+WhS264PU9jWfA0Ft9uN999/\nH0lJSQ8cCCRvlFt+guqHTTAu+hQsTM87xy/q13ejYUOGyEgF0tIQEA8GpHDcj1eEvVVbaGdNhzWu\nIypW1GHOHMet+wCj+4AE5esdiV27duHYsWP44IMPoNFooNVqs79evHjRV43Bx2ZD2IihcLzyGuxN\nm/OuIcSnzEOGQ0hPg2bpIt4pJA/ytVKoXbs23LTfcaFp582CePkSMlZ9BRRgbyNCAomndBnY2neC\ndt5M2Dp3BaIieSeRh5D+flZBRjx3Fto5M2Dt2RfuipV45xDiF5aBQyHYbNAsnM87hTwCDQU/CxuV\nAE+RojAPGsY7hRC/8RQvAWuX96FZOB9C6k3eOeQhaCj4kfLH/0L102aYJkwBdIFzrAYh3mDpNxBg\nDOo5s3inkIegoeAvFgvCPhwGR936cLzdmHcNIX7HihaFtUcvqBYvBP79l3cOeQAaCn6inTMd4tV/\nYUpMpjeXSciy9uoHplQBCQl0am2JoqHgB7Izp6CdNxuWvv3hLv847xxCuGHhEbAmTgVWrIBq7mze\nOeQ+AuPkM4GMMYSNGApP8RKwfDCYdw0h3DnaxEF3+Ty0Y0fBWao0HI3f5Z1E7kJDwceU32+CcusW\nZHy+BtBqeecQIg3jx8Nx9DgMfbojPbYkXNVf4F1EbqGXj3zJbEbY6OGwv/EmHG824l1DiHSIIszz\nFsJV9WmEd2gD8cJ53kXkFhoKPqSbmQzx5g2YJibxTiFEetRqZHz2BVhYGMLbtYSQkc67iICGgs/I\nTp6AZsFcWD4YBE/ZcrxzCJEkVrQoMlavg3jtKgxdOwJOJ++kkEdDwRcYQ9jwwfCUfAyWvgN41xAi\nae7HK8L46Soodu9A2LCBtKsqZzQUfED17QYot/+edUyCWs07hxDJc778CjJnzIVm1efQzKUjnnmi\nvY+8TDBlQjdmJOyNGsNR/w3eOYQEDHvrOJjPnkHYxI/gLlsWjibNeCeFJBoKXqZNngIxIx2mCYm8\nUwgJOJY08C1eAAATo0lEQVSEDyE7d/bOrqrPv8g7KeTQy0deJDt6BJpFH8M8aBg8pUrzziEk8AgC\nMmfNh+uZ5xDesQ3E8+d4F4UcGgrecuvNZXfZcrD27Mu7hpDApVYjY/lqsDA9wuNaQEhP410UUmgo\neIlq3ZdQ7toB05TpgErFO4eQgMaKFEHGF+sgXk+BoVtHwOHgnRQyaCh4gZCRjrCxo2Br2hzO1+ry\nziEkKLgrVITxsy+g2L2TdlX1IxoKXqCdOhmC2QzzuEm8UwgJKs6aLyNz1nxoVq+AZs4M3jkhgfY+\nKiTZwX+gWboI5tHj4YktyTuHkKBjb9kma1fVSePgKVMW9nff450U1GgoFIbHA/3wwXBXrARr9168\nawgJWpahIyA7dxb6fj3hjn0Mrhdf4p0UtOjlo0JQfbkaij1/ZL25rFDwziEkeAkCMmfOg/PZ6gjv\n1AbiubO8i4IWDYUCEtLTEDZ+NGzNW8L58iu8cwgJfioVjMtXwRMeQbuq+hANhQLSTR4POJz05jIh\nfsSiisC4+iuIN2/A0KU97arqAzQUCkC+fx/Uny2DJWEkPMWK884hJKS4yz+etavqnj+gH9KfdlX1\nMhoK+eXxICxhENxPVoG1a3feNYSEJGeNWsic/THUa1ZBO2sa75ygku+h8OOPP6J48eKIi4vzRY/k\nqVd+BsXf+5A5ZTogp523COHF/l4rmIeNhC5xAlQbvuKdEzTy9aiWnJyMZcuWoVKlSr7qkTTh5k3o\nJo2FrXUcXDVq8s4hJORZBidAdvYM9P17w12yFFwv1eCdFPDytVLQaDT4888/UaFCBV/1SJpu0ljA\n7YFpzATeKYQQIGtX1Rlz4az+QtauqmdO8y4KePkaCn379oVer/dVi6TJ//oTmpWfwTxiNFh0NO8c\nQshtKhWMn66EJzIK4e1aQkhL5V0U0OiN5nu43cDEiUq8/roaffoANlvWhWEJg+Gs9gxsnbvxTiR+\nsHatHA0bqvHuu8ClSwLvnIfavFmGRo3UaNQIOHFC2q2+wiKjkLHqK4hpqbSraiH55J1SURQgijnv\nnDKZmP1VLpfuLFq0SI45c7JOfb13LyCTKTGl1CdQHDwA409bIVdJ78jlu7et1AVC6969Ivr1U4Ox\nrPvwpUtq/PyzlXPV/Z0+LaBbNw2czqzWo0fV2LdPmq23+ew+UKkiTCu/hP7dt2EY8gEs8z8BhMIN\nyUC4v97mrVafDIWoKB0E4d6h4AYAGAwaGAw6X9ysV5y95+j564fToP18PBAfD8Prdbg05ZXBoOGd\nkGdSbr1yJeeu78eOiYiMlOZ9NiUFcDrvfH/unAiNRge1ml9TXvnkPvBWA2D5cqji4qB66glg9Giv\nXK2U76/3KmyrT4ZCaqo510rBbM569mI0WuF2y3xxs15Rt64MS5fe+Y0aZRwCjyiDMWEUWJqZY9mD\nyWQiDAbNrW3r4Z3zUIHQ+vTTAvR6DTIzs+7Db77pQlqanXPV/VWqBMTEaJCSkvXssE4dN6xWG6wS\nXiz4/D7wZhOoR46GZswYmIo/BmeLVgW+qkC4v972qNa8PrHxyVDweBg8npxHGd6OdLs9cLmku3Hf\neMODNWs82LlTjqZF/kSVj1Ygc9psOMOjAAl3A9LftneTcmtsLPD99xZs2KBAmTJKtG1rl2xreHhW\n6xdfKBETo0T79jbJtt7Ll/cBU/8hEE6fhq5vT6QXL1noXcilfH+9V2FbBcbyfoy4RqOBIAhw3lqv\nyuVyCIIAi8WS4+euX8/M9XctFhPKlo3FuXNXoNWGFTjYX+TwILLBq3DJlUj74RdAlO5rinJ51ssb\naWlmyd9xqdU3qPU+HA6Et24G+dHDSP/hF7jLP57vqwim7Rodnbc9R/P1SGe1WmGxWOB0OuF0OrO/\nD0aqxZ8ABw/CkjxD0gOBEPIASmXWrqpRRWCIawkh9SbvooBAj3b3odi6BZrEiUDPnnA/+xzvHEJI\nAbGIyKxdVTPSs3ZVtUvzvSEpoaFwFyHTiLDBHyCidTO4qlcHJk/mnUQIKSRPufLI+GwNFPv+gn5g\nXzqr6iPQULhF8duviHy1BlQb1iEzeRZMGzYBERG8swghXuB68SVkzl0I9bovoZ02hXeOpIX8aT6F\nTCN0Y0dDs+JTOF6pg8yZc+EpXQbyQh70QgiRFvu778F87ix0k8fDXbYc7C3b8E6SpJAeCorffoV+\nUD8IaWnInDoTtk5dC30EJCFEuiz9B0M8ewb6gX3heawUnDVf5p0kOSH58lHWewf9EdHqXbjLlUfa\n77uyzmlEA4GQ4CYIMCXPgvPFGjB0joPs9EneRZITckNB8ftWRL5WE+r1a5E5dSYy1m2Ep3QZ3lmE\nEH9RKmFctgKeotFZu6repF1V7xYyQyF7ddCyKdzlyiN1225aHRASorJ3Vc00IrxzHO2qepeQGAo5\nVgdJM5Dx1be0OiAkxHnKlkPG52sgP/A39P17066qtwT1UBBMmQgbMiDn6qBLPB2hTAgBALiefxHG\neZ9AveEraKfScUlAEO99pPh9K/QD+0JMTUVm0oysPYtoGBBC7uFo0gymUWcRNnFs1q6qreN4J3EV\ndENBMGVCN24MNJ8thaP2q0j/+nt4ypTlnUUIkTBrv4GQnT0D/aB+8JQqDWet2ryTuAmqp86Kbb9l\nvXfw1RpkTpmetWcRDQRCyKMIAkxTZ8JZ4+WsXVVPhe6uqkExFARTJsKGDkREiyZwly6D1N93wdb1\nfXq5iBCSdwoFjMs+hyemGMLjWoTsrqoB/6ip2P57ztXB+k20OiCEFAgLj0DGqq8gmEwI79QWsNl4\nJ/ldwA6F7NXBe41pdUAI8RpPmbLIWLEG8n/2Q9e3J+CR9ofreFtAPoLmWB0kTqPVASHEq1zVX4Bx\n/iIoN6wDOnaEbN/ekDmOIbCGgsmEsGF3rQ5+2wlbt+60OiCEeJ2j8bswz5gDbNkCQ4PXEFnzOWiT\nEyGeOc07zacC5tFUsf13RNWpCfXaL+6sDsqW451FCAlijs5dgUuXkLl+I1wvvATNgnkoUuNZRLxZ\nF5rFCyCkpPBO9DrpDwWTCWEJg7JWB4+VQupvu2h1QAjxH5kMrrr1kDl3IW4ePgXj4uXwRMdA99GH\nKPJ0ZYS3aQ7VV2sAk4l3qVdI+uA1xf+2QT+gD8Qb15GZmAxbF3ojmRDCkUYDe9PmsDdtDiH1JlQb\nv4F6/VoY+nQH02phf7MR7O+1gqNOfUCh4F1bINJ8hL29Omj+TtbqYOtO2Lr1oIFACJEMFlUEts7d\nkL7pR9z86yDMA4dCfvgQwtu1QpFqlRCWMAjyP/8IuDeoJfcoq/jfNkTVqQX1l6uRmZiMjA3fwVOu\nPO8sQgh5IE/pMrD2H4y0bX8gdcv/YGvdDsrNPyDyndcR9eLT0E6ZANnJE7wz80Q6Q+Hu1UHJkrQ6\nIIQEHkGAu2o1mMdOROq+w0jf8B0ctV+FZskiRL38PCIavArNgnkQr/7Lu/SBJPGIq9ix/c7qYPJU\nZHz9Pa0OCCGBTSaDs/arMM2ch5uHTiJj2Up4SpWGbtJYRD3zJMJbNIVqzSoImUbepTnwHQomE8KG\nD0ZEs7fvrA7ie9LqgBASXNRqON5pAuOnK3Hz0EmYps0G3C7o+/dGkSqPQx/fCcr/fg84HLxL+e19\npNixHfr+fSDeSEHm5KmwdaXdTAkhwY9FRMLWvhNs7TtBvHIZqg3roF6/FupObeGJiIC9cTPYW7SC\n86WaXB4T/X+LZjPCRgyh1QEhJOR5YkvC2rc/0rbuQOrvu2Hr2BXKrb8goulbiHq+KnQTPoLs6BG/\nNvn9kdjwZn2oV6+AaVISvXdACCG3uJ98CuZRY5H610Gkb9wMR73XoV7xKaJeq4HIOrWgmTsL4uVL\nPu/w21AQbmQdDu4pXhypW3fC+n4vWh0QQsi9RBHOGrVgmjYLNw+dQsbna+CqWAm65MmIeq4Kwt9t\nBPWK5RDS03xz8z651vtgRWMAAKa1G+ApX8FfN0sIIYFLqYTjzUbIXLwcNw+fQubsjwG5AmFDB6DI\nfyrC0LkdlJu+9ernPvj/qbpAqwNCCMkvpjfA3qYdMtZ9i9QDx2AeNRbi5UsI79YBRf5TEdp+vYFD\nhwp9O/QITQghAcZTrDisPfsi/effkbrjL1jje0C+YzuQnFzo6/b6LqmMMZjNmRAEIcflFos5x1ep\nk8lEyGRumM1WuN3S/uQlavUNavUNavWykrHI/GAAZAMHwWDQwGw03bfVaGTQ6/W5HpvvJTDm3bM1\nGY1GhIeHe/MqCSGEeEFGRgYMBsNDf8brQ4ExhvPn/73vSuGppyriyJGT0Gp13rxJn5DJRBgMGhiN\nEn6GcAu1+ga1+ga1+sajWiMjdXlaKXj95SNBEKDT6R/451qtDlptmLdv1uvkchEGgw5utwwul7Tv\nDNTqG9TqG9TqG49qNRge/Lh8N6+vFB7k9stKeVm+EEII4cNvQ4ExhszMzDwtXwghhPDht6FACCFE\n+ug4BUIIIdloKBBCCMlGQ4EQQkg2GgqEEEKy0VAghBCSjYbCffz4448oXrw44uLieKc80oULF9C8\neXMULVoUJUqUQJcuXWA0SuuDwAHgwIEDaNCgASIiIlCiRAm0adMG165d4531SAMHDoQo8c/9EEUR\nGo0GWq02+2v//v15Zz3QpEmTEBsbC71ejzfeeAPnz5/nnZTL9u3bs7fl7f/UajVkMhnvtPvav38/\n6tevj8jISMTGxqJDhw64ceNGga5L2vd2DpKTkzFgwABUqlSJd0qeNG7cGFFRUbh48SL27t2Lw4cP\nY8iQIbyzcnA4HGjYsCHq1auH69ev49ChQ7h27Rp69+7NO+2h9u/fjxUrVkj+uBpBEHDixAlYLBZY\nrVZYLBbMnj2bd9Z9zZ8/H6tXr8a2bdvw77//4qmnnsLMmTN5Z+XyyiuvZG/L2/999NFHaN26Ne+0\nXNxuN95++23UqlUL169fx+HDh5GSkoI+ffoU7AoZyWHu3LnMaDSyzp07s7Zt2/LOeaj09HTWrVs3\nlpKSkn3ZvHnzWOXKlTlW5ZaWlsaWLl3K3G539mVz5sxhlSpV4lj1cB6Ph9WoUYNNnjyZiaLIO+eh\nBEFg58+f552RJ+XLl2fffPMN74x8O3/+PCtatCi7dOkS75RcLl68yARBYMeOHcu+bOHChaxixYoF\nuj5aKdyjb9++0Ovzdo4Q3sLDw7FkyRJER0dnX3bhwgWULFmSY1VuERER6Nq1a/bLMMePH8fy5cvR\npk0bzmUPtnDhQmg0moB4CREAEhISUKZMGURFRaFHjx4wm6V3ivorV67g7NmzuHnzJqpUqYKiRYui\nZcuWBX6Zw5/GjBmD+Ph4yf1uAUDJkiXx7LPPYtGiRTCbzUhJScH69evRuHHjgl2ht6ZVsAmElcK9\n9uzZw7RaLfv11195p9zX+fPnmVKpZDKZjPXq1Yt5PB7eSfd19epVFhMTw44dO8bOnTsn+ZVCrVq1\n2LJly5jD4WDHjh1j1apVY506deKdlcsff/zBBEFgb731Frt69Sq7fPkyq1GjBmvWrBnvtIc6e/Ys\nMxgM7Nq1a7xTHujMmTOsfPnyTBRFJooiq1evHrPZbAW6LlopBIkdO3agYcOGmDp1KurWrcs7575K\nly4Nu92O48eP4/jx42jfvj3vpPsaPHgwunXrhsqVK/NOyZMdO3agS5cuUCgUqFy5MpKSkrB69Wo4\nnU7eaTmwW2fUSUhIQLFixRAbG4tx48Zh48aNcDgcnOsebP78+WjevDliYmJ4p9yXw+FA48aN0bp1\na2RkZODy5cswGAwFX+V6d14Fj0BaKWzcuJGFh4ezlStX8k7Js127djFBENiNGzd4p+Twyy+/sHLl\nyjGLxcIYy3qWKPWVwr2OHj3KRFGU3Ovf586dY4IgsL///jv7suPHjzNRFNnFixc5lj1c2bJl2YYN\nG3hnPNAPP/zAdDpdjssOHDjABEFgaWlp+b4+WikEuJ07d6Jz585Yv3492rVrxzvnvrZu3Yonnngi\nx2WCIEAQBCiVSk5V97dq1SqkpKSgdOnSiI6ORvXq1cEYQ0xMDNauXcs7L5f9+/fn2tvsyJEjUKlU\niI2N5VR1f4899hgMBgP279+ffdnZs2ehUCgk13rbgQMHcOHCBbz++uu8Ux7I7XbD4/HA47nzGQo2\nm63ge815aVgFnUBYKbhcLvbUU0+xxYsX8055qIyMDFaiRAk2bNgwZrFYWEpKCnvrrbdYnTp1eKfl\nkp6ezi5fvpz93+7du5kgCOzKlSvMarXyzsvl8uXLTK/Xs6SkJGa329nx48dZlSpV2IABA3in3deg\nQYPY448/zk6dOsWuXbvGXn75ZRYfH88764E+/fRTFh0dzTvjoW7evMmio6PZqFGjmMViYTdu3GBN\nmzZldevWLdD10VC4h1qtZhqNhsnlciaXy7O/l6Lt27czURSZRqPJ7rz99cKFC7zzcjh06BCrU6cO\n0+l0rFixYiwuLo5duXKFd9YjBcIbzdu3b2e1atVier2eRUdHs6FDhzK73c47677sdjvr27cvi4qK\nYgaDgXXt2pWZzWbeWQ+UmJjIqlatyjvjkfbt28fq1q3LoqKiWIkSJVjbtm3Zv//+W6Dros9TIIQQ\nko3eUyCEEJKNhgIhhJBsNBQIIYRko6FACCEkGw0FQggh2WgoEEIIyUZDgRBCSDYaCoQQQrLRUCCE\nEJKNhgIhhJBsNBQIIYRko6FACCEk2/8BjvTuBRzBrW8AAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "lst_plt = list_plot(tbl4_lst, zorder=2)\n", "prd_plt = list_plot([_p(q, y, 8)*tbl4.sum() for y in (0..8)], plotjoined=True, rgbcolor=\"red\")\n", "(lst_plt + prd_plt).show(figsize=4)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.05 8.36578947368\n" ] } ], "source": [ "# x = 4 の平均と分散を求める（平均に対して分散が大きいことが分かる）\n", "d4 = d[d.x==4]\n", "print d4.y.mean(), d4.y.var()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.99396217995198\n" ] } ], "source": [ "# 二項分布から期待されている分散は２なので、8.4はかなり大きい\n", "N = 8\n", "vr = N*q*(1-q) \n", "print vr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

一般化線形混合分布（GLMM）を使った分析

\n", "\t

\n", "\t\t個体差を考慮した分析手法として一般化線形混合分布（GLMM）が出てきます。\n", "\t

\n", "\t

\n", "\t\tここでは、個体$i$の個体差を$r_i$とし、以下の様なリンク関数を定義します。\n", "$$\n", "\t\tlogit(q_i) = \\beta_1 + \\beta_2 x_i + r_i\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tここで$r_i$をどのような分布にするかですが、久保本では平均０，分散$s$の正規分布と\n", "\t\tしています。（統計モデルとして便利だからと理由です）\n", "\t

\n", "\t

\n", "\t\tglmmMLを使った分析では、個体毎に異なる独立パラメータをclusterオプションで指定します。\n", "\t\tここでは、idが個体を識別する指標として使われています。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ " [1] \"glmmML\" \"jsonlite\" \"ggplot2\" \"stats\" \"graphics\" \"grDevices\" \"utils\" \"datasets\" \n", " [9] \"methods\" \"base\" " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#r(\"install.packages('glmmML')\")\n", "r('library(glmmML)')" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "\n", "Call: glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = d, cluster = id) \n", "\n", "\n", " coef se(coef) z Pr(>|z|)\n", "(Intercept) -4.190 0.8777 -4.774 1.81e-06\n", "x 1.005 0.2075 4.843 1.28e-06\n", "\n", "Scale parameter in mixing distribution: 2.408 gaussian \n", "Std. Error: 0.2202 \n", "\n", " LR p-value for H_0: sigma = 0: 2.136e-55 \n", "\n", "Residual deviance: 269.4 on 97 degrees of freedom \tAIC: 275.4 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# python版ではない、R一般線形混合モデルGLMMを使う\n", "# 個体差は、平均０の分散sの正規分布に従うと仮定して先の問題を扱う\n", "r('fit <- glmmML(cbind(y, N - y) ~ x, data = d, family = binomial, cluster = id)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

結果の可視化

\n", "\t

\n", "\t\tGLMとの違いは、分析結果から予測される二項分布を見ると明らかです。\n", "\t\tGLMのようなfitted.valuesがGLMMにはないので、分析結果の二項分布をmyfuncで定義して、\n", "\t\tstat_function使ってプロットしてみました。\n", "\t

\n", "\t

\n", "\t\tきれいに二項分布が表示されました。（これが図7.10の(A)）\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAAAAACDQ6CjAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC\nAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAJiS0dE\nAACqjSMyAAAACXBIWXMAAABIAAAASABGyWs+AAAh5UlEQVR42u2dT4jcSJaHFRlZpbSrlTVtucvj\nbdTucXvpqCILs1Namu3tXdMgaLwUXrwwechl8Cw9ixGL6TnUxXWxwODDgH1pfNGtroY9+BAw7MHH\nORgffNCpD3krTJM0dSl8MAy5Ear/VZIynjJClcp8D1uVtjLqp/e+jFBEpCKeNUSbSbPO+wLQzscQ\n/Iwagp9RQ/AzaiPAvyuwnZ/fge0XeJHdSlR+KVMGXqREzH4uEebCmO2qgE8KbLefgG0ALzKEFymh\nMihTBl6kRMz6JcJcGLN3CP54EQSP4A3KIHhNTmhTQfAI3qQMgtfkhDYVBI/gTcogeE1OaFNB8Aje\npAyC1+SENhUEj+BNyiB4TU5oU0HwCN6kDILX5IQ2FQSP4E3KIHhNTmhTQfAI3qTMpIN/s/nm/MEH\ni42mFxsHH3nO9S+ajY+6BmVuzZFWlPiO+7t2o/E5TAUG/ltKqON8eZUSi5CGE2W8pQD81srWypvz\nBu85jZXFizQ2DL5Lu4+JdfXruXlmTMYjnz7+2Gqy+PdW44FPFkEqIPCfWN6vGtZNIrBbFqFzJDz7\nngLwYW/YC88ZfGTb4qKdjmsYPI0SdmneG7wgNDYkE5N/Fsc5kiSXGs04ibJo5BsEPLceRHbStjxC\nVslFSv65MX/2TQXgt1fWVraHwxs3buwW2N/e74Ltg/I7V79piuOd1ebQpMrublv8+fJ37Q8f2l8+\nNCTzsCGPvya7u/RfV4UIXYWovAeEWSitfrP70FpstFYb7fbcl6uNs2/ayQcvqnsoavzz58/7Bfb+\nXR9sO8rv7H3WfNXvL/utoUmVfp/0+5evLDk7O81LPxiSeUK2xPGi1e/Pe86T/ivyHUTlHSDMT6xX\nvc/6vnW90VxuXGySmx+Rs2/6OR/8ra3hVu+cm3pOP/F4SOZ9w00985JwzvIH19q2MZl5ypMHhETJ\nbUs0vYuEQ1RA93jaiulX1sJFi3xukQb5bSOj41LQ1PO13go/785dZDcIIcx0r54z6kghsgC9xavL\n8Kboa30SUsemsq8dqZZLDQSeCwGLLDRE104YIV9kvKdoOLe9NQnjeP5YsjA+jueiAsb/+9qoDH8s\nazmPd/sx9PMFHMc/fixkhpzHf4r5nzKbFpzAOVEEJ3AQvEEZBK/JCW0qCB7Bm5RB8Jqc0KaC4BG8\nSRkEr8kJbSoIHsGblEHwmpzQpoLgEbxJGQSvyQltKggewZuUQfCanNCmguARvEkZBK/JCW0qCB7B\nm5RB8Jqc0KaC4BG8SRkEr8kJbSoIHsGblEHwmpzQpoLgEbxJGQSvyQltKggewZuUQfCanNCmguAR\nvEkZBK/JCW0qCB7Bm5RB8Jqc0KaC4BG8SRkEr8kJbSoIHsGblEHwmpzQpoLgEbxJGQSvyQltKgge\nwZuUQfCanNCmguARvEkZBK/JCW0qCB7Bm5RB8Jqc0KaC4BG8SRkEr8kJbSoIHsGblEHwmpzQpoLg\n1cBzn7GX6RWJF4BtuCcTPL88f+HbJAkZC8zIhIvz7UvMi5LdJ4z5MB04+N+2Wh8XbIk/DnhOWRg0\nn4gXtnhB1RM4TST4iDQ6V6xPfLsbup4JGb9BO4Rc7NLunWY3nCdCx1VWAYO/RD71GwXJEMYB7/ni\n8NCR9V28iKn+WB2ZefAXGjwZfN1oyFg5gERBysmIKBUtpGsHMW2+Ev/6OBA6ypUFCj4kYTLk5JPc\nNyiBH2RbWx4+tMSLF/LV0sZA0d6rvvHIhvAiQJVGRxR537ggX6939MtsXOkM2s/uuZ3Bxx/vDDY6\nG0LjnrLOzg7Mmw6VMevM577hlzFqvCPTLOw2k8RNP7i2ctaFiazx5Koo8pq05Gvf1y8TXnQTr+sv\nsGRuoZ+EXsD2m0olg9Z4n8iYLc3nvmGcpt53RbO4uix8ssUL39Eeq2NmHvx10Ti+bhMqWvkIknxO\nWcam4s7eoHHgtL5LEtqIRHuvrAMFH4vP8fA/yWruG8bq1TPq2UuvxItAvHDUsypNJHjuWJSQbkRd\nF9BPBchEc4QSyxGB+rHleA3qegAdcOfuNiFz1tX88+ON43kY741JeQhJqjSR4EUt+faPskwESgEJ\nSqn1bZgGarcvNUA6Jcbxd+4UVUWcwDlRZJomcIJHRWcR/Iki0wM+oNatovMI/kSRaQEvsHsxztWr\nF5kO8KEtsOOXNJAi0wA+dFLsCB5SpP7gD7EjeEiRuoOPXcs5nBRC8OpF6g0+9qzjXy8hePUidQbP\nWTrffGQIXr1IjcH7hAYn/wfBqxepLfiQEnZ6hhbBqxepKfhIdOXPTswjePUitQTPvWNd+WOG4NWL\n1BG8f6pPd2gIXr1I/cDLm3vOKQSvXqRu4GM36+a+bwhevUjNwDNS9DwwglcvUivwopUPisogePUi\nNQJf2MqnhuDVi9QHvE/sUc/sIXj1InUBH9nEH1kGwasXqQd4zixX4Yl8BK9epBbgQ6r2QD6CVy9S\nA/DctZja2hUEr15k8sF3R3fqDgzBqxeZdPBiDKe8zhLBA4pMOHgxhgMsVUPw6kUmGnzsAKp7guAh\nRSYZPKy6JwgeUmRywb8CVvcEwUOKTCz4HrS6JwgeUmRCwYux+zJcBsGrF5lM8CGhTyZqg8MyTpSM\n1QyD58xifLJ2toQ7UTpWsws+skk4aVuagp0oH6uZBe9bcm8xBG9WZeLAc2f/e3cEb1Rl0sB3Dwdx\nCN6oyoSBZ0dzNgjeqMpEgY9tcvS8BYI3qjJJ4ANyfK9QBG9UZXLAc+/k1DyCN6oyMeD3Bu/HDMEb\nVZkU8CebeWkI3qjKZIA/3cxLqxT89iZXAx8zBtrxGRIr7jOfVwA+YJ7H1s2CDxlLH4cvBh+fbual\nlQIf7OtlWQH47bWna1sq4H0aBDbkKQFArGLKQka5cfCue4OQG+1rJmVcp+unT8QXgu+SrJ3/y4C/\n4u7rZVkB+HBzuP1GAXxM5YVCnhMAIEmznvieafC+y2kckUE7MCcTyJxTacqmIvAs+0GbEuD9K4nM\nrZVztgB8b23t1vZwePv27aKsNx929jLqmMjfI4ykx6bpZERLzzauiePmJsALqEwnTdZ0bUPGLO89\nr9vkXuYJaDIi6VOcHnNSRBUkI+ptDp/2lMF//bWBWO2Df20c/LV7G0uDQfv59ybBp0iXnhWAf0Za\nL7LPlAB/7b48tp9lny0AHz5NwY9s6jmNEkheHVAjzGTD5zHTTX3X5jTs2q9bsLQkIJk0ZVPXTvKb\n+iB/xXuJpr7bOtDLsqLO3UpvReUen4TUY6d3VdQTK/m0mcNs13yvntGrpHGVrJuU8SnzZCXJA+9Z\nfm7ZMp27m0Ivt+tVOI7fUhzH8zBUT0EFilUic/bIazc+nItDYS/MDueERvozEzy3SUH/uNRw7kAv\ny3AC50SR85vAiYhdVHlw5s6oyvmBF7f3wiII3qjKuYH3rKC4CII3qnJO4HnWJO1JQ/BGVc4HfERH\n52NG8EZVzgV8QNzRYyIEb1TlPMAzpVWwCN6oSvXguVe8LemBIXijKpWDL561OWYI3qhK1eBFt05x\nyhPBG1WpGHxXpVu3ZwjeqEq14APA5iYI3qhKpeBHztYdNwRvVKVC8NwdOVt33BC8UZXqwHObgjYz\nQvBGVSoDP+JL2LOG4I2qVAX+B+LBuCN4syoVgb8z4sv3DEPwRlWqAc+s78BlELxRlUrAe+SOkfzx\nZwzBqxcxD14M4yIj+ePPGoJXL2IcfDqMQ/CanNCmYhx8ZMthHILX5IQ2FdPgo72VsAhekxPaVAyD\n75K9YRyC1+SENhWz4A+fnUfwmpzQpmIUfHD4bRyC1+SENhWT4I99C4vgNTmhTcUgeO/YNpUIXpMT\n2lSMgefe8YcqEbwmJ7SpmAJ/6mFaBK/JCW0qhsBz++QaKQSvyQltKmbAx/appy4QvCYntKkYAX/2\naRsEr8kJbSomwEdnH55H8Jqc0KZiAHxEzj5tg+A1OaFNRT/4bgZ3BK/LCW0q2sFnb22D4DU5oU1F\nN/icRVIIXpMT2lQ0gw9yFkkheE1OaFPRCz6PO4LX5YQ2Fa3g8xdFInhNTmhT0Qm+YDEsgtfkhDYV\njeC9gsWwCF6TE9pU9IH3iva2qTN47rOCnZhLxSrPeMC6yX0WAJcZlgAfdHRtW17EnfurP4BV4OC7\n7G7R6ZLgI8oCJ52bMA4+FlIfkyuBN3r7xzFUEvnlqbd+Dbh6OUemiLsI3R3H+KJJ4UxwvciZkuDT\nVf1OkFQA3ukmCWvdlAmDDKoI83zR1DMwkrMyvLCdT1fSpKGDGBS8cCYZ+gXOKIHvn7YfLsvjE3l8\n/64Pth3Im1vy79aiOBKDKnu/fmen3xzbmVct8iT/7TJ079+loYPYO2CYpRvDImd+LgV+77p7n/Ur\nAC8vvvXnCsA3JfhXY4Mv5p6G7v27NHQQg4KXsRoWRUwJfEZ7JXt2zsjkecqtY4HJVpHRm2n2OXMq\nSZrySDT1kMSJWTIj96oUoRNNfTeBGbSpl34Mi5wp37nzafprzXfubM9faCz6LrTfBe7cOe7XbWfM\nzh23R12mCN03TfDHC9y5c1z/iqO/cyfGWP7eB7uCcXyXRckjBq0jJYZz3c7GmM6M5i5D980TsEqZ\n4dyjotM4gXOiyJgTOCrck3pP4IzjxNSCV+SO4HU5oU1lPPCq3BG8Lie0qYwFXpk7gtflhDaVccCr\nc0fwupzQpjIGeAB3BK/LCW0q5cFDuCN4XU5oUykNHsQdwetyQptKWfAw7ghelxPaVEqCB3JH8Lqc\n0KZSDjyUO4LX5YQ2lVLgwdwRvC4ntKmUAf8azB3B63JCm0oJ8LwFf0wPwWtyQpsKHDy3W2DuCF6X\nE9pUwODF/f01XAbBa3JCmwoUvOzXjZ0/XskQvFEVIPi0P4/g1W1KwO+N4xC8uk0HeG6TuJwMgtfk\nhDYVCPiD56gRvLpNBXh3//l5BK9u0wD+cH0cgle3KQB/tC4Swatb/cEfWw+L4NWt9uCPr4NG8OpW\nd/An1r8jeHWrOfiT+x4geHWrN/hT+10geHWrNXjPOrn+HcGrW53Bs9P7FiJ4dasx+LP70yJ4dasv\n+Ix9iRG8utUWfNZ+1Ahe3eoKPnMfcgSvbjUFn73/PIJXt3qCz8k7gODVrZbgc/LLIHiA1RF8Vv64\nsjIIXpMT2lTywedyR/AAqx/4fO4IHmC1A1/AHcEDrG7gi7gjeIDVDHwhdwQPsHqBj0nhQuipBP90\nNPi333WBuWJgsYq63ZFOjK9yUORF94w3oza8gMoIh972k243ApUyD/6nY9x7ayPBh81lRn3gFanH\nijsOc+W269WAX28yRk/mixu50QlQhlHmNf/L9hgsIZF58Ndbtw7Y89HgOf2xLw5hAjL1WKW5gWRq\nikrAx83XMpHAcdCjN7iByXQdcXhIAnF0A0C5Cpr6+Gbzyl35YvvWUIK/ffv2INc2rn3YGQzWOwOQ\nvVd+Z+u1PJLBYAhTgKkcWOdfZJnOvaP/ed3auwJtMkvPxOEDeSGOL9qAcjs78AAUxuyXjHv8/ZbV\nvDkcrodbK3wU+I5h8OmxMvD/fgq8AnegTHsPvDxOGPhH18mVuz/9JCr9Zq/3m82RTf1L0QbZ5pp6\n+cWITDhXSVMftV5Llw67d0obmcFkApk870dLBoz5gHLmm/qbd9Nb/KO/yOPozl3Q/Cqwoal1AJ07\n2ws8SaKazl2nFfg0OPynQxS63kAZzwlY8z8I67oOpNjkjeNffsVgAxNgrIK9pLIVDeeeMXY0nPNU\nuINlQua/7XMfmFxp8sBP7wSOGvfpnMCZZfCK3BE8wOoAXpU7ggdYDcArc0fwAJt88J6lPEJF8Oo2\n8eD97AdqdckgeE1OaFPZBx8AuCN4gE04eBB3BG82VhWCh3FH8GZjVR14IHcEbzZWlYGHckfwZmNV\nFfjA8s3LIHhNTmhTGWxYoAeiSsogeE1OaFN5RsDcEbzZWFXzIAa5VokzCF6TE5pUBPcx8scDDMFr\nckKPSkS8MfLHQwzBa3JCi4pcKIXgZw98ukAOwc8c+L2FkQh+1sBHewsjEfyMgef7C2IR/GyBP1w4\ngeBnCvzRghkEP0vgjy2UQvAzBP74AjkEPzvgTyyMRPAzA/7kglgEPyvgTy2ERvAzAv70AngEPxvg\nz2x8gOBnAvzZDS8Q/CyAz9joBMHPAPisDW4Q/PSD5zY9u7ERgp968NzOWgCP4KcdfDZ3BD/t4HO4\nI/gpB5/HHcFPN/hc7gh+qsGLcVzexkYIforBF21Qi+CnF3zhxsQIfmrBF29IjeCnFfyIjcgRvAr4\nKIz3nYhD2Ib14A2fQ64HfBF3IZJsbIzen35cZ4S9/SHcC566AcDH4X7fdS9mOUqlwXPXZpSl4MVP\n14b4AVzVRl1G74Kjm6FSwD2iHms0ljrHdqs34oy0oPl3DXKVgnb5VwcvaDhO6qbkJrM5ZSqVBp/m\njnB9Ab4roxlAtt2HxYqKDzBvgnOcnVUp4M7lyJ6SF4Nj+SmMOJPIPHYvL/sR5U4AKKQMPk2B4acb\ne0huglCSZCmVBk9SHxwB3k3beUiVB8UqTJ24C82CcValsJ335F+/M0gCc+k29o35u80k8cIIUleU\nwTspByoPw4MXWUpK4LOS2TTTZD3twYedNL/OXnolRQPl79lIMx3dB+Y7OqtSmFhIimx01lffD+5B\nhcA5jzrrH1riuGEmGVE7dTJN4jQ8eJGl9EvZGu+k+XQ8UeN9WUciaqqSxOn35ovAfEdnVKLC/nxM\neMLpwsYgcWApQ0rU+NDZbYacxAzStijX+PS3hmkVT2t8ePB/p6x0Ux9RFnoilqJz57hhAMo5CEzc\nRIOucx0a3lMqUXGC2FRkjqzfg+V/hDsjzWt/1SBfeAoZro5MGTy3vZDRtF8vuR1gOmPlh3OcMZ/v\nDeeC4yl89McqYqw75nBuFHeZI4iFUacDre+lhnPfL3uXQcnHQMM5n7E9X9OYiYYlU2k2JnAi4ilV\nL5zAmS7wkeq+hQh+qsArc0fwUwU+VN+nFMFPEfgAsC8xgp8e8IEFGC4j+KkBD8wvg+CnBDw0vwyC\nnw7wDJpfBsFPBXgPnF8GwU8DeI9Av9lB8FMAnl9Tzgt8dGEIvvbg85fLFF0Ygq87eG7bz0pcGIKv\nOfjYtvlY+eONOoPgNTlxxtKv3xF80cmpBB8Rl4+VP96wMwhekxOnrLv3dRyCLzo5heAPvo5D8EUn\npw/84fQ8gi86OXXgj6ZpEXzRyWkD75HD52QRfNHJ6QLPnWPTdQi+6ORUgT+5uw2CLzo5TeAjejLr\nQIkLQ/A1BL83bTOWCoKvIfgzT9Mi+KKTUwPeP/M0LYIvOjkt4DOeskLwRSenAzx3M56yQvBFJ6cC\nPLdp1ubzJS4MwdcJfM7idwRfdHIKwHdzFr8j+KKT9QefuzgOwRedrD34/EUTCL7oZM3Bcyd/0QSC\nLzpZb/BRfs4BBD/F4EPiFO1BXuLCEHwdwI/Y6wLBF52sMfhRa2ERfNHJ2oIv6taVVkHwEw8+oiO3\ny0bwRSdrCr576qELTSoIfsLBM5WtrBB80ck6gucuUdlrGMEXnSwC/2aTjwLPw7dmwcdhkubW2b8G\nmVgnyvwSdiyVRObsEYfXG7LjEP2joQQ7PNxLFFQV+KggTgXgt1a21p4Wg2fUay7Dr0g5VtyxHeuC\nY0d74H3q2e5thds7TCWRHybHJSFrLlHvQYNQ66L6XvLqMoxcsOiCcKYa8FHTdfOrSAF4Ud+3eoXg\nmSc+vZd98CUpx8oOYvonGgV2Cj7NetRS3qkSUhVlnoWIXHw9SK6SJhcvF/U707UfU7n1vl0N+DSt\nUkTzPsGF9/jtW6LGP3/+vJ9j5FW///6vzT7UdhTf96TVX/6u/91y//IPQ/HPy0/6rxzSeKVZRVjv\nM3mkyzs7/Z71RLxcJqoq6jLOk+Wv+v1m/7Mf3r8Dh+wduMhXyzJmIn7Z9nMB+M21LXG8cePGbo41\nxd+/vSe7UPug+L6Hn++ufr/7cHV39eFQ/LP98sdm68f2W80qwr5flce51Q8fdr+3XoqXd6iqirpM\n+61wZrctnPnbe3DI3oOLrN6RMftmNef0Tj74p+ujevWu6Fzv9lxwu6XaOnLKAzfxfPFTXgNbsDyu\nnvUI0oVMW8TGp6JXf8eS6W9aDe3OJIwJZ0TbS+NKmvrIljHL/fayoKlfX1lb2ywEH1E/+kZHKsA8\nC+yuMzfXtX15j+dXrYWoS5WzxkA6d8wJI/cq+fpZQJfI5TsXlMaLQBlO2Ud2o2P7FXXuvMUodHI7\nROON42PmrL4CXxEASeg5bdcL5XBOjOIe+A5T37wONJzrum6Q8E6bxcm/zTfagC3y1GW477Q/8rqV\nDefuum7+57c2EziB6ihuHBWcwJk08Pw6JN9AWRUEP2ngIxu8HzGCnwLwAXF+Mq+SIPjJAs89yx8z\n4aByEQQ/OeAjKiecETxcpt7gfSvtzSN4uEydwXOHBKOdGFflqAiCnwzwXXIw5Yjg4TK1Bc+PPWGF\n4OEydQV/YvCO4OEyNQXPrONztAgeLlNL8LG936tTcaK0ypkiCP6cwfvEOfltL4KHy9QPfOxYPsSJ\ncipZRRD8eYI/GsQpOlFKJbMIgj8/8NzN+gYWwcNl6gU+IJmPiSF4uEydwIu7e/YDFwgeLlMj8D7J\neyoUwcNlagM+dkju81UIHi5TF/Ds9Nhd2QmISnERBF85+NAmflkn1FVGFUHwFYPnJ2fmgU6oqowu\nguCrBd+loxbIIHi4zMSDl2O4UYslEDxcZtLBi07d6BVLCB4uM9ngQ3ry+9cyToxWUS2C4CsCH7uj\nW/nRToxSUS+C4CsBz5Va+dFOFKtAiiD4KsAHaq38aCcKVUBFELx58JFCX17RiQIVYBEEbxo89ywH\nspUGgofLTCJ4Rihs3TOCh8tMHvh1WjgvD3Yi2xB80clzAB/aBHBzV3Ei2xB80cnKwYs+nfdCsxPZ\nhuCLTlYMPhZ9ushM/vgzhuCLTlYKnjPLCcvFCsHDZSYFPBdd+aBsrBA8XGYywB9hR/AzBP44dgQ/\nM+AF9hPPzyJ4cJE6gk+xnxi4I3hwkfqBl4386fkaBA8uMoHgX8VJ/jzcyXv7GLHavwZDWYIOiwwS\nDp1ULAseJlQaPOfZQuOC5w3LskhONu/YszKwA2IVUofuPXadXgMjDvWV/S4D/hm1bdC3hmXBc1e4\nBvieqiR4IdOw7Kzd98cF3yCv7GvWr8KMpDcCux2MFatY7mnpO/tOJMzjCbeNJCrYt9fNUGZHAJUp\nB95lMvuR+rb4JcE7vssi+vuMZFRjgo+slz0v+cxKmH/qTOjszdKNEysWyKN9kHcu/XDFjqrfJYh8\n3ZFHB/SNcSnwcfrhCtS3Yi8HPrKlUMC6Z/OtK4Ef5NqG9eEPncEfrcG9zon/X29Z1zZyS70fqFln\n4/Aor6GV/mdbsbCyyjFb/f5I1qDMh52Na/LnM2VfBjs7cJnhYKMjhTY6GUK/jFfjY+sfntjJHEmO\nJ8HgjBKv6EapWkl8+UHl9GAv27TBCpRzH5WoiutL8khBN/lSNZ4T6dSZZjLfytV4TmWaJc/PEBr3\nHu9Y7hKx/okdNcCxd3b8VjZWNoujvXw68hpC2uWBOpUynbv2oaB6GbiMuMczO4qZrd6xL3mPZ87n\n85/MZwmNPZz7e2KRecc/+GdX3NoDbbHizNlvStJriDzHU6+NpYZzzBl9+WPLyOFc13Ugj6OUHc51\n3fm5hSwhrRM43KeWp9AzwgkccJEJnMA5vKJIjNrVPsQIHlxkYsHzwFZo48vHCsHDZSoAH3mEMKN3\nXwQPlzENnvuAyl42VggeLmMW/I+fWZDKXjZWCB4uYxj8ZfX8u2PECsHDZc77QQwtsULwcBkEr14E\nwSN4gzIIXpMT2lQQPII3KYPgNTmhTQXBI3iTMghekxPaVBA8gjcpg+A1OaFNBcEjeJMyCF6TE9pU\nEDyCNymD4DU5oU0FwSN4kzIIXpMT2lQQPII3KYPgNTmhTQXBI3iTMghekxPaVBA8gjcpg+A1OaFN\nBcEjeJMyCF6TE9pUEDyCNymD4DU5oU0FwSN4kzIIXpMT2lQQPII3KYPgNTmhTQXBI3iTMghekxPa\nVBA8gjcpg+A1OaFNBcEjeJMyCF6TE9pUEDyCNymD4DU5oU0FwSN4kzKTDv7N5psi8HeIZVlLbz1L\n/JR/rirvwq0eq9hJkx21GlLLIsSm6pvJQ4hwRm0aJoPXHmmQeVB+e2WZyKHUS39xNeDDlk2uXyXU\nycx8UwD+zdrm2pt88I8t0k0WCLUesF9Z5Ft6dU5/thhOfdr4dH6+QcSn7HKLNBn3lMlDwHuCSETj\nQcuPaYdyppwGQ12G0+5BZqVKwMf0UZI0FpKkS7I+xgXge0+Hm2E++IW08hHry4Rcv2AlXdeeV02q\npIzEZ97HQbLQII2rF254DZ+S/UwlWlVkkOQxYBtLMgES8xNb/47MfvqJTatfJeBZMExCRyb0yMyA\nUgR+a7jVGw5brVbmjeCic18erQfDxZu/toY/ieOjoWa7eX9x+dFw2W4utS7832LzwRVxKYs/6VYZ\nDh9dkcd48dHN4c1Hw/vpUbsvj46OVdhiLPy6eeWRPGac3k2PReCHOTX+EvVE+2VZnyZkdc5KgqtU\nOV2Mcl0MXK/NEqdByHLrEiO/IVTUeNXCgBq/14z47FlbVkw3UG9W1GWCNAUYra7G+0Nx95KOeEHG\n6YIaH27KP7nguUWu/6HRuGB94VOLXGosNA3cfW2PktaiaOqXiTW/ZBE/tn39KjJpT5z4lA+WvJgu\n2LELyEcEyKzE939xJeA5uZXEjWbMWWYCxQLw22u9te188Ek8TyzCdu9SQqiok01ff6xEb/uC+NXz\nSwu0ISUuOI76XtmgcVbXcUT/TiYjWphfOMqtpFOGM2c/wVU1vXp+XURLasJzy27t/yz69TiOhxeZ\n9HH8oWl2AsEjeD1OaFNB8AjepAyC1+SENhUEj+BNyiB4TU5oU0HwCN6kDILX5IQ2FQSP4E3KIHhN\nTmhTQfAI3qQMgtfkhDYVBI/gTcogeE1OaFNB8AjepAyC1+SENhUEf9L6BfbfL/tgG8CLfFeJyp//\nXIlMiZhtb8NlCmP2iwr4IrvxvHxZgLXG/xUKdvt2JTKTEzMEv2cIXt3+56+VOFENkefVEJmcmI0B\nHq3OVh78m01exQXyg/Wbpu1pBRq9Xq8Kb7YV0JQGv7WytVZBrLZuba1UQr63Zl7jzVol4LfXnq5t\njXpTafDiQ7W/zMqoCZneSCc0GK8C/FZva7sCX8LN4fbIz9cY9/jtW1W0jsP1lQqCtX1rWAH4cKW3\nUsH9sbe2dmtkzMqD3xzdnGjyIzSvsR5uVYFE9CQqaCV7mwoypcE/XTfvgXTiTSWx2uz1frNpXOVp\nNc6ET02CX19ZWzMfquFT0TpW062v4h6/1luronOnErOJH8dvV3RDmSpnFGQmHjyaGUPwM2oIfkYN\nwc+oIfgZNQQv7f6t4V+unPdFVGsIPrVWvPjovK+hWkPwqcXk5nlfQsWG4FP7C4KfTVt81IrP+xqq\nNQQv7ebNYVzNQ50TYwh+Rg3Bz6gh+Bk1BD+jhuBn1BD8jBqCn1FD8DNq/w9ZrTU/tSqPTgAAACV0\nRVh0ZGF0ZTpjcmVhdGUAMjAxNi0wNy0xN1QxMTowODowMCswMDowMGunxTUAAAAldEVYdGRhdGU6\nbW9kaWZ5ADIwMTYtMDctMTdUMTE6MDg6MDArMDA6MDAa+n2JAAAAIHRFWHRwZGY6SGlSZXNCb3Vu\nZGluZ0JveAA1MDR4NTA0KzArMKV3vKMAAAAUdEVYdHBkZjpWZXJzaW9uAFBERi0xLjQgHEc6eAAA\nAEp0RVh0c2lnbmF0dXJlAGFlMGNhNDIxMDIzMTcxN2Y2YzU0OGZjNmQ2YzBkZGE5Y2RjYzQ0MjAx\nNDhlZjY5MzNmNmMzNTQ3ZmQ0YzZhNDMqQYCqAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 結果からロジスティック関数をmyfunc <- 1/(1 + exp(-(-4.190 + 1.005*z))*Nで定義\n", "r('myfunc <- function(z) { 8*(1/(1 + exp(-(-4.190 + 1.005*z))))}')\n", "# 結果をプロット\n", "graph = preGraph(\"images/fig7.10-R.pdf\")\n", "r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=position_jitter(width=.2, height=0), shape=21, size=2.5) + stat_function(fun=myfunc)')\n", "r('plot(p)')\n", "postGraph(graph)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

$r_i$の分布

\n", "\t

\n", "\t\tもう一つGLMMから求められた$r_i$の正規分布をプロットすると以下の様になります。\n", "\t\t-10, 10ではほぼ０になっています。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEBCAYAAACaHMnBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XlcVOX+B/DPOXNmYdgEFbymqGRqCkqodU1tcc0FUHMN\nb4lmZZplpt5+XbesvJaVlpaapWma1yUFb1hamZlZbqm5lQtqZi4oCgwwyznn98dcLBJlm5kzy+f9\nevECDsPMM8+cOZ95vs9ZBFVVVRARUUATtW4AERFpj2FAREQMAyIiYhgQEREYBkREBIYBERGBYUBE\nRGAYEAEAVFVFbm4ueNgNBSqpMv908WJeqctFUUBkZDAuX7ZAUfimcjX2r/tYLHlo0OAWZGX9huDg\nUK2b43e47rpPWX1bs2b51meXjgxEUYAgCBBFwZV3S//D/nUfQRBKfCfX4rrrPq7qW5aJiIiocmUi\nIn+gqsDPP4v47jsdCgqcb4XFiyW0bi3ijjsUSHx3UADh6k4BR5aBpUv1eO89PY4e1UGSVEiSFQDw\n8ssG2O3BqF1bwahRNqSm2hEUpHGDiTyAZSIKKD/+KKJbNzPGjzchLk7BsmUFyMrKx5EjBQCAI0cK\n8OmnFrRtK2PiRCPatAnG5s06jVtN5H4MAwoYixbp0b27GXY78OmnFsyfX4TOnWUYjX/cRqcDWrdW\nMHduEbZts6BRIwUDBwZh/nw9uNcp+TOGAfk9VQWmTjViwgQT0tLs2LSpAK1bK2X+X2ysihUrCjFy\npA0TJ5owYYIRsuyBBhNpgHMG5PdmzDBg7lwDpk0rwuOP2yv0v6IITJpkQ2ysiueeM0IUgenTreAe\nqORvGAbk1+bN0+ONN4yYNKniQfBngwfboarA2LEm1Kql4plnbC5sJZH2GAbktzZv1mHKFCNGjbJi\n1KjKB0Gxf/zDjvPnBbzyihF16ijo29fhglYSeQfOGZBfOnVKwOOPB+H++2W88ILrPsWPHWtDv352\njBtnwokTrBWR/2AYkN9RFGD0aBPCwlS8+24hdC7cM1QQgBkzilCzpoonnwyCveoDDiKvwDAgv7N4\nsR7bt0uYNasI1aq5/v5DQoB33y3Evn0iZs40uP4BiDTAMCC/cvq0gBdfNOKRR2xo1859+4G2bKlg\n/HgbZs0y4PvveVAa+T6GAfmN4r19IiJUTJpkdfvjjR5tQ+vWMsaMMcHq/ocjciuGAfmNjz+WsGWL\nhNdfL0KoBy5JoNMBr79uxenTAt55h+Ui8m0MA/ILV68CU6ea0L+/HR06eO4w4caNFQwbZsfs2Qac\nP8+9i8h3MQzIL7z9tgFWKzBxoufrNWPHWmEyqZg+naMD8l0MA/J5v/8uYMECA554woboaM+fTS48\nHBg3zoaPP9bjp5/4liLfxDWXNHXmzK9ITe2HJk3qo1WreEybNvmGt7VYLBgx4lFER4fj+PGj15bP\nnGmAw3EVb71VG/XqRSMmJgoxMVHo0KGdJ54CAODhh+1o2FDB5MlGnt2UfBLDgDSVljYYtWvXwa5d\nB7BqVToyM9dj/vy5193u/Plz6Nz5Huj1+hLXKT56VMSyZXrUrbsEs2ZNx6lT53H69AWcPn0BX331\nrceeh17vLFF9+62Ebdu4qyn5HoYBaWbv3j04dOgAJk2aipCQEDRoEIsnnhiFpUsXX3fb7OxsTJ78\nEsaNex7qnz56v/yyAXXqqKhVK73Eci107SqjeXOZB6KRT2IYkGb279+HunVjEBoadm1Z8+YtcOzY\nUVgslhK3bdYsDl27diuxbM8eEZmZekyYYIUo2rFu3Rq0b38nYmNvQb9+KTh5Mssjz6OYIDjPXfTd\ndxK++46jA/ItlTprqSgKEMXrd6PT6cQS38m1/K1/r17NQbVqEZCkP55PjRrVAQC5uTkID7/+YIHi\n2+p0IubMMaJhQwX9+yv48cfbYTYHY+HCRVAUBePHj8WgQQ9i+/ZdkMpxZfs/9+2f21NRPXsqiI+X\n8frrRtxzT1Gl78ff+Nu6601c1beVCoPIyOASddu/CgvjFcTdyV/6NyjIAJ1OQERE8LVlly+bAQDh\n4eYSy4vl5jr/fulSOD79VMKCBUCNGsF47735JW63ePEHiIyMxIEDe3D//feX2RadznlsQlhYEMLC\nrn/cipg6FejTBzhwIBjt21fprvyOv6y73qiqfVupMLh82XLDkUFYWBBycwshy2VfVpAqxt/612wO\nw8WL2cjJ+aMklJV1BoIgQJLMJZYXu3rVeeH6efPCEB2toGfPQuTklHbvAiIiInD0aBYSEu4ssy0W\nSyEA/K9vq1biueceoGnTIEycqGLtWo4OAP9bd71JWX1b2oeq0lQqDBRFhaLceLJOlhU4HHzB3cVf\n+jc+PgFnzvyKixezERERCQDYtWsnGjVqAoPBVOpzdC77GzIzI/DPf9qh0ym4ciUP06ZNxrPPTkB0\ndDQA4NKlS8jOzkadOjHl6qviN5Gr+vbZZ6149NEg7N4NtGjh+6+Vq/jLuuuNqtq3LOCRZuLjmyMh\nIREvvTQF+fl5OHr0F8ybNxdpaY8CANq2bYUdO34o8T+qqkJVR8NgUPHII86L1oSEhGL37l34v/8b\nhytXcnDlSg4mTHgWcXHN0br1XR5+Vk7duzsQE6Ng3jzuWUS+gWFAmvrgg6X4/feziIu7DX369MTA\ngakYMmQYAOD48WOwWPIBAG+++RpiYqLQtm1HAE+gqGgW4uOjMWvWTADAkiUfQ1VVtGmTiMTEOMiy\njGXLVmr1tCBJwPDhNqSnS/j9d56ziLyfoFZi5+yLF/NKXS5JIiIigpGTY+FQ0A3Yv8CcOXr8+99G\n7N5tcempJwoK8lG/fm2cPHkWZnOIS+4zLw9o0SIEQ4fa8K9/ue7Sm76I6677lNW3NWuW7xS+HBmQ\nz3A4gIULDejXz67JOYgqKjQUGDzYjiVLDLBcPxdO5FUYBuQzNm6UcPasiKFDfefCw48+akNuLrBi\nhV7rphDdFMOAfMaiRXq0aiUjPt53ygwxMSq6d3dg8WI9T2BHXo1hQD7hxAkBW7ZISEvzvdr7kCF2\n/PyzDj/8wFNUkPdiGJBPWLzYgMhIBUlJDq2bUmHt2smIjVWweDFLReS9GAbk9QoKnDX3hx6yw2TS\nujUVJ4rAww/b8N//SsjO5m6m5J0YBuT10tMlXL3qvICMrxowwAFBAFasqNRB/0RuxzAgr7dokQEd\nOsioX993Z2CrV1eRlOTAkiUGKL4z/00BhGFAXm3vXhF79+p8cuL4rx55xI6TJ0V88w0nksn7MAzI\nqy1dqscttyjo2FHWuilVduedMm6/XeZEMnklhgF5rfx84JNP9Bg0yA6dH3yYFgTnvMfGjRIuXOBE\nMnkXhgF5rfXrJRQUAIMG+e7E8V/16WOHKAKrV3MimbwLw4C81sqVerRrJ6NuXd+dOP6riAigWzcH\nVqzgEcnkXRgG5JV+/VXAtm0S+vf3n1FBsUGD7DhyRIe9e/n2I+/BtZG80urVepjNKnr08L0jjsty\n770yatdWsHw5J5LJezAMyOuoqrNE1KOHAyGuubSAV9HpgP797Vi7Vo8iXiKZvATDgLzOnj0ijh8X\n/bJEVKx/fztycwVs2sSJZPIODAPyOitX6vG3vylo1873jy24kYYNVdxxh8y9ishrMAzIq1itwLp1\nevTt6x/HFtxM3752fPGFhMuXtW4JEcOAvMwXX0jIyRHQr5//TRz/VUqKA4oCZGRwIpm0xzAgr7Jy\npYQWLWQ0aeL/Z3OLilJx330sFZF3YBiQ17h0ScAXX/jnsQU30q+fHTt2SDh5kqenIG0xDMhrrFsn\nQVWBXr38v0RU7IEHHAgOVrFmDUtFpC2GAXmNVav06NhRRs2agXOeBrMZ6NHDgdWreXoK0hbDgLzC\niRMC9uzRoV+/wCkRFevb147jx0WenoI0xbWPvML69c7TT3TqFDglomLt28uIjlawahVLRaQdhgF5\nhfR0CV26OGA2a90Sz9PpgD59HFi3ToIj8LKQvATDgDR3/LiAAwd0SE4O3C1h7952ZGeL2LbNz4+0\nI6/FMCDNZWQ4S0QdOwZuGLRooaBePQXp6TzmgLTBMCDNZWRI6NrVgaAgrVuiHUFwjg7++189bDat\nW0OBiGFAmjp2TMDBg4FdIiqWkuLAlSsCvvmGpSLyPIYBaSojQ4/gYBUdOjAMmjZVcNttMtat415F\n5HkMA9JUejpLRMUEwXn0dWamxIvekMcxDEgzR4+KOHxYh5QUjgqK9erlQH6+gK++4kQyeRbDgDST\nkSEhJETF/fczDIrddpuCpk1l7lVEHscwIM0U70VkMmndEu/Sq5cDn38uwWLRuiUUSBgGpImffy4u\nEQXeuYjKkpJiR0GB83TeRJ7CMCBNZGRICA11XtyFSmrQQEVCgox16xgG5DkMA9LE+vUsEd1MSood\nX34pIS9P65ZQoGAYkMcdOSLiyBGWiG4mJcWBoiIBn33G0QF5BsOAPI4lorLVqaOidWsZ6ek8AI08\ng2FAHpeRIaFbNweMRq1b4t1SUuz4+msdcnO1bgkFAoYBedSRIyJ++UWH5GSWiMrSs6cDNhtLReQZ\nDAPyqPR0CWFhKu69lyWistSu7SwVrV/PUhG5H8OAPEZVWSKqqORkOzZvZqmI3I9hQB5z+LCIo0e5\nF1FFJCWxVESewTAgj8nIkBAeruKee1giKq/atVW0asVSEbkfw4A84s8lIoNB69b4FpaKyBMYBuQR\nhw6JOHaMJaLKKC4Vff45S0XkPgwD8oiMDAnVqqlo354looq65RYVLVvKWL+eYUDuwzAgt1NVID1d\nzxJRFThLRTxXEbkPw4Dc7uBBESdOiCwRVUFSkgNWK0tF5D4MA3K7qpSIzpz5Famp/dCkSX20ahWP\nadMm3/C2FosFI0Y8iujocBw/frQqTfY6deo4S0UZGQwDcg+GAblVcYmoRw879JXYOzItbTBq166D\nXbsOYNWqdGRmrsf8+XOvu9358+fQufM90Ov1EATBBS33PklJLBWR+zAMyK0OHBCRlSUiKani1zne\nu3cPDh06gEmTpiIkJAQNGsTiiSdGYenSxdfdNjs7G5Mnv4Rx456HqqouaLn3KS4VbdzI0QG5HsOA\n3CojQ0JEROVKRPv370PdujEIDQ27tqx58xY4duwoLH+5QHCzZnHo2rVbldvrzerWVZGYyFIRuUel\n1ipRFCCK1w/FdTqxxHdyLV/rX+eBZnr07OlAUFDF23z1ag6qVYuAJP3xvzVqVAcA5ObmIDw89Lr/\nKb6tTieW+L+y/LlvK/J/npaS4sArrxhQWCgi9Pqn77V8bd31Ja7q20qFQWRk8E3rsmFhQZVuEJXN\nV/p3zx4gKwuYP19ERETFJwyCggzQ6QRERARfW3b5shkAEB5uLrG8WG7uzf9+Izqdc+QSFhaEsLDy\n/5+nPfwwMHkysG1bMAYN0ro1Fecr664vqmrfVioMLl+23HBkEBYWhNzcQsiyUqWG0fV8rX+XLNEj\nMlKPhIQC5ORU/P/N5jBcvJiNnJw/SkJZWWcgCAIkyVxiebGrVwuufS/t7zdisRQCwP/6VlfxxnpI\neDiQmGjC8uUqHnjAqnVzys3X1l1fUlbflvdDUaXCQFFUKMqNJ+lkWYHDwRfcXXyhf1UVWLtWQo8e\ndgAKHBWfP0Z8fALOnPkVFy9mIyIiEgCwa9dONGrUBAaDqdQ+cDgUCIJQ4T4qfhP5Qt8mJdkxY4YR\nV64oCAnRujUV4wv966uq2rcs4JFb7N8v4vRpEcnJlUiB/4mPb46EhES89NIU5Ofn4ejRXzBv3lyk\npT0KAGjbthV27PihxP+oquq3exMVS0pyoKhIwKZNnEgm12EYkFukp0uoXl1B27ZVOxfRBx8sxe+/\nn0Vc3G3o06cnBg5MxZAhwwAAx48fg8WSDwB4883XEBMThXbtWkMQBNx/f1vUqxeNWbNmVvm5eJuY\nGBV33MG9isi1uDaRyxXvRdSjhwNSFdewWrX+huXLV5f6t3Pnrlz7ecyYcRgzZlzVHsyHJCXZ8eqr\nRuTnw+dKReSdODIgl9u7t+olIro5lorI1RgG5HIZGXrUqKHg7rt5ump3qVdPRUICS0XkOgwDcqni\nK5q5okREN5eU5MCXX0rIz9e6JeQPGAbkUj/+KOLXX0WkpLBE5G7JyXYUFQn44gumLlUdw4BcKj3d\nWSJq04YlInerV09FixYsFZFrMAzIZVQVWL9eQs+eDui89yBev1JcKrKU/2BrolIxDMhl9uwRceYM\nS0SelJxsR2EhS0VUdQwDcpn0dD1q1lTw97+zROQp9euraN6cpSKqOoYBuYSisESkleRkB774gqUi\nqhqGAbnEzp06/PabiF69WCLytKQkZ6noyy85OqDKYxiQS6SnS6hVS8Fdd7FE5GkNGqiIj2epiKqG\nYUBVJsvOA81SUhwQuUZporhUVFCgdUvIV/GtS1W2fbsOFy6ISEmxa92UgJWUZEdBAUtFVHkMA6qy\ndesk1K2roGVLXrREK7GxKuLiWCqiymMYUJXY7cCnn0pISbHjJpfFJg9ITnZg0yaWiqhyGAZUJVu3\n6nDpEvci8gbJySwVUeUxDKhK0tP1aNBAQXw8S0Rai4117lW0di3DgCqOYUCVZrUCmZkSevdmichb\n9O5tx6ZNEvLytG4J+RqGAVXa11/rcPWqwHMReZFevRywWgVkZnJ0QBXDMKBKW7dOj8aNZdx+O0tE\n3qJOHRV//7sDa9fqtW4K+RiGAVVKYSHw2WcSJ469UO/eDmzZokN2Nmt3VH4MA6oU5zn0BfTqxQPN\nvE1SkjOg169nqYjKj2FAlbJunYS4OBm33qpq3RT6ixo1VNx7L/cqoophGFCF5ecDmzaxROTNeve2\n4/vvJZw5w1IRlQ/DgCps40YJhYUCkpNZIvJW3bs7YDKpWLeOowMqH4YBVdjq1Xq0bi2jfn2WiLxV\naCjQuTP3KqLyYxhQhVy4IGDzZh369uWowNv17u3ATz/pcOwYS0VUNoYBVUh6ugRRBE9X7QM6dXIg\nNFTFJ59wdEBlYxhQhaxerUfHjg5ERmrdEiqLyeScO1i7Vg+VFT0qA8OAyu3YMQE//qhDv37ci8hX\n9O5tx/HjIn76iW91ujmuIVRuq1frERamonNnhoGvuOceGTVqKFi1iqUiujmGAZWLqjrDIDnZDpNJ\n69ZQeUkS8OCDDqxZI8HOaR66CYYBlcuOHTqcPi2ib1+OCnzNgAF2ZGeL+OorndZNIS/GMKByWbVK\nQp06Cv7+d1nrplAFxcUpaNZMxn/+w1IR3RjDgMpktQIZGXo8+KAdItcYnzRggB0bN0rIydG6JeSt\n+NamMn35pYQrVwSWiHxYnz4OyDJ4RDLdEMOAyrRqlYT4eBmNG/MiNr4qKkpFx44yVq5kGFDpGAZ0\nU1euOM9QytNP+L4BA+zYs0eHo0f5tqfrca2gm/rkEz1k2VlmIN/WpYsD1aqp+M9/eCZTuh7DgG5q\n+XI9Ond2IDqa5zPwdUYj0KuXHatWOQOe6M8YBnRDP/0kYv9+HR56iCUifzFggB2//y5i61Yec0Al\nMQzohpYv1yMqSkGnTvwY6S8SExU0bMhjDuh6DAMqVWGh8/QTAwfaIbHE7DcEAejf34HMTAl5eVq3\nhrwJw4BKlZkp4epVgSUiP9Svnx1FRc4DCYmKMQyoVMuX69GmjQOxsZw49je33KLivvtkfPQRw4D+\nwDCg65w8KWDrVomjAj/28MN27N6tw4ED3ASQE9cEus6KFXqEhqpISuKxBf6qSxcHoqIUjg7oGoYB\nlSDLzjDo08cOs1nr1pC76PXAQw85jzkoKNC6NeQNGAZUwtdf63D2rIjUVJaI/F1qqh15eQIyMri7\nGDEM6C8++kiPpk1ltGjBk9L5u3r1VNx3nwMffmjQuinkBRgGdM25cwI+/1zC4MF2CILWrSFP+Mc/\nnBPJBw9yUxDouAbQNUuX6mEwAP37s0QUKB54wIGaNTmRTAwD+h+73RkGffvaERbmucc9c+ZXpKb2\nQ5Mm9dGqVTymTZt8w9u+9967uPvulmjYsC6Skx/A/v17r/2tV6/uuOWW6qhXLxoxMVGIiYlChw7t\nPPEUfBonkqkYw4AAAJ99JuHcORFpaZ4dFaSlDUbt2nWwa9cBrFqVjszM9Zg/f+51t/v88w2YOfPf\neOed93Dw4DF07vwAUlP7o7CwEAAgCALefHMOTp06j9OnL+D06Qv46qtvPfpcfNXgwXbk5QFr1nB0\nEMgYBgQAWLRIj7vucqBZM89NHO/duweHDh3ApElTERISggYNYvHEE6OwdOni6267dOkiDBw4GAkJ\niTAajRg16mkIgoCNGzdcu42q8mjpyqhXT0XXrg4sXKgHuzBwMQwIP/8s4ttvJQwd6tlRwf79+1C3\nbgxCQ/+oSzVv3gLHjh2FxWIpcdt9+/aiefMW134XBAFxcfH48cc915atW7cG7dvfidjYW9CvXwpO\nnsxy/5PwE8OH23H4sA7btvHU1oGqUjsYi6IAUbx+dxOdTizxnVzLXf374YcGREUpSElRIEmee+2u\nXs1BtWoRJR6zRo3qAIDc3ByEh4deW56TcxmRkZElbhsZGYmcnMuQJBFNmjRBcHAIFi5cBEVRMH78\nWAwa9CC2b98FqRynXf1z33qyD7zFffepaNJEwfvvG3DffVaX3z+3De7jqr6tVBhERgZDuMm+h2Fh\nQZVuEJXNlf2bkwN8/DEwZgwQHR3ssvstj6AgA3Q6ARERfzzu5cvOw57Dw80llgNASIixxDKDQYIk\nSYiICMbChQtK3Hbx4g8QGRmJAwf24P777y+zLTqd85oNYWFBCAvzbD94izFjgBEjRFy5IqFBA/c8\nBrcN7lPVvq1UGFy+bLnhyCAsLAi5uYWQZR605Gru6N/Zs/Ww2/VITS1ETo5nC8ZmcxguXsxGTs4f\nJaGsrDMQBAGSZC6xvHr1Gjh9+myJZefOXUDTps1KLPuDgIiICBw9moWEhDvLbIvF4pyIdvZtYJZK\nevQAwsLMeOMNO1580bUlQ24b3Kesvv3rh6obqVQYKIoKRbnxhkOWFTgcfMHdxVX9a7MB8+dL6NfP\njurVZTg8fF66+PgEnDnzKy5ezEZERCQAYNeunWjUqAkMBlOJ59iixR348cc9ePDBAQAARVGwb99e\npKY+gitXrmLatMl49tkJiI6OBgBcunQJ2dnZqFMnplx9VfwmCuR112AABg+2YelSA8aOtSLYDQOk\nQO5fd6tq37KAF8DWrnXuTvrEE9ocZBYf3xwJCYl46aUpyM/Pw9Gjv2DevLlIS3sUAHD33S2xY8cP\nAIAhQ4Zh5cqPsXv3ThQWFuKNN16FyWRCp05dEBISit27d+H//m8crlzJwZUrOZgw4VnExTVH69Z3\nafLcfFVamnM309WruZtpoGEYBChVBd5914BOnRxo3Fi7T2offLAUv/9+FnFxt6FPn54YODAVQ4YM\nAwCcOHEcFks+AKBDh0544YUpGD58CJo0qY+tW7dg+fLVMBqNAIAlSz6Gqqpo0yYRiYlxkGUZy5at\n1Ox5+aq6dVV06+bAggV6KPwAH1AEtRI7Z1+8WPrFUyVJREREMHJyLBwKuoEr+3fLFh369TNjzZoC\ntG/PC94XFOSjfv3aOHnyLMzmEK2bo6mdO0X06BGMxYsL0b27a2qH3Da4T1l9W7NmaCn/dT2ODALU\nO+8YEB8vo107BgGV1Lq1gjZtHHjrLQMPQgsgDIMAdPiwiM2bJYwYYePZSalUo0fbsGePDt99F5h7\nVgUihkEAevddA2rXVpCSwstaUuk6dJDRrJmMt97itQ4CBcMgwJw7J2DNGgnDh9ug5w4jdAOC4Bwd\nbN4sYf9+biYCAV/lADN3rgFms/OiJkQ3k5TkQL16Ct5+m6ODQMAwCCAXLghYskSPxx6zefSaBeSb\nJAkYOdKG9eslnDjBySV/xzAIIO+8Y4AkAcOH27RuCvmIgQPtqF5dxZw5HB34O4ZBgMjOFrB4sR7D\nh9tQrZrWrSFfYTIBTz5pw4oVeo4O/BzDIEDMm6eHIACPPcZRAVXM0KF21Kih4tVXjVo3hdyIYRAA\nLl4U8P77BgwdakNkpNatIV8TFASMHWvD2rUSDh7kJsNf8ZUNAG++aYBO55wMJKqMhx6yo149FTNm\ncO7AXzEM/FxWloAPP9Rj9GiOCqjy9HpgwgQrPvtMj507udnwR3xV/dy//21EjRoqHn2UowKqmt69\nHbj9dhmvvGLkOYv8EMPAj+3bJ2LtWj3Gj7fBbNa6NeTrRBF4/nkrtm2TsGULz1nkbxgGfmzaNCMa\nNZIxYACPNibX6NpVRsuWHB34I4aBn/r6ax2++UbCCy/YIFXq4qZE1xME4IUXrNi7V4f167li+ROG\ngR+SZeDFF41o3VrGAw/wzKTkWu3ayejSxYEpU4woKNC6NeQqDAM/tHSpHgcO6DB1ahGvV0Bu8eKL\nRbhwQeBJ7PwIw8DPXL4MTJ9uxKBBdrRqxcsLknvExqoYMcKGOXMMOHWKnzj8AcPAz0yfboQsO+u6\nRO709NM2REaqmDKFp6nwBwwDP7Jvn4glS/QYP96KqCju6kHuFRICTJ5sxaef6rF5M3c19XUMAz/h\ncABjx5rQtKmCoUO5Kyl5Ru/eDrRt68C4cSZYLFq3hqqCYeAnFizQ46efRLz+ehF3JSWPEQTg9ded\nk8kzZrBc5MsYBn7gyBER06cb8dhjdiQmctKYPCs2VsX48VYsWKDH7t3cpPgqvnI+zmYDRo40oX59\nhZPGpJknnrAjPl7BmDEm2HgaLJ/EMPBxr79uwOHDIt55pwgmk9atoUAlScCbbxbh2DERr77KYw98\nEcPAh+3cKWL2bAPGjbMhPp7lIdJWXJyCf/7ThrffNuDbb7l3ka9hGPgoiwUYNSoId9yh4KmnOC4n\n7zBypA1t28oYOdKEy5e1bg1VBMPAR02ZYsT58wLmzi3k3kPkNXQ6YM6cIhQVCRgzxsQzm/oQhoEP\nWr9ehw8/NGDqVCtiY/luI+9Su7aKN98swoYNeixerNe6OVRODAMfc+AAMHKkEcnJdjz8MA8uI+/U\nvbsDaWn2dwIdAAALv0lEQVQ2/OtfRuzYwc2ML+Cr5EOys4GkJKB+fQWzZvGMpOTdpk2zIjFRRlpa\nEH77jSurt2MY+AibDXjkERMKCoDly60ICdG6RUQ3ZzAA779fBIMBePhhIwoLtW4R3QzDwAeoKjB+\nvBG7d4tYtw6oU4fzBOQboqJUfPhhIQ4fFvHYY+CEshdjGPiAefP0WL7cgNmzbWjTRuvWEFVM8+YK\n3n7bio8+AqZO5YSyt+JOiV5u5UoJU6YY8dRTVgwY4ADAk4GR73nwQRkWC/DMMwaEhakYPZrHxngb\nhoEXS0+XMHq0Campdrzwgg0cyJEve/pp4MwZG156yYiICBX/+Af3hvMmDAMvtWGDhBEjTOjd24HX\nXrNCZA6QH3j+eTsuXwaee86IkBAVvXs7tG4S/Q/DwAtt3KjD8OEmdOvmwNtvF0HH07yQnxAE4JVX\nrMjPFzBihAlWaxEGDmQgeAOGgZdZvVrCU0+Z0LWrA+++ywvVkP8RReCtt4pgNBoxenQQLJYiDBvG\nkpHWuKnxEqoKvPWWAS+/bMTAgXa88QaDgPyXKAIzZ1phNgPPP2/CyZMipkyxchSsIW5uvEBhITBm\njAmffKLH2LFWjBtn4xwB+T1BcB6l7LwwkxGnTgl4550iHlCpEW5yNHb+vIDevc3IzJSwYEEhJkxg\nEFBgGTbMjo8+KsTWrRKSk804dYqnrtACNzsa2rxZh44dzTh7VkBGRgF69eJEGgWmTp1k/Pe/BcjL\nE9CpUzAyM1m08DSGgQby8oCxY40YMMCMJk0UbNpUgIQEXqmMAluzZgq++MKCtm0dGDIkCE8/bUJu\nrtatChwMAw/bskWHe+8Nxpo1erz6ahFWrSpEdDRP2EIEAOHhwKJFRZg1qxDr10u4995gbN7MWWVP\nYBh4yPnzAsaMMaJfPzMaNFDwzTcWDBli52moif5CEICHHnJgyxYLYmMVDBhgxpAhJpw8yTeLOzEM\n3Cw/H5gxw4C77gpGZqYeM2Y4RwMxMRwNEN1M3boqVq8uxLx5hdi7V4f27YPxyisG5Odr3TL/xDBw\nE4sFWLhQjzvvDMacOQYMG2bDjh35SEuzc28honISBKBPHwe2bbNg1Cgb5s0zoGXLEMycacCVK1q3\nzr9ws+Ri588LmD7dgMTEEEycaETHjjK2b7dg4kQbwsO1bh2RbwoOBiZMsGH7dgv69LHjrbec77Fp\n0ww4c4blI1dgGLiA3Q589pkOQ4ea0LJlMBYsMKB/fzt++MGCt98u4sVoiFzklltUTJ9uxa5dFqSl\n2bBokQEtWwZj0KAgfPqpBDvPalFp3Jm3kux2YMcOHTIzJXzyiYRLl0TExcn417+sGDTIzlEAkRtF\nRamYONGGMWNsSE/X46OP9EhLC0KNGgq6dXOgZ08H2rWToee1dMqNYVABZ88K+OorCV9+qcOWLRLy\n8wVERSno39+B/v3taNaMxwoQeVJICJCaakdqqh0HD4pYs0bC+vV6LF1qQHi4irvvduCee2S0by/j\nttsU7r13EwyDG8jPBw4e1OHAARG7d+uwY4cOp0+LEEUVLVsqGDXKho4dHYiPVzghTOQFmjVT0KyZ\nDRMn2nDggIgNGyRs3arDxIlGOBwCatVS0KaNjBYtZCQkKIiPlxEaqnWrvUdAh4EsOz/tnzwpIitL\nxMmTArKyRBw6pENWlnMLr9eraNZMwQMPOHDnnTLat3cgIkLjhhPRDQkCEB+vID7ehvHjnR/sduzQ\n4ZtvJOzcqcNnn0koLBQgCCrq1lXRsKGCW29VEBuroGFD5/e//U0NuLMG+83TVRSgoADIyxOQmysg\nL8/589WrAi5e/POXiAsXnD9fuCDAbneOG0VRRZ06KurVU9CliwPNmsmIi1PQqJECg0HjJ0dElRYS\nAnToIKNDBxkA4HAAR4+K2LdPxM8/63D8uICvv9Zh8WL9te2BIKioUUNFdLSKWrVUREcriI5WUbOm\nivBw51doKBAeriIszPl7cDB8ugwlqKpaoV1dVFXFqVO/QyjlWZ8+rcPmzSYUFtohyypk2bmRlmXn\n+fpL/ixc+7n4dn/+u6oCdruAoiLnZK3VKsBqBWw258/O738sdx6IUvorYTA4X9iaNZ1f1as7v0dF\nqYiJcQZAnTqq10826XQiwsKCkJtbCFnm/IQrFRRY0LTpbTh06CjM5mCtm+N3fGHdlWXgzBkBJ06I\nOH9ewPnzzg+MxV/nzzs/RMpy6dsZQVARFAQYjYDR6PzZYABMJvXaMpMJkCRApyv+rkKnc17foXj5\nH19qid8F4Y8vUfzz7wJ69tQjJqb0vo2ICEZoaGip2+wS7a9oGOTm5iKcu8oQEfmMq1evIiws7Ka3\ncenIwBfS35exf92HIwP34rrrPmX1bXlHBhWeMxAEAcHBpU/BS5KIsLBgyLIODgdfcFdj/7qf2RwM\ns5mX2nI1rrvuU1bfhoWVb5cp7hRJREQVLxMR+aPiubDy1FaJ/BHDgAjOubC8vLxy1VaJ/BHDgIiI\nOGdAREQMAyIiAsOAiIjAMCAiIjAMiIgILggDh8OB5557DjqdDhs3bizxN1VV8cILL+DWW29F9erV\n0b17d2RlZVX1IQPShx9+CJ1OB7PZDLPZjKCgIJjNZuzatUvrphHdkCiK19bV4u9PP/201s3yWZ9/\n/jlq1aqFhx566Lq/ffXVV7jrrrsQHh6O+Ph4LF++vEL3XaUwKCgoQLt27ZCTk1Pq3+fMmYMVK1Zg\nw4YNOH36NBo2bIjevXtX5SED2r333ouCggIUFBSgsLAQBQUFaNWqldbNIrohQRDwyy+/lFhnZ8+e\nrXWzfNJrr72GZ555Bo0aNbrub+fOnUNKSgqefPJJXLx4EbNmzcLw4cOxZ8+ect9/lcIgPz8fw4YN\nw/vvv4/SDldYsGABnn32WTRq1AjBwcF45ZVXcOjQIezYsaMqD0tEPkJV1VK3DVRxQUFB2LFjB269\n9dbr/rZs2TI0btwYjzzyCAwGAzp27Ijk5GQsXLiw3PdfpTCIiorC8OHDS/1bUVERDh06hDvuuOPa\nspCQENx2223YuXNnVR42YJ0+fRpdunRBZGQkGjZsiGXLlmndJKIyTZgwAfXq1UNkZCQef/xxWCwW\nrZvkk0aNGoXQG1ync/fu3UhMTCyxLDExsULbWrdNIOfk5EBVVUT85RqRkZGRyM7OdtfD+q2aNWui\ncePGmDlzJs6fP4+XX34ZaWlp+Prrr7VuGtENtWnTBl26dMGxY8ewfft2fP/99xg5cqTWzfI7ly5d\nqvK21u2XveQQ0TW6d++O7t27X/t9wIABWLt2LRYtWoT77rtPu4YR3cS2bduu/dy4cWPMmDEDycnJ\neO+996D39ksL+piqbmsrNDL46KOPEBQUdG2vgJuJjIyEKIq4dOlSieWXLl1CVFRUxVsaYMrT1/Xr\n18fZs2c93DKiyqtfvz5kWcaFCxe0bopfqVmzZpW3tRUKg8GDB6OwsPDaXgE3YzQaERcXh927d19b\nduXKFRw7dgx33XVXRR42IP21r+fPn49Vq1aVuM3hw4cRGxurUQuJbm7v3r147rnnSiw7dOgQjEYj\nateurVGr/FOrVq1KbGsBYOfOnRXa1rr1oLMRI0Zg9uzZ+Pnnn5GXl4cJEyagZcuW1010UNmsViue\neuop7N69Gw6HAx9//DE2bNiAESNGaN00olJFRUVhwYIFePXVV2Gz2fDLL79g0qRJePzxx3macBdL\nTU3FyZMn8cEHH8BqtSIzMxMbNmzA448/Xv47Uatg6dKlqslkUoOCglRRFFWj0agGBQWpjz322LXb\nTJkyRY2OjlaDg4PVnj17qr/99ltVHjKgvfzyy2qDBg3UoKAgtWnTpmpmZqbWTSK6qa1bt6p33323\nGhoaqtasWVMdN26carVatW6WTyre1kqSpEqSdO33Ylu3blUTEhJUk8mkNmnSRF23bl2F7p/XMyAi\nIp6biIiIGAZERASGARERgWFARERgGBARERgGREQEhgEREYFhQEREYBgQEREYBkREBIYBEREB+H9n\nW/Rx8Pr+NgAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# σ=2.4の正規分布\n", "T = RealDistribution('gaussian', 2.4, seed=100)\n", "plot(lambda x : T.distribution_function(x), [x, -10, 10], figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

GLMMから求められた尤度

\n", "\t

\n", "\t\t個体毎の尤度$L_i$は、以下の様に$r_i$を積分して求めることができます。\n", "$$\n", "\t\tL_i = \\int_{-\\infty}^{\\infty} p(y_i | \\beta_1, \\beta_2, r_i) p(r_i | s) dr_i\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tGLMMの結果からqを求める関数_qを定義し、上記の積分をSageの数値積分を使って計算する関数_Liを以下の様に\n", "\t\t作りました。rの積分範囲は、sの分布から-10から10としました。\n", "\t

\n", "\t

\n", "\t\tチェックのために、久保本図7.8にあるx=4, r={-2.20, -0.60, 1.00, 2.60}のqの値と、\n", "\t\ty={0, 4, 7}の尤度$L_i$を出力してみました。\n", "\t

\n", "\t

\n", "\t\t最後に、Liに20（個体数）を掛けた分布をx=4のヒストグラムに重ね合わせてみました。\n", "\t\tきれいに、図7.10(B)の通り、生存種子数を表現できているのが分かります。\n", "\t

\n", "\t

\n", "\t\tこのような積分を含む処理でも数式処理システムSageを使えば簡単に分析できるのが、\n", "\t\tご理解頂けたと思います。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 回帰分析の結果からqを求める\n", "def _q(x, r):\n", " return _logistic(-4.190 + 1.005*x + r)\n", "# 指定されたyの尤度を求める(上記分布から積分範囲は-10から10とした)\n", "def _Li(y):\n", " return numerical_integral(lambda r: _p(_q(4, r), y, 8)*T.distribution_function(r), -10, 10)[0]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.0854891394348065, 0.316479106263684, 0.696354929823834, 0.919086532784535)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "_q(4, -2.2), _q(4, -0.6), _q(4, 1.0), _q(4, 2.6)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.18814448426796998, 0.07913801041038665, 0.10840844640974827)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "_Li(0), _Li(4), _Li(7)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEDCAYAAADayhiNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl8E3X+x/HXTI42aZseUKHlFBEExAsRKIgCKiKHbhWp\nqNygHCI37IpKURBUPBCERUVWEVFA8Nh1lRWUe/H4gQtyeAAtoG0ppUfSNk0yvz8KgdoCPZJOUj7P\nx4NH6WQy827S5p3vzGRG0TRNQwghhABUvQMIIYQIHFIKQgghvKQUhBBCeEkpCCGE8JJSEEII4SWl\nIIQQwktKQQghhJeUghBCCC8pBSGEEF5SCkIIIbyqrRTGjx+PqkoHCSFEIKuWV+ldu3bx7rvvoihK\nue+jaRo5OTnIqZmEEKL6KP4+IZ6maSQkJNCnTx+mT5+O2+0ucXtGRm6Z97Pbc7n88nocOnSMsLAI\nf0asElVViIkJ4+RJOx5PYBZYMGQEyelLwZARJKcvXSxjbGz5Xkf9PlJYvHgxFouF/v37V+h+Z0YV\nFRld6EFVFRRFQVUDN+fx4ypbtijY7YGbEYLjsYTgyBkMGUFy+pKvMhp9lKdMaWlpzJgxg02bNvlz\nNeIC1q83MGSIhcJCaNTIwmefOahTJzDf6Qgh9OfXUpg4cSJDhw6lefPmHDlypMx5VLXsZjMYVO9X\nozFwd1CfmzMQvfRSCIWFxY/vkSMqK1aYmTy5SOdUZQv0x/KMYMgZDBlBcvqSrzL6rRS++uortm3b\nxhtvvAFw3h3GMTFhZW4iMhiK9z3YbBZstjB/xfQZm82id4QyhYeX/D462kx0tFmfMOUUqI/lnwVD\nzmDICJLTl6qa0W+l8N5775Genk7Dhg0B8Hg8aJrGZZddxoIFC7j//vsBOHnSXuZIwW7PByAnJx+3\n2+CvmFVmMKjYbJbTOT16xyllxgyVfv1CSU9X6NDBTVJSAVlZeqcqW6A/lmcEQ85gyAiS05culjE6\nunxvrv1WCi+//DLPPvus9/vU1FQ6dOjA7t27iY6O9k73eLQy95Sf+aHcbg8uV2A+CecK1JytWnnY\ns8eBooRhMBTgcnlwufROdWGB+lj+WTDkDIaMIDl9qaoZ/VYKkZGRREZGer8vKipCURTi4uL8tUpx\nHkYjREcTsCMEIUTgqLa9Jo0aNSr1GQUhhBCBJXB3pQshhKh2UgpCCCG8pBSEEEJ4SSkIIYTwklIQ\nQgjhJaUghBDCS0pBCCGEl5SCEEIILykFIYQQXgFbCqZPP9Y7ghBCXHICthSMe/cAoORk65xECCEu\nHQFbCgVDhgEQsmypzkmEEOLSEbCloF1WB4CQpW+i5OXqnEYIIS4NAVsKZyh5uYQufVPvGEIIcUkI\n+FJw9u2HddF8sNv1jiKEEDVewJdCweixKNnZWN55W+8oQghR4wV8KXjqN6CgbxLWBa9Afr7ecYQQ\nokYL+FIAcDw+ESXzBKHv/UPvKEIIUaMFRSl4mlxBYWJfrK+9AoWFescRQogaKyhKAcAxfjLqH78T\n+v5yvaMIIUSNFTSl4L6yGYV3/wXr/JfA6dQ7jhBC1EgVLoXdu3dz2223ERUVRVxcHElJSaSlpZWa\nLzk5GaPRiNVqxWq1YrFYsFqtZGRkVDqsY/wUDEdTCV21stLLEEIIcX4VKgWn00n37t3p2rUrGRkZ\n7Nmzh7S0NEaNGlXm/AMGDMDhcOBwOMjPz8fhcBAbG1vpsO4WLSns2QfrKy+Cy1Xp5QghhChbhUrB\n4XAwe/Zspk2bhslkolatWiQmJrJnzx5/5SudYcJkDEcOE7Lmw2pbpxBCXCoqVApRUVEMGTIEVS2+\n24EDB1i2bBlJSUllzr979246duxIZGQkrVu3Zv369VUO7Gp9LYXdexSPFtzuKi9PCCHEWZXa0ZyS\nkkJISAitWrWiXbt2zJgxo9Q89evXp2nTpixfvpy0tDSGDh1Kr169+Pnnn6uaGceEKRh//YWQjz+q\n8rKEEEKcpWiaplX2zr/++isjRoygbt26vPfeexedv3379nTv3p3k5GTvtMzMPFRVKTWv3Z5HgwZ1\nSU39g7Cw8FK3h/f9C+rRVHK27gRVv4OoDAYVm81CTk4+brdHtxwXEgwZQXL6UjBkBMnpSxfLGB0d\nVq7lVKkUAHbs2EFCQgIZGRnUqlXrgvMmJSURERHBG2+84Z2maRqKUroUcnJyiIyMJDs7G5vNVnph\n27ZBx46wahXcd19VfgQhhBCnGSsy88aNGxk5ciT79+/3TlMUBUVRMJvNJeadNWsWCQkJdOnSxTtt\n3759pfY/nDxpP89Iofg8R8WtZygdpsW1hN9yK0ryTHK73gllFEt1qAnvIAKF5PSdYMgIktOXfDVS\nqFAptGnThpycHKZOncqMGTPIy8sjOTmZzp07ExERwVVXXcXSpUtJSEggMzOT0aNHs27dOho1asSC\nBQv49ddfGThwYIllejwaHk/pwcqZH8rt9uBylf0k2MdPIeqeu1D/+U+cd95VkR/F5y6UM1AEQ0aQ\nnL4UDBlBcvpSVTNWaGO8zWZj/fr17Ny5k9jYWFq3bk1UVBQrVqwA4OeffyYvLw+AOXPm0KNHD7p1\n60ZMTAwffPABGzZsID4+vtJh/6wooRPODh2xzpsLVdsKJoQQggqOFABatWrFxo0by7zNfc4homaz\nmXnz5jFv3rzKpysHx4QpRPW9G/OG9Ti73eHXdQkhRE0XNOc+Op+izrdS1KYt1hdltCCEEFUV9KWA\nouCYNBXT999i2vS13mmEECKoBX8pAM6ut1N03fXF+xaEEEJUWo0oBRQFx4SpmHdsw7Rti95phBAi\naNWMUgCc3XvgatUa67zn9Y4ihBBBq8aUAoqCfcIUzJu/xrjzv3qnEUKIoFRzSgFw9uyN66oWhL0k\n+xaEEKIyalQpoKo4xk/GvOE/GH/4Tu80QggRdGpWKQCFff6Cq+mVWF9+Qe8oQggRdGpcKWAw4Hh8\nIiFffI7xf7v1TiOEEEGl5pUCUHjv/bgbNcb6kowWhBCiImpkKWA04hg3iZB/foJh3096pxFCiKBR\nM0sBKOibhLtBQ6wvy+cWhBCivGpsKWA243hsPCEfr8Xw80G90wghRFCouaUAFDzwEJ66cXIkkhBC\nlFONLgVCQnA8No6Qj1ah/var3mmEECLg1exSAAoeHIindizWV/17sR8hhKgJanwpYLGQP/pxQlet\nRD1yWO80QggR0Gp+KQD5AwajRUZinf+y3lGEEMI/NA327q3yYi6JUiAsDMfIsYSuXI567KjeaYQQ\nwudC/v46XH01SuaJKi3n0igFoGDIMLTwcKwLXtE7ihBC+JThwH4syU/B2LFotWpXaVl+LYXdu3dz\n2223ERUVRVxcHElJSaSlpflzleelhUeQ/8hoQpf/AzXtD10yCCGEzzmdRIwajqdhI5gzp8qL81sp\nOJ1OunfvTteuXcnIyGDPnj2kpaUxatQof63yovKHPYIWEoplwau6ZRBCCF+yvjgH47692JcsBYul\nysvzWyk4HA5mz57NtGnTMJlM1KpVi8TERPbs2eOvVV6UZoskf/ijWN5ZipKRoVsOIYTwBeN/d2Cd\n/xKOKX/Dfe11Plmm30ohKiqKIUOGoKrFqzhw4ADLli0jKSnpovf96SeVxx4LAeC33xSf5sofMRJN\nNWBd9JpPlyuEqJm+/trAgAEhPPooZGbqneYsJS8X25gRuNq0JW/UOObPN3HPPTB/vglNq/xyjb6L\nWLaUlBSuvPJK3G43I0aMYMaMGSVuV1UFVT37wp+TA/fdZ+XECTcADz0Uyvffq4SE+ChQbG0KRzyK\nZckinI+Pq/JOGYNBLfE1EAVDRpCcvhQMGSHwc+7fr/DQQxaczuLXqL17Q/nkkwKdUxWzPv0E6okM\n8j76hL+/ZWXGDDMAH39sxmTSGDnSVbkFa9Xkl19+0bp27ar179+/xHSPx1Pi+//7P00rPuA2WwM0\nyNZ++83HYdLTNc1q1bQnnvDxgoUQNcny5Wdej4r/hYbqnei0jz8uDvTGG5qmaVrfviVzJiVVftF+\nHymcccUVVzBr1iwSEhKYP38+tWrVAuDkSXuJkUJMDMTFWfj99+LvGzb0YLXaycryYRijFcvgYYTM\nf43soY+iRUVXelEGg4rNZiEnJx+32+PDkL4TDBlBcvpSMGSEwM/ZrJlCaKiFgoLi16iEBDdZWfqO\nFJSMdGxDh+K68y7siUmQZadNGyOrVp3dnHLjjYVkZZUcKURHh5Vr+X4rhY0bNzJy5Ej279/vnaYo\nCoqiYDabvdM8Hg2P5+wGMIsF1q518MorRaxcCStXFqCqBlyVHAmdT96jjxHy1hJMi17HMfmvVV6e\n2+3B5Qq8X+pzBUNGkJy+FAwZIXBzNm4Mq1bls2KFifr1TYweXaBvTk3D9vgYAHJenI/mLt6gMmSI\nE5MJdu0K4frrC3noIWelXzP9tiGvTZs25OTkMHXqVPLz88nIyCA5OZnOnTsTERFxwfs2aaIxZ44T\ngPj4KuwxuQCtTh3yHx6EZckilNwcv6xDCBH82rVzs3Chk+efB5tN3yyhK94l5N//Infea2iXXVbi\ntsGDXbzzDgwaVLV30H4rBZvNxvr169m5cyexsbG0bt2a6OhoVqxY4a9VVlj+mHEo+Q4sby3RO4oQ\nQlyQeug3wp+YSv6DA3D26Om39fh1n0KrVq3YuHGjP1dRJZ64eAr6P4xl8QIcwx6F8HC9IwkhRGlu\nN7Yxj+CpHYv9mef8uqrAPA6sGjnGTkDJzcWy7C29owghRJksC17B+P235CxcghZ+4c3vVXXJl4Kn\nfgMK+vXH+vp8cDj0jiOEECUY/7ebsLmzyH9sPK527f2+vku+FOD0aCHrJJbly/SOIoQQZ+XnEzFq\nOK6rWmL3wVGS5SGlAHgaX07hff2wvPYKFATGpxWFECJsdjKGw4fIff0NOOdQfn+SUjjNMW4iakY6\noSve1TuKEEJg2vQ11r+/jn36DNxXtai29UopnOa+4koK77kX62svg9OpdxwhxCVMOZVFxNiROG++\nhfzhI6t13VIK53CMn4x6/BihHwTOZymEEJee8GmTUPLyyJ2/CNTqfZmWUjiHu/lVFPa+B+ur86Co\nSO84QohLUMja1YR+tIq8ufPw1Ktf7euXUvgTx/jJGFKOELLmQ72jCCEuMervxwmfMoGCexIpTOyr\nTwZd1hrA3K2uprBHL6wvv4DPz8InhBDn4/EQMXYkmsVC3tyXQPHtBcbKS0qhDI4JkzEe+o2QdWv0\njiKEuESELl2C+ZuN5M5fhBYdo1sOKYUyuK69nsLbu2N95UVwu/WOI4So4QwHDxA+8ykcwx6h6Nau\numaRUjgPx4QpGA8eIOSzj/WOIoSoyZxOIkYNx92gIfbpyXqnkVI4H1ebtjhv6YL1pRfAE3gX/xBC\n1AzWl+Zi/GkPuQuXgNWqdxwphQuxT5yGcd9ezJ//U+8oQogayPjtf7G+Mg/HpGm4rrtB7ziAlMIF\nudp3wNnxZqwvPV98PWwhhPCVvDxso0fgur4NjrET9E7jJaVwEY6JUzH9bzfm9f/WO4oQogYJf/oJ\n1PQ0chYuAaNfr3dWIVIKF1HU8WaKbmovowUhhM+Yv/wcy7tvkzfzOTxNrtA7TglSChejKNgnTsX0\nw/eYvt6gdxohRJBTTpwgYtwYCm/vTsHDg/SOU4qUQjkU3dqVohvaEDZvrowWhBCVp2lETBwLmofc\nlxbo9qnlC6lwKaSkpJCYmEjt2rWJi4tj8ODB5OTklJovOTkZo9GI1WrFarVisViwWq1kZGT4JHi1\nUpTifQs7d2DaulnvNEKIIBWy8j1CPv+M3Bfno9Wpo3ecMlW4FHr37k1MTAypqal8//337N27l0mT\nJpU574ABA3A4HDgcDvLz83E4HMTGxlY5tB6ct3Wn6JrrsM6bq3cUIUQQUg8fIvxvU8h/4CGcPXvr\nHee8KlQK2dnZtG3blueeew6LxUJ8fDwDBw5k06ZN/soXOBQFx/jJmLduxrRjm95phBDBxO3GNuYR\ntFq1sD87R+80F1ShUoiMjOTNN98s8W4/JSWFevXqlTn/7t276dixI5GRkbRu3Zr169dXLa3OnD16\n4mrRSkYLQogKsSx8FeO3/yV3wd/RImx6x7mgKu1o/u6771iwYAHTp08vdVv9+vVp2rQpy5cvJy0t\njaFDh9KrVy9+/vnnqqxSX6qKY8JkzN9sxPj9t3qnEUIEAeP/dhM2dxb5Y8ZR1D5B7zgXpWha5Q6n\n2bp1K3369GHmzJmMHj26XPdp37493bt3Jzn57EmfMjPzUNXSe+Dt9jwaNKhLauofhIWFVyaif7jd\n2DrehLvx5dhXrsZgULHZLOTk5ON2B+Y5koIhI0hOXwqGjHAJ5CwowNalE5rJRO76ryEkRLeM0dFh\n5VpOpT5G9+mnn/Lwww+zcOFCHnzwwXLfr3Hjxhw/frzEtJiYMJQyDssyGIpPWW2zWbDZyvfDVJun\nnsTw0EOYf9sPbdoAxTkDXTBkBMnpS8GQEWpwzglPwqHf4Pvvia5bPddIqOpjWeFS2LZtG4MGDWLN\nmjV069btvPPNmjWLhIQEunTp4p22b98+kpKSSsx38qT9PCOFfIDTrWeoaEz/uqMXtiua4n56BgUr\nPgz4dzo1/t1YNQuGnMGQEWp2TuM3G4l4+WUczzxHYb3LIcuua0a/jBTcbjfDhw9n7ty5ZRZCixYt\neOutt0hISCAzM5PRo0ezbt06GjVqxIIFC/j1118ZOHBgift4PBoeT+ktWGd+KLfbg8sVaL8sKvbH\nJ2IbO5KCH3+ETu0CNGdJwZARJKcvBUNGqHk5lexT2EY/irNTZ+zDR0I1/mxVfSwrtKN5+/bt7N+/\nn7Fjx3o/jHbma0pKCgcPHiQvLw+AOXPm0KNHD7p160ZMTAwffPABGzZsID4+vtJhA0nhvffjbtiY\n0BflSCQhREnh0yah5OaSO38RqMF14ogKjRQ6deqE+wKXpzz3NrPZzLx585g3b17l0wUykwnH4xMI\nn/Q47NgBzVvrnUgIEQBC1q0hdM2H5Cxcgqd+A73jVFhwVViAKejXH/f1N0DXrpjWrNI7jhBCZ+rv\nxwmfMp6CPn+h8L5+esepFCmFqjCbyf3kc0hMJHz4YMJmPgUXGEkJIWowj4eIsSPRQkLJe/6lgDzZ\nXXkEzpUdgpXVCu++i+Oqq7E8/QTGvf8j5+9L0aKi9U4mhKhGoW+/gfmbjZxa+RFaTC2941SajBR8\nQVEoHP0Y2R+sxbjrB6LvuBXD/n16pxJCVBPDzwcJT36S/CHDKep6m95xqkRKwYeKbulC1hdfo1ms\nRN/ZFfNnn+gdSQjhb0VFRIwajrt+A/KeekbvNFUmpeBjnsaXk/Wv/1B42x1EDnkI65xnwRP4x18L\nISrHOm8uxr3/I/f1N4o3Jwc52afgD2Fh5L6xDNc11xI2Kxnjnh/Jff0NNFuk3smEED5k/G4n1lde\nxDH5r7iuu0HvOD4hIwV/URTyx04gZ8UqTDu2E3VnVwy/BPEZYoUQJeXlETF6BK7rb8Dx+ES90/iM\nlIKfObvdwakvN4KqEtW9C+YvP9c7khDCB8JnTMeQ9ge5C5eAseZsdJFSqAbuJk059e8NFHXqjO3h\nJKwvPS/7GYQIYub1/8byzlLykmfjbtJU7zg+JaVQTbTwCHLeXo5j8l8Jm/MstqEDUPJy9Y4lhKgg\n5cQJIsaNofC2OygYMFjvOD4npVCdVBXHpGlkv7MS0zcbibrrNtRDv+mdSghRXppGxMSx4HaR+/LC\noP3U8oVIKejAeeddnPr3BnA6ib7jVkwb/qN3JCFEOZjff4+Qzz8j98X5aHXq6B3HL6QUdOJu1pxT\nX2ykqO1NRPa/D8uCV6FyV0YVQlSHQ4ewTptMQdKDOHv10TuN30gp6EiLjCLn3Q+KT8E980kiHh0C\nDofesYQQf+Z2w8CBeKKjyZtVs6+hIqWgN4MBx1+fIvutdwj54t9E9boDNeWI3qmEEOcIfWUebNmC\nY9EStAib3nH8SkohQDh730PWv/6DmpND9B23YNr8jd6RhLjkqYcPYXvofiyzZsK0abgSOukdye+k\nFAKIu2UrstZ/jav1tUTefw+WJa/LfgYh9JCfj/X52cTcfBPGvXvIe/tdmDVL71TVQkohwGjRMWS/\nv4b8R0YTPn0aEY89Cvn5escS4tKgaZj//S9ibm6H9dV55D86hpNbvqXo7r/UyMNPy1JzPptdkxiN\n2Gc8i6v1NUSMH4Ph4H5y3n4PT736eicTosZSf/uV8OlTCfnPlzi7dCP7gzW4r7hS71jVTkYKAazw\n3vs59dmXqCdOEH37LRh3bNc7khA1j8OBdc4zxHRuh3H/PrLffo/slR9dkoUAfi6FlJQUEhMTqV27\nNnFxcQwePJicnBx/rrLGcV1zHVlffoOr+VVEJfYkdNlbsp9BCF/QNMz/+oyYm2/CuuBVHKPHcnLL\ntzh79r5kNhWVxa+l0Lt3b2JiYkhNTeX7779n7969TJo0yZ+rrJG02rXJ/nAd+YOHETFlPOGTHofC\nQr1jCRG0DL/9QuQD9xI5qD+uZs3J2rQDx1+fqhEXyakqv5VCdnY2bdu25bnnnsNisRAfH8/AgQPZ\ntGmTv1ZZs5lM2Gc9T878RYR++D5Rf+mJmvaH3qmECC52O9bZM4nu3B7DLz+T/Y/3yVmxusad6bQq\n/FYKkZGRvPnmm8TGxnqnpaSkUK9ePX+t8pJQmPQgpz7+HPXYUaJu64zxu516R/KZ9HQ4dkzvFKI6\nnTgBKSnVsCJNw/zpx8R0aot10Ws4HhvPyc07cfboeUlvKipLte1o/u6771iwYAHTp0+vrlXWWK4b\nbiTry2/wNGxE1D13EbriXb0jVdmrr5pp0cJK/frwxBNmveOIavD22yZatLDSqBGMHu2/59zwy89E\n3n8PkUMfxtWyFSc3/RfH1CfAYvHbOoOZomn+32u5detW+vTpw8yZMxk9enSJ2zIz81DV0k1tt+fR\noEFdUlP/ICws3N8RK81gULHZLOTk5ON2V/OFc5xOrNMmEbJsKQXDRpA/ay6YTIGVsRwyM6FZMyua\ndvb3YPt2B82bB+YO9UB/PCHwMzqd0KCBlaKis8/5v/6VT/v2Psyal4dl3vOEvP4anvh65D/3PEV3\n3lWpRQX64wkXzxgdHVa+BWl+9sknn2iRkZHa8uXLy7zd4/GUOT07O1sDtOzsbH/GqxkWL9Y0k0nT\nOnfWtLQ0vdNUWHq6phUfUnX23549eqcS/lRYqGkGQ8nn/JtvfLRwj0fTPvxQ0+rX17TQUE2bMUPT\nHA4fLbzm8+uH17Zt28agQYNYs2YN3bp1K3Oekyft5xkpFH+Kt7j1DP6MWSUB8Q7i/ocwNLyC8EEP\nwg1tyHv3fdzXXR9YGS/AaISJE03Mm1e8CWHAABfx8YVkZekc7DwC/fGE4Mj41FNGZswwo2kKd9/t\n5uqrC6r8nKsH9mOdNgnTN1/jvPMu8mfPxdP4cijwQIG90ssNhsfTVyMFv20+crvdXHPNNYwfP55h\nw4add76MjLIvSelw5NG4cTyHDx/Hag3czUdGo0p0dBhZWXZcLn1/WdTfj2Mb/CDGn/aSO28+hX2T\nAi7jhRw9asBqtXLZZYGdMxgez2DICPD77waMRit169qr9mKbl0fYS89jWbwAT7365M1+Huftd/os\nZzA8nhfLGBsbUa7l+G1H8/bt29m/fz9jx47FYrFgtVq9X1NTU/212kuaJy6eU+s+p/Cee7GNHkHY\nU38Dl0vvWOXWuLFG8+Z6pxDVqUEDjZYtq3AAkKYRsm4NMR1vxPLmYhwTpxYfVeTDQrjU+G3zUadO\nnXC73f5avDif0FByX32domuvI3z6NIx79+B4+x9Q3p1MQgQJw4H9hP9tMubN31DYoxd5zzyHp2Ej\nvWMFPTn3UU2kKBQMfYTs1Z9g/Ol/RHTtDLt3651KCJ9Q8nIJe/oJorskoB5N5dTKNeT8Y4UUgo9I\nKdRgRR1vJuvLb9CiouDGGwl7oC8ha1eDvfI73ITQjaYR8tEqoju0wbLsTRyT/0rWpv9S1PV2vZPV\nKFIKNZynQUNy/7UeXn4Z9eRJbI8MoXarpkSMGo75qy+hqEjviEJclGHfT0T+pSe2R4fiuvEmTm75\nFsf4yRASone0Gkeup3ApsFphzBhyHxyM55dfCf1oFSFrPiR09Qd4atemsM9fKEi8H1fbm+Qj/yKg\nKLk5WJ9/Dsubi3E3vpxTH6ylqEvZh7cL35CRwiXG0/hyHBOmkLXlW05+tYWC+/tj/vyfRPe6nZi2\n12CdPRPD/n16xxSXOk0jZNXK4k1F776N/a9PkvX1dimEaiClcKlSFNytr8E+41lO/t9PnFr7T5y3\ndMHy9pvEdG5HdJeOWF57BfXYUb2TikuMYe8eIu/ugW30CIraJ3By63fkj50gm4qqiZSCAFWlqOPN\n5M2bT+ae4tMJu5peSdgLs6l1fUsi7+5B6Dtvo2Sd1DupqMGUnGzCpk8l+rabUU9kcGrVx+S++Q+5\nDG01k1IQJYWE4OzRk9w3lpG59xdyXlsMISGETxlPrauvxPZwP0LWrQGHQ++koqbQNMwrVxDToQ2W\n5e9g/9vTxZuKbumid7JLkuxoFuelRdgo7Nefwn79UdLTCf14DSEfrcI2YjCesHCcd/Wi4N6+FHXu\nUnwCIyEqQMnJxvTtDljwCmFbt1JwTyL2GbPwxMs1V/Qkf8miXLTLLiN/+Ejyh49EPfSb9wimqFUr\ni49gujuRgsS+uG6UI5hE2ZRTWZh2bMe0dTOm7Vsx7vkRxeOBVq3IXfcZBQmd9Y4okFIQleC5vAmO\niVNxTJiCcc+PhKz+kJC1q7G8tQR3w8YU3HsfhYn3425+ld5RhY6UzExM27di2r4F87atGH7ag6Jp\nuOvVp6hDRwoGD8Nz881E3tAa1ykHBOiJ5i41Ugqi8hQFV+trcbW+FvtTMzFt30rIR6uwLH2TsJdf\npOjqayhM7Eth4n2ySeASoKSnY9qxFfO2LcUjgX0/AeBu2JiihI44HhlFUYeOxaejOD2aNBpVGVkG\nGCkF4RsoKZ3kAAARA0lEQVQGA0WdOlPUqTN5z72I+av1hK75kLC5zxL2zFMUdehI4b33U9irD1p0\njN5phQ+oaX9g2rYF07atmLZtxvjzQQBclzehKKETjtGPU5TQCU/9BjonFRUhpSB8LyQE5129cN7V\nCyU3B/M/PyV0zYeETx5H+LSJOLvdXlwQt99Z/GlrERTUY0eLS2D7VkzbtmD87VcAXE2vpKhDJxwT\nphSXQFy8zklFVUgpCL/SImwUJj1IYdKDKGlpZ49gGj6o+Aimnr0pSOxLUedbwei/i7eLilNTjnhL\nwLxtC4YjhwFwNb+Kolu64Pjrkzjbd0SrU0ffoMKnpBREtdHq1CF/xCjyR4xC/e3sOZiiPnwfT+1Y\nnIn3wt29UWPqQFx9GUVUJ01DPXwI8+lRgGnbFgxHiy+G5Wp5NYW3d6eoQyeKOnREq11b57DCn6QU\nhC48Ta7AMWkajolTMf64i5A1qwhdtxqWLCbyzDy1a+Nu0BB3g0Z4GjTE3aAhngYNcDdohLt+AwgP\n3Mu0BjxNw/DbL5i2bvGOBgy/H0dTFFxXX0Nhz97FJdC+A1pMLb3TimokpSD0pSi4rr0e17XXU/js\nbKLtWeT+bx/aocMYjqaipqZgSE3BtPv/UI8dRTnn8qKeWrWKi6J+cWG4GzY8XR6N8DRogBZevmvS\nXhI0DcPBA6cLoHjnsCE9DU1VcV1zLYX33EtRx04UteuAFhmld1qhIykFEThUFRo2xBVRC9dNCaVv\nd7tR0/5ATUnBkHoEQ2rK6dJIxfzvf2I4mopyzvUhPNHRJUYZ7oZnC8TTsCFahK0af7hqomngcKAU\nOOBIDiFffoV1S3ERqCdOoBkMuK67gcJ+/SlK6EjRTe1r5uMgKk1KQQQPgwFPfD088fVwte9Q+naP\n52xpHE05WxopRzB/+XlxaTidZ2ePiioujfoNSowyzmym8vs7Zo+n+AXcbke156LY7af/5cE5/1fz\n8lDsdrDnlZhHzTsz3znTHXYUTfOuwmIy4bq+DfkPDaKoQ0eK2raTzW7igqQURM2hqnji4vHExeNq\n17707R4PakY6akrJUYYh9Qjmr9ZjSE1BKSw8O7st8k+jjOL9GUrjRmCzYPw9AzWn5Iu5Yrej5OWV\nfKE+33THxS+LqoWGooWFoYWFF3+1hqGFh6OFheOOqVXytrDw07eFodoiCG8Qx6nLm+Myh/rwQRY1\nXYVL4YsvvmDgwIF07dqVFStWnHe+5ORknnnmGczm4sMMNU1DURSOHDlCbGxs5RMLUVmqiqdOXTx1\n6uJq26707R4PSkbG2VFGSvFXQ+oRzF9vKC6N/Hzv7H/eY6FZLKdfuE+/SIeffbH21I4t+QJ++oW9\neFrZL+yaNazSJxo0GlWIDoMsu5w+QlRIhX7jXnjhBZYuXUqzZs3KNf+AAQNYunRppYIJUe1UFa1O\nHVx16uBq07b07ZqGcuIE5t+PYosKI9ut4rKEnX0BNxiqP7MQPlahUrBYLOzcuZOxY8dSeM4wW4hL\ngqKgxcbijqsD0WF4sux45F24qGEqdJGdMWPGEBFR/sP8du/eTceOHYmMjKR169asX7++wgGFEEJU\nH79dea1+/fo0bdqU5cuXk5aWxtChQ+nVqxc///yzv1YphBCiivx29NHQoUMZOnSo9/tx48axcuVK\nli9fTnJysne6qiqoaulT5xoMqver0Ri4Vw09N2egCoaMIDl9KRgyguT0JV9lrNZDUhs3bszx48dL\nTIuJCUMp43zqBoMbAJvNgs0WVi35qsJms+gd4aKCISNITl8KhowgOX2pqhn9VgqzZs0iISGBLl3O\nXnx73759JCUllZjv5El7mSMFu7340L+cnHzc7sA9qsNgULHZLKdzBuZOx2DICJLTl4IhI0hOX7pY\nxujo8r259mkptGjRgrfeeouEhAQyMzMZPXo069ato1GjRixYsIBff/2VgQMHlriPx6Ph8WillnXm\nh3K7PbiC4AiPYMgZDBlBcvpSMGQEyelLVc1Y4UNSFUWh6PT5ZdauXYuiKDgcDgAOHjxIXl4eAHPm\nzEFRFLp168bJkydp1aoVGzZsID5eLsAhhBCBqkKlkH/OpznL4na7vf83m83MmzePefPmVS6ZEEKI\nahe4u9KFEEJUOykFIYQQXlIKQgghvKQUhBBCeEkpCCGE8JJSEEII4SWlIIQQwktKQQghhJeUghBC\nCC8pBSGEEF5SCkIIIbykFIQQQnhJKQghhPCSUhBCCOElpSCEEMJLSkEIIYSXlIIQQggvKQUhhBBe\nUgpCCCG8pBSEEEJ4SSkIIYTw8nspfPHFF9StW5f+/fv7e1VCCCGqyOjPhb/wwgssXbqUZs2a+XM1\nQgghfMSvIwWLxcLOnTu54oor/LkaIYQQPuLXkcKYMWP8uXjdff21gS++MNG6NTz4oN5pRHX4738N\nfPyxiaZNYdAgUANwr5ymwT/+YeTgQUhIMNC9u0fvSCKI+LUUarIdOwwkJVnweBQADh40M2NGgc6p\nhD/t2aNy770WnM7i53z3bjOvvhp4z/lrr5l59tkQABYtCmXZMo277nLpnEoEC91LQVUVVFUpNd1g\nUL1fjcbAezu2ZYvRWwgA33xjCMicUPKxDGSBnnPnTqO3EAC++cYYkM/5pk0l/6w3bzbSp09gjhYC\n/Tk/Ixhy+iqj7qUQExOGopRVCm4AbDYLNltYdce6qISEkt/fcINKdHTg5TyXzWbRO0K5BGrODh1A\nUYo3zwBcf70SkM9527awadPZ79u1MxEdbdIvUDkE6nP+Z8GQs6oZdS+FkyftZY4U7PZ8AHJy8nG7\nDdUd66I6d4YXXzTy6adGWrY08MQT+WRlBe67MZvNcvqxDMyMEPg5W7eG11838sEHRpo0MfDkk4H5\nnE+aBIWFZvbtM3HzzUXcd5+TrCy9U5Ut0J/zM4Ih58UylvcNjO6l4PFoeDxaqelnfii324PLFZhP\nwoABToYMcREdHUZWVuDmPCOQH8tzBXLOe+910q9fYD/nBgPMnOkkOtpEVpYzIDP+WSA/5+cKhpxV\nzahomlb6FdlHLBYLiqJQVFQEgNFoRFEUHA7HRe+bk5NDZGQk2dnZ2Gw2f0UUQghxDr+WQlVomkZu\nbi4RERFl7nMQQgjhewFbCkIIIapf4B5fJYQQotpJKQghhPCSUhBCCOElpSCEEMJLSkEIIYSXlIIQ\nQggvKQUhhBBeUgpCCCG8pBSEEEJ4SSkIIYTw0vUsqWfObySEEML/ynMuOV1LITc3l8jISD0jCCHE\nJaM8Z53W9YR4MlIQQojqU56RgpwlVQghhJfsaBZCCOElpSCEEMJLSkEIIYSXlIIQQggvKQUhhBBe\nUgpCCCG8pBSEEEJ4SSkIIYTwCshSSElJoVevXtSuXZvLL7+cadOm6R2pTF988QV169alf//+ekc5\nr5SUFBITE6lduzZxcXEMHjyYnJwcvWOVsnv3bm677TaioqKIi4sjKSmJtLQ0vWOd1/jx41HVgPzz\nQVVVLBYLVqvV+/Xxxx/XO1aZZs2aRXx8PBEREdxxxx0cOXJE70glbN682fsYnvkXGhqKwWDQO1oJ\nu3btolu3bkRHRxMfH8/DDz/MiRMnKrWsgPytTkxMpEGDBhw+fJj//Oc/rF27lldeeUXvWCW88MIL\njBs3jmbNmukd5YJ69+5NTEwMqampfP/99+zdu5dJkybpHasEp9NJ9+7d6dq1KxkZGezZs4e0tDRG\njRqld7Qy7dq1i3ffffeipwvQi6IoHDx4EIfDQX5+Pg6Hg1dffVXvWKUsXLiQFStWsGnTJn7//Xda\ntmzJyy+/rHesEm6++WbvY3jm39NPP02/fv30jubldrvp2bMnCQkJZGRksHfvXtLT0xk9enTlFqgF\nmG+//VYzmUxadna2d9rixYu1Fi1a6JiqtNdee03LycnRBg0apD3wwAN6xynTqVOntKFDh2rp6ene\naQsWLNCaN2+uY6rSsrKytLfeektzu93eafPnz9eaNWumY6qyeTwerX379trs2bM1VVX1jlMmRVG0\nI0eO6B3jopo0aaKtW7dO7xgVcuTIEa127dra0aNH9Y7ilZqaqimKou3fv987bfHixdqVV15ZqeUF\n3Ejhhx9+oHHjxiXO5HfDDTdw4MAB7Ha7jslKGjNmDBEREXrHuKDIyEjefPNNYmNjvdNSUlKoV6+e\njqlKi4qKYsiQId7NMQcOHGDZsmUkJSXpnKy0xYsXY7FYAnqTIcDUqVNp1KgRMTExPPLIIwH1twNw\n/PhxDh06RGZmJq1ataJ27dr07du30ps8qstTTz3FsGHDAupvqF69elx//fUsWbIEu91Oeno6a9as\noXfv3pVaXsCVQmZmJtHR0SWmxcTEAAT8L0yg++6771iwYAHTp0/XO0qZUlJSCAkJoVWrVrRr144Z\nM2boHamEtLQ0ZsyYwaJFi/SOckEdOnTgjjvu4JdffmH79u3s2LGj8psS/OTo0aMArF69mg0bNvDj\njz9y9OhRRowYoXOy8zt8+DBr165l/PjxekcpQVEUVq9ezbp167DZbMTFxeF2u5k9e3allhdwpQDF\np9QWvrV161a6d+/O888/T5cuXfSOU6aGDRtSWFjIgQMHOHDgAA899JDekUqYOHEiQ4cOpXnz5npH\nuaCtW7cyePBgTCYTzZs3Z+7cuaxYsYKioiK9o3md+RufOnUqderUIT4+nuTkZD755BOcTqfO6cq2\ncOFCEhMTueyyy/SOUoLT6aR3797069eP7Oxsjh07hs1mq/RoNuBKITY2lszMzBLTMjMzURSlxGYQ\nUX6ffvopPXv2ZP78+QH3jrEsV1xxBbNmzeL9998v9bugl6+++opt27bx5JNPAsH1xqVx48a43W7S\n09P1juJVt25dgBIX2WrcuDGapgVUznOtXr2aPn366B2jlK+++orDhw8ze/ZswsPDqVu3LsnJyaxd\nu5ZTp05VeHkBVwo33ngjKSkpnDx50jtt586dtGzZEqvVqmOy4LRt2zYGDRrEmjVrePDBB/WOU6aN\nGzdy1VVXlZimKAqKomA2m3VKVdJ7771Heno6DRs2JDY2ljZt2qBpGpdddhkffvih3vG8du3aVero\nsp9++omQkBDi4+N1SlVa/fr1sdls7Nq1yzvt0KFDmEymgMp5xu7du0lJSeH222/XO0opbrcbj8eD\nx+PxTisoKKj80XE+2Pntcx06dNCGDx+u5eTkaPv27dOaNGmiLVq0SO9YZQrko49cLpfWsmVL7Y03\n3tA7ygVlZ2drcXFx2pQpUzSHw6Glp6drPXr00G699Va9o3mdOnVKO3bsmPffjh07NEVRtOPHj2v5\n+fl6x/M6duyYFhERoc2dO1crLCzUDhw4oLVq1UobN26c3tFKmTBhgta0aVPtl19+0dLS0rSOHTtq\nw4YN0ztWmd5++20tNjZW7xhlyszM1GJjY7Xp06drDodDO3HihHb33XdrXbp0qdTyArIUjh07pt11\n112a1WrV4uLitJkzZ+odqZTQ0FDNYrFoRqNRMxqN3u8DyebNmzVVVTWLxeLNd+ZrSkqK3vFK2LNn\nj3brrbdqYWFhWp06dbT+/ftrx48f1zvWeR0+fDhgD0ndvHmzlpCQoEVERGixsbHa5MmTtcLCQr1j\nlVJYWKiNGTNGi4mJ0Ww2mzZkyBDNbrfrHatMzz33nNa6dWu9Y5zXDz/8oHXp0kWLiYnR4uLitAce\neED7/fffK7UsuRynEEIIr4DbpyCEEEI/UgpCCCG8pBSEEEJ4SSkIIYTwklIQQgjhJaUghBDCS0pB\nCCGEl5SCEEIILykFIYQQXlIKQgghvKQUhBBCeEkpCCGE8Pp/hFSdgIk09xcAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "prd2_plt = list_plot([_Li(x)*20 for x in range(9)], plotjoined=True, rgbcolor=\"red\")\n", "(lst_plt + prd2_plt).show(figsize=4)" ] }, { "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 }