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

入門機械学習による異常検出

\n", "\t

\n", "\t\t井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。\n", "\t

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

2章 正規分布に従うデータからの異常検出

\n", "\t

\n", "\t\t井出本の基本は、データに対するモデルを使って負の対数尤度を求め、それを異常検出関数として使うことです。\n", "\t

\n", "\t

多変数正規分布に基づく異常検知

\n", "\t

\n", "\t\tM次元、N個の観測データ$\\mathcal{D}$が以下の様な多次元正規分布に従っていると、仮定します。\n", "$$\n", "\t\\mathcal{N}(x|\\mu, \\Sigma) = \\frac{|\\Sigma|^{-1/2}}{(2 \\pi )^{M/2}} exp \\left \\{ -\\frac{1}{2}(x - \\mu)^T \\Sigma^{-1} (x - \\mu) \\right \\}\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\t対数尤度は、以下の様に計算されます。\n", "$$\n", "\tL(\\mu, \\Sigma | \\mathcal{D}) = ln \\prod_{n=1}^{N} \\mathcal{N} (x^{(n)} | \\mu, \\Sigma) = \\sum_{n=1}^N ln \\mathcal{N} (x^{(n)} | \\mu, \\Sigma)\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tこれに多次元正規分布$\\mathcal{N}$の定義を代入すると、対数尤度Lは以下の様になります。\n", "$$\n", "\tL(\\mu, \\Sigma | \\mathcal{D}) = -\\frac{M N}{2}ln(2 \\pi) - \\frac{N}{2} ln | \\Sigma | \n", "\t\t-\\frac{1}{2} \\sum_{n=1}^N (x^{(n)} - \\mu)^T \\Sigma^{-1} (x^{(n)} - \\mu)\n", "$$\t\t\n", "\t\t\n", "\t

\n", "\t

\n", "\t\tこれを$\\mu, \\Sigma$で偏微分して、尤度Lを最大にする$\\mu, \\Sigma$を求めます。\n", "\t

\n", "\t

\n", "\t\t最初にLを$\\mu$で偏微分すると、以下の様になります。\n", "$$\n", "\t\\frac{\\partial L(\\mu, \\Sigma | \\mathcal{D})}{\\partial \\mu} = \\sum_{n=1}^N \\Sigma^{-1} (x^{(n)} - \\mu)\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\t上記の式が0になる$\\hat{\\mu}$は以下の様に求まります。\n", "$$\n", "\t\\hat{\\mu} = \\frac{1}{N} \\sum_{n=1}^N x^{(n)}\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\t次に、$\\Sigma$での偏微分です。これには、以下の関係式を使います。\n", "$$\n", "\t-ln | \\Sigma | = ln | \\Sigma^{-1} |\n", "$$\t\n", "\tと、\n", "$$\n", "\t(x^{(n)} - \\mu)^T \\Sigma^{-1} (x^{(n)} - \\mu) = Tr \\left \\{ \\Sigma^{-1} (x^{(n)} - \\mu) (x^{(n)} - \\mu)^T \\right \\}\n", "$$\t\n", "\tと、以下の公式を使います。\n", "$$\n", "\t\\frac{\\partial}{\\partial A} Tr(A B) = \\frac{\\partial}{\\partial A} Tr(B A) = B^T\n", "$$\t\n", "$$\n", "\t\\frac{\\partial}{\\partial A} ln |A| = (A^{-1})^T\n", "$$\n", "\t

\n", "\t

\n", "\t\t偏微分の結果は、以下の様になります。共分散行列$\\Sigma$が対象行列であるので、$\\Sigma = \\Sigma^T$を使っています。\n", "$$\n", "\\begin{eqnarray}\n", "\t\\frac{\\partial L(\\mu, \\Sigma | \\mathcal{D})}{\\partial (\\Sigma^{-1})} & = & \\frac{N}{2}ln | \\Sigma^{-1} |\n", "\t\t- \\frac{1}{2} \\Sigma^{-1} (x^{(n)} - \\mu) (x^{(n)} - \\mu)^T \\\\\n", "\t\t& = & \\frac{N}{2} \\Sigma^T - \\frac{1}{2}\\sum_{n=1}^N (x^{(n)} - \\mu) (x^{(n)} - \\mu)^T \\\\\n", "\t\t& = & \\frac{N}{2} \\Sigma - \\frac{1}{2}\\sum_{n=1}^N (x^{(n)} - \\mu) (x^{(n)} - \\mu)^T \n", "\\end{eqnarray}\t\t\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tこの値が0になるような$\\hat{\\Sigma}$は、以下の様に求まります。\n", "$$\n", "\t\\hat{\\Sigma} = \\frac{1}{N} \\sum_{n=1}^N (x^{(n)} - \\mu) (x^{(n)} - \\mu)^T \n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\t負の対数尤度を基に、異常度$a(x')$を以下の様に定義します。これは、マハラノビス距離と呼ばれ、「ばらつきが大きい方向の変動は大目にみる」効果があるそうです。\n", "$$\n", "\ta(x') = ({x'} - \\mu)^T \\Sigma^{-1} ({x'} - \\mu)\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tしかし、実際のRでのプログラムをみると、ちょっと違っています。実行例2.5(Davisに対する異常度の計算部分)は、以下の計算をしています。\n", "$$\n", "\ta(x') = \\sum_{m=1}^M ({x'} - \\mu)^2 \\Sigma^{-1}\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\t値としては、Wikipediaにある共分散行列が対角行列であるときの以下の式に近いと思われます。\n", "$$\n", "\td(\\vec{x}, \\vec{y}) = \\sqrt{\\sum_{i=1}^P \\frac{(x_i - y_i)^2}{\\sigma_i^2}}\n", "$$\t\t\n", "\t

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

前準備

\n", "\t

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

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

データを使って異常度を計算する

\n", "\t

\n", "\t\tRのcarパッケージに入っているデータDavisを使って、多次元データの異常度を求めてみましょう。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ " [1] \"car\" \"jsonlite\" \"ggplot2\" \"stats\" \"graphics\" \"grDevices\" \"utils\" \"datasets\" \n", " [9] \"methods\" \"base\" " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 2.2 carパッケージのDavisを使って例題を試す\n", "# carパッケージのインストールには、時間が掛かりますので注意してください。\n", "r('library(car)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

データの性格を知る

\n", "\t

\n", "\t\tRのデータをpandaのデータフレームに変換して、データの性格を調べます。\n", "\t\t最初にinfoとdescribeで個数や欠損値、平均、分散などの統計情報を求めます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Int64Index: 200 entries, 0 to 199\n", "Data columns (total 5 columns):\n", "height 200 non-null int64\n", "repht 183 non-null float64\n", "repwt 183 non-null float64\n", "sex 200 non-null object\n", "weight 200 non-null int64\n", "dtypes: float64(2), int64(2), object(1)\n", "memory usage: 9.4+ KB\n" ] } ], "source": [ "# Rのデータをpandaのデータフレームに変換する\n", "davis = RDf2PandaDf(\"Davis\")\n", "# データの個数と欠損値の有無をみる\n", "davis.info()" ] }, { "cell_type": "code", "execution_count": 7, "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", "
heightrephtrepwtweight
count200.000000183.000000183.000000200.000000
mean170.020000168.49726865.62295165.800000
std12.0079379.46704813.77666915.095009
min57.000000148.00000041.00000039.000000
25%164.000000NaNNaN55.000000
50%169.500000NaNNaN63.000000
75%177.250000NaNNaN74.000000
max197.000000200.000000124.000000166.000000
\n", "
" ], "text/plain": [ " height repht repwt weight\n", "count 200.000000 183.000000 183.000000 200.000000\n", "mean 170.020000 168.497268 65.622951 65.800000\n", "std 12.007937 9.467048 13.776669 15.095009\n", "min 57.000000 148.000000 41.000000 39.000000\n", "25% 164.000000 NaN NaN 55.000000\n", "50% 169.500000 NaN NaN 63.000000\n", "75% 177.250000 NaN NaN 74.000000\n", "max 197.000000 200.000000 124.000000 166.000000" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# データの情報を取り出す\n", "davis.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t次にDavisの体重のヒストグラムと体重(weight)と身長(height)の散布図をプロットし、おおざっぱな性質を調べます。\n", "\t

\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEDCAYAAADOc0QpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEmNJREFUeJzt3X+Q3HV9x/Hn5WJjLxdjwlwwjZWUCG9lrH9YO7bRmkgp\n1oqgBlsHzKCxyjg4Q4tDR50iNWp1oGFEqKODgDEDU2x1kJQRJSNi64+q7Yw/qr5R0iSVRO7infHO\nhJjeXf/4LnqBC7e79927vU+ej3+y+9397vc138u+vt/9fL/73Z7JyUkkSWVZNN8BJEn1s9wlqUCW\nuyQVyHKXpAJZ7pJUIMtdkgq0uJknRcTFwJXAMeBdwHeAHVQbhwPA5sw81qmQkqTWzLjnHhErqQp9\nPXAe8EpgK3BDZm4AHgS2dDKkJKk1zQzLnAPcm5mHM/PhzLwU2AjsbDy+s/EcSVKXaGZYZi2wNCI+\nAzwVeDfQN2UYZhBY3Zl4kqR2NFPuPcBK4FVURX9fY9rUxyVJXaSZcn8Y+EpmTgC7I2IUOBYRSzLz\nKLAG2P9ELzA5OTnZ0+M2QJJa1HZxNlPunwdujYhrqPbg+4F7gAuB24BNjfsnTtfTw9DQaLsZ58zA\nwDJz1sic9TJnfRZCRqhytmvGA6qZuR/4F+BrwN3AZcDVwCURcT+wAtjedgJJUu2aOs89M28CbnrM\n5HPrjyNJqoPfUJWkAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtd\nkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWp\nQJa7JBXIcpekAlnuklQgy12SCmS5S1KBFs/0hIjYAPwz8F2gB/g2cC2wg2rjcADYnJnHOpizCOPj\n4+zZs7vl+dauPZ3e3t4OJJJUqhnLveGLmfnnj96JiFuAGzLz0xHxPmAL8NFOBCzJnj27ufzau+hb\nvqrpeQ4fGuT6K89n3bozOphMUmmaLfeex9zfCFzauL0TeBuWe1P6lq+if8Wa+Y4hqXDNlvtZEXEn\nsBLYCvRNGYYZBFZ3IpwkqT3NHFD9IfB3mflK4PXAzRy/UXjsXr0kaZ7NuOeemfupDqiSmbsj4ifA\n8yNiSWYeBdYA+2d6nYGBZbPNOic6mXNkpL+t+Vau7H9cLtdnvcxZr4WQcyFknI1mzpa5CFidmdsi\n4mnAqcCtwIXAbcAm4J6ZXmdoaHSWUTtvYGBZR3MOD4+1Pd/UXJ3OWRdz1suc9VkIGWF2G6Bmxtzv\nAm6PiAuAJ1EdSP0W8ImIeDOwF9jedgJJUu2aGZYZA86f5qFz648jSaqD31CVpAJZ7pJUIMtdkgpk\nuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQJa7\nJBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtdkgpkuUtS\ngSx3SSrQ4maeFBFPBr4LbAW+AOyg2jAcADZn5rGOJZQktazZPfergJ82bm8FbsjMDcCDwJZOBJMk\ntW/Gco+IAJ4F3A30ABuAnY2HdwLndCydJKktzey5bwOuoCp2gKVThmEGgdWdCCZJat8TjrlHxGbg\nK5m5t9qBf5ye6SZOZ2BgWYvR5kcnc46M9Lc138qV/Y/L5fqslznrtRByLoSMszHTAdWXA78TEa8A\n1gC/BMYiYklmHm1M29/MgoaGRmcVdC4MDCzraM7h4bG255uaq9M562LOepmzPgshI8xuA/SE5Z6Z\nr330dkS8C9gDrAcuBG4DNgH3tL10SVJHNHUqZMOjQzBXAzsi4s3AXmB77an0K5MTE+zbt/e4aSMj\n/TN+Cli79nR6e3s7GU1SF2u63DPz3VPuntuBLJrGkdEhtt1xkL7lB5qe5/ChQa6/8nzWrTujg8kk\ndbNW9tw1T/qWr6J/xZr5jiFpAfHyA5JUIMtdkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QC\nWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQJa7JBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDl\nLkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtdkgq0eL4DqH6TExPs27e3rXnXrj2d3t7emhNJmmszlntE\n/CbwceBUYAnwXuBbwA6qPf8DwObMPNa5mGrFkdEhtt1xkL7lB1qa7/ChQa6/8nzWrTujQ8kkzZVm\n9txfAXwjM/8hIp4B3At8GbgxMz8VEe8DtgAf7WBOtahv+Sr6V6yZ7xiS5smM5Z6Zn5xy9xnA/wIb\ngEsb03YCb8Nyl6Su0fSYe0R8GVhDtSd/75RhmEFgdQeySZLa1HS5Z+YLI+K5wG1Az5SHek4wy3EG\nBpa1GG1+dDLnyEh/x167LitX9te6Dvy718uc9VkIGWejmQOqzwMGM/PHmfntiOgFRiNiSWYepdqb\n3z/T6wwNjc4+bYcNDCzraM7h4bGOvXZdhofHalsHnV6fdTFnvRZCzoWQEWa3AWrmPPcXU42pExGn\nAv3ALuDCxuObgHvaTiBJql0zwzIfAW6OiC8BTwbeAvwnsCMi3gzsBbZ3LqIkqVXNnC3zCHDxNA+d\nW38cSVIdvPyAJBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJU\nIMtdkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy\n3CWpQJa7JBXIcpekAlnuklQgy12SCrS4mSdFxDXAi4Be4APAN4AdVBuHA8DmzDzWqZCSpNbMuOce\nERuBszJzPfAy4IPAVuDGzNwAPAhs6WRISVJrmhmWuR94TeP2z4ClwAbgrsa0ncA59UeTJLVrxmGZ\nzJwEjjTuvhG4G3jplGGYQWB1Z+JpLk1OTLBv396W51u79nR6e3s7kEhSu5oacweIiAuohl/OBX40\n5aGeZuYfGFjWWrJ50smcIyP9HXvtOhwZHWLbHQfpW36g6XkOHxpkx/sv4swzz5z2cf/u9TJnfRZC\nxtlo9oDqS4F3UO2xj0bEaEQsycyjwBpg/0yvMTQ0Orukc2BgYFlHcw4Pj3XstevSt3wV/SvWtDTP\n8PDYtOut0+uzLuas10LIuRAywuw2QM0cUH0KcA1wXmYeakzeBWxq3N4E3NN2AklS7ZrZc/8L4BTg\nkxHRA0wClwA3R8SlwF5ge+ciSpJa1cwB1ZuAm6Z56Nz640iS6uA3VCWpQJa7JBXIcpekAlnuklQg\ny12SCmS5S1KBmr78gDSdJ7oezchI/wm/lev1aKTOstw1K+1ej+b6K89n3bozOphMOrlZ7pq1dq5H\nI6mzLPc2jI+Ps2fP7pbna+dyupLUDsu9DXv27Obya++ib/mqlub76Y+/zylPf3aHUknSr1nubWpn\nKOLwoYc7lEaSjuepkJJUIMtdkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIK\nZLlLUoEsd0kqkOUuSQWy3CWpQJa7JBXIcpekAjX1Yx0R8RzgTuC6zPxwRDwd2EG1cTgAbM7MY52L\nKUlqxYx77hHRB3wI2DVl8lbghszcADwIbOlMPElSO5oZlnkEeBnVHvqjNgI7G7d3AufUG0uSNBsz\nlntmTmTm0cdMXjplGGYQWF17MklS2+r4geyeGl5DJ5HJiQn27dvb1rxr155Ob29vzYmk8rRb7qMR\nsaSxR78G2D/TDAMDy9pc1NxqJufISP8cJCnXkdEhtt1xkL7lB2Z+8hSHDw2y4/0XceaZZ9aeqaT/\nn91gIeRcCBlno91y3wVsAm5v/HvPTDMMDY22uai5MzCwrKmcw8Njc5CmbH3LV9G/Yk3L8w0Pj9X+\nf6nZv/t8M2d9FkJGmN0GaMZyj4jnAduA04BjEXEhcDGwPSIuBfYC29tOIEmq3Yzlnpn/BbxkmofO\nrT+OJKkOfkNVkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIKZLlLUoEsd0kq\nkOUuSQWy3CWpQHX8EpM0J9r9BSd/vUknI8sdeN1b3s7i/tX0Ll7E+P9NzPj8sUMHWbSi/l8D0hNr\n5xecDh8a5Porz2fdujM6mEzqPpY78KRlq+lZ+btM0tw41aJFD3U6kk6g3V9wkk42jrlLUoEsd0kq\nkOUuSQWy3CWpQJa7JBXIs2VUtGbOjR8Z6Wd4eOy4aZ4br4XOclfRPDdeJyvLXcXz3HidjBxzl6QC\nWe6SVCDLXZIKZLlLUoHaPqAaEdcBfwBMAH+Vmd+sLZU0j+by0sLj4+Ps2bO7qedOPWXTUzXbNz4+\nzgMPPPC4019nstDWeVvlHhEvBp6Zmesj4lnALcD6WpNJ82QuT5/cs2c3l197F33LV3V8WaqcLOu8\n3T33PwbuBMjMH0TEUyOiPzNb2xRKXWouT5/0VM25dzKs83bH3J8GDE25f7AxTZLUBer6ElNPTa8z\nL46NHmDxJE3/EtPEoYM8suipLS/nyOgwra6quZpnLpdVYr7DhwbbGqfft28vhw8Nzsmy5sp0l3Po\nJu2u84WmZ3JysuWZIuJqYH9m3tS4/yDw3Mz8Rc35JEltaHdY5vPAhQAR8TzgIYtdkrpHW3vuABHx\n98AGYBy4LDO/U2cwSVL72i53SVL38huqklQgy12SCmS5S1KBOvJjHRFxDfAioBf4APANYAfVxuQA\nsDkzj3Vi2a2KiCcD3wW2Al+gC3NGxMXAlcAx4F3Ad+iynBGxFPgEsAL4Dar1+T26JGdEPIfqW9XX\nZeaHI+Lp02VrrOvLqU4UuCkzb5nnnL9NdXmPJwG/BF6XmYPdlnPK9JcCn83MRY37XZUzIhYD24Fn\nAj8HLszMQ12Y88XA+6je82NU/z9byln7nntEbATOysz1wMuAD1K90W/MzA3Ag8CWupc7C1cBP23c\n3grc0E05I2IlVaGvB84DXkkX5gReD/wgM88GXgNcT5f83SOiD/gQsGvK5Metw8bzrgLOBl4C/HVE\ntP5ttXpzvgf4SGZupHrzX9GlOYmIJcDbgf1TntdtOd8EDGbmC4A7gD/q0pzbgDc03k9fBS5tNWcn\nhmXup3pzA/wMWEp1yuRdjWk7gXM6sNyWRUQAzwLupvpK4gaqfNA9Oc8B7s3Mw5n5cGZeCmyk+3Ie\nBE5p3F5JdXmKbvm7P0K1ozH1SmAbOX4d/gnwAuDrmTmWmY8A/w68cJ5zvgX4dOP2ENU67sacAO8E\nbqT6hAHdmfMVwG0AmfmxzPzXLs05BAw0bq+gen+1lLP2cs/Mycw80rj7RqriXDrl4/ggsLru5bZp\nG3AFv/6ueTfmXAssjYjPRMT9EXE20NdtOTPzDuC0iPgh8EWqYaSuWJ+ZOZGZRx8zebpsp3L8NZOG\nmMPM0+XMzCOZORkRi4DLgNt5/LWd5j1nRJxJ9S31T02Z3HU5qd5PfxYR90XE7RGxgu7MeQVwZ0R8\nn2qI++O0mLNjB1Qj4gKqj+Fv5fgLdXTFdWgiYjPwlcw80UU6uiInVY6VwKuANwC30p3r82Jgb2ae\nQfWx8R8f85SuyHkCJ8rWFZkbxb4D2JWZ903zlG7IeR1VIUF3r88e4PuZ+RLgv4F3nOA58+0G4ILM\nfDbVHvpl0zznCXN2pNwbB1XeAfxpZo4Co43xOIA1NMbk5tnLgQsi4qtUnzCuAsa6MOfDVBuhiczc\nDXTr+nwh8DmAxreVVwO/6MKcj3rsOnyIKt/UPaFuyXwrkJn53sb9rsoZEb8FBHBb4/20OiLuo1qn\nXZOz4SfAlxq3PwecRXfmfG5mfq1xexfwe7SYsxMHVJ8CXAOcl5mHpoTb1Li9Cbin7uW2KjNfm5kv\nyMw/BD5GdYBtF41r5tAlOamu43N2RPRExClAP92Z80dUv8xFRJxGtRG6l+7L+ajp/k9+HXh+RDwl\nIvqpDmL/2zzlA371iehoZm6dMvk/6J6cPZm5PzPPyMz1jffTgcaecdetT+CzVOPbUBVm0p05DzR+\nCAng94Ef0mLO2i8/EBFvAq4GHqD62DAJXALcDCwB9lIdBR6vdcGz0LjK5f9Qbcl30GU5G+v0L6nW\n5XuAb9JlORunQt5CNW7dC/wt1RvnE8xzzsbF7bYBp1GdWvYQcDHVKXHHZYuIVwN/Q/XzkR/KzH+a\n55yrqA64jVL9/b+XmW/twpyvzsyfNR7fnZmnN253W86LqM5MWU21Ti/JzKEuzPlO4B+oDk4PA1sy\n8+et5PTaMpJUIL+hKkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtdkgpkuUtSgSx3SSrQ/wOvCzbZE05m\nJwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# weightの分布をみる\n", "davis.weight.hist(bins=20)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAFhCAYAAADXxalcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmQnHd95/F3d899yLKk0WFJvmL7azsgbIwRWIBtca0J\nFASbEOJ4vUASNoFdCmqzZUicOE5Sm2UX18YQKikTE8dFCGyiJRgCGNnGkY1tmdhIBktfyZeuGc2M\nNJI1V8/08ewfz9Otnp57NN1Pz8znVaVS99PPdH+nNfr0b37XkwiCABERqa5k3AWIiCxFCl8RkRgo\nfEVEYqDwFRGJgcJXRCQGCl8RkRjUVfoFzOwLwFuAFPAXwNPA14B6YBT4TXfvMbObgU8DOeAed7+3\n0rWJiMQlUcl5vmZ2HfDf3P29ZrYCeBZ4GPieu/+Tmf0ecC5wJ/AM8AYgSxjQb3X3kxUrTkQkRpVu\n+T4KPBXdPgm0AL8LjETHeoErgc3ATncfADCzx4AtwPcqXJ+ISCwqGr7uHgDD0d3fAv7V3dMAZpYE\nPgn8CbCWMIgLeoF1laxNRCROVRlwM7P3Ax8FPhXdTwL3A9vd/ZEJviRRjbpEROJSjQG3dwOfA97t\n7v3R4a8B7u5/Ft3vZGxLdz3wxFTPGwRBkEgoo0UkVnMOoUoPuC0DdgBvd/dj0bGbgWvd/XdKzmsC\ndhMOuOWBnwJXl4T1RILe3qkerq6OjnZqqR6ovZpUz/RqrSbVM7WOjvY5h2+lW74fBlYC3zIzCD8l\nNgInzewRIACed/dPmdltwIOE4XvHNMErIrKgVXrA7R7gnhmeuw3YVsl6RERqhVa4iYjEQOErIhID\nha+ISAwUviIiMVD4iojEQOErIhIDha+ISAwUviIiMVD4iojEQOErIhIDha+ISAwUviIiMVD4iojE\nQOErIhIDha+ISAwUviIiMVD4iojEQOErIhIDha+ISAwUviIiMVD4iojEQOErIhIDha+ISAwUviIi\nMVD4iojEoK7SL2BmXwDeAqSAvwCeBu4nDP4u4BZ3z5jZzcCngRxwj7vfW+naRETiUtGWr5ldB1zu\n7tcANwD/B7gT+LK7Xwu8CHzMzFqA24GtwPXAZ8xseSVrExGJU6W7HR4FPhTdPgm0AtcC34mOPQC8\nE9gM7HT3AXdPA48BWypcm4hIbCra7eDuATAc3f048D3g3e6eiY71AOuANUBvyZf2RsdFRBalivf5\nApjZ+4GPAe8CXih5KDHJl0x2fIyOjvYzrGx+1Vo9UHs1qZ7p1VpNqqcyqjHg9m7gc4Qt3n4z6zez\nRncfAdYDR4BOxrZ01wNPTPfcvb39lSh5Tjo62muqHqi9mlTP9GqtJtUztTP5IKj0gNsy4AvAe939\n1ejwduDG6PaNwA+AncAbzGyZmbUB1wA7KlmbiEicKt3y/TCwEviWmSWAALgV+Fsz+wRwALjP3XNm\ndhvwIJAH7nD32vl4k2nlg4DHd3dxuHeQDR2tbNm0jmRiRr1HIktSpQfc7gHumeChd01w7jZgWyXr\nkcp5fHcXDz97BIB9h08C8NbXnRNnSSI1TSvcZF4c7h2c8r6IjKXwlXmxoaN1yvsiMlZVpprJ4rdl\nUzhZpbTPV0Qmp/CVeZFMJNTHKzIL6nYQEYmBwldEJAYKXxGRGCh8RURioAE3qQnzvUJOK+6k1il8\npSbM9wo5rbiTWqduB6kJ871CTivupNYpfKUmzPcKOa24k1qnbgepCROtkDuTftvC8x3qHWA4neVQ\nzwA7dnWq71dqhsJXasJEK+R27Ooc028bBAGJRGJGYVx4vtLn2H8k3FJafb9SCxS+UrPK+2l37u1h\nMJ0FZj6Ipr5fqVXq812i8kHAjl2dfGP7fnbs6iQfBHGXNM50/bSHegem/R7U9yu1Si3fJWohTMUq\n7wcOgEeimgGG09lpvwfttia1SuG7RE3163itLFAo7wfOBwEJTgfpoZ6BMedP1KWg3dakVil8l6gN\nHa3F1mLhfkEttIrzQcBjuzrZubcHgDdetoa3bFo3po4duzqLg2gw9y6FwofNoZ4BhkeyNDfVsbGj\nTTMjpKIUvkvUVL+O18Ig1eO7u3jgJwfoHxoFoLtvmARjPwTmq0uh8GEzMJShf2iU9pYG9h/WzAip\nLIXvEjXVr+NTtYrnS7G1Gc3DbW6sY+Pq063Nw72DjGZz5PIBATCYznCod2w3w1Tfw2y6TgofLqPZ\nXMnf9ZoZIRWl8JVxqjFINWFrs2Qe7oaOVvL5gHw+nMGQyeYZjqaZzeb5Yfquk8KHTUNdipHRHA11\nqeJxkUpR+MqsTNaiLD2+vqMVgoAjx4YmbXWWtjZz+YBTUfdCYRBty6Z1PPn8UV7u6gegtame5saZ\n/7jOpOsknw+n2x3qGWDDqlaaGlOkR3Jj+nxFKkXhu4iUB+MHtl4yp+d5bHcXDzz+CqPZsBUYAG+L\nWo2TtShLW7I7dneSSiZYtbx5wlZnPggYSmfoO5VmNJMrtm77h0YZHglbt8lEgjddvpahkVzx6zau\nbpvx9zCTrpOHnj5Y/F4Atl65Xn28UjUVD18zew3wbeAud/+Kmb0N+HMgAwwAt7j7q2Z2M/BpIAfc\n4+73Vrq2xaY8GNvbm7jiwhWzfp6de7qLA10jozl27ukuhu9kLcrDvYPFLoRsPoAgYGAoQ1vL+L7T\nQl9vQ12K9GiO+rokDfVJGuvraG46/SN5Jt0fM/naV46emvB7EamGioavmbUAdwPbSw5/EfiIu79g\nZp8DPmFmXwZuB94AZIGnzWybu58c96QyqfLweOXoqTmF71Qma1Fu6GjlyeePApAASCSKA1flrc5D\nPQMMDmejlnWSRAJWntUMwMaO063bieb57tjVOau9HaZy/tpl7NrXO+57EamGSrd808ANwG0lx3qB\nDuAF4GxgL7AZ2OnuAwBm9hiwBfhehetbVMqD8fy1y+b0PG+8dDXdfcPFboc3Xrq6+NhkLcotm9ax\n79BJdr90nPq6JASwdmULmy9bM67VOTySLbasAc5d08YF686atnU73/OP3371ufT3p7X6TWJR0fB1\n9zwwYmalhz8LPGpmfcAJwmD+MGEoF/QC+p8wS+XB+Parz+X48YFpvmq8azatY//hVznUM8DG1W28\n6bVriy3OwmBauWQiwUd/5bIZTe9qbqqjrbmeoah/t74+xYffftG0Cxrme/5xMhn/6rf56qeXhSeO\nAbcvAe939yfN7AvAJ4FjZedoWdEclP+qnUzO7W184rmjHD42SCKZ4PCxQe7/vnP4WBh0z0S/pre1\n1I9rfc50Ke/Gjjae3XesONDW0zfM47u7pv3aasw/rrb56qeXhSeO8N3k7k9Gt7cDvwHcC7yv5Jz1\nwBPTPVFHR/v8V3cGaq0emFlN2WyeL//fn/Fy16tcsO4smptSYddBpOvEUPF+Np8HKN4/Pjg66Wvk\n8wEPPX2QV46e4vy1y3j7yjY6Otr5wNZLePaF47zc9SqN9SnaWxs4NjDKz17qO33u1eeO+/D4wNZL\naG9vmvKc2Yr73+z44OiY9/qVo6d45+bzYqxovLjfo3K1Vs9cxRG+XWZ2qbvvBa4G9gM7ga+a2TIg\nD1xDOPNhSr29/RUtdDY6Otprqh4Ia+ruOTXhvgVvfu3asIXbO8jLXa9ysDvsnjjQ1U9bcz11qSSt\nzXUkEgnWr2xh14vHGc3mSQBtzfVksmEIr2htYNt2H9dv+vjuLp7a00133zCtTSke2nmQr/9wD2e3\nNXLe2mW0NqZY3tYIhOHfd2KI3fvDVvWufb3096fZsmnduG6MKy5cUWwZzqVLpfz9ifvfbGVrQ/G9\nhLCfPu6aStXCe1SqFuuZq0rPdng94eyG84CMmd0E/GfCoB0F+oCPuXvazG4DHiQM3zvcvXbe4QVs\nsn0L9h06WexK6Dw2WOwCyOcDhkayLGtpoK25ns2XrcEPniA9miMg7A86u72RC88JB8iCIODhn3UC\njOkSePjZI/SdSjMymqN/KFyhNjSSpe/VEbqOD7HyrCY2drTR0lQ/6Q5ltbDBT6XNVz+9LDyVHnB7\nBrh+gofeMsG524BtlaxnKZpo34KBIdj94nEa6lO0tdTTUJdkeDR8PACCaEDtnFWtvPV157D93w+P\nec7RbJ6PvONiAL6xff+ErwdQn0oylM+Si4KdIAzvTDZPIpHgxMAILU31QNh/W75DWS1s8FNp89VP\nLwuPrmSxyBUGpQr7FeTzAf1Do+SD8O+BoQwrz2riwnXtNDekSCTCgOwfGi3updBYnySfDwiivRYa\n65Pjnr/0fuFYIpq9kIoCpTCZob4uycBQhu6+YfYdPhm2bhMJtl65nks2LGfrlevZsmmdrkIhi5qW\nFy8Qc93gvHgV36jP92BPP32nRgiCsJ9xJJPj0nOXc+t7LuVbD73As/uPFef3FvZSOG9tO13HhxjN\n5mmoS3Le2vZxzz/RXNmHohZza3Md3X1DBAEsb2vgly9cwd4DJ8O5vkPh40d6B4ut6Zk8t8hCp/Bd\nIGbb/5kPAn701AF+8UIvB7r7GcnkOXd1GxtWtdLdN0wuH5CLWrI/e+EYJ761i7PbGmltrqMtcbor\nYMeuTrqOD9HUUMeq5eHxc1e3j/swKJ+nW6it0N8cBHBWWwO5POw9cJJTg6OMZnKMZvLF1yqnq1DI\nYqbwXSBm2//5+O4udjzXRdexQYbSWZLJBN194UBXe0sDJwdGAMjlA4bSWV7qPMWaFS1jBsECwvAs\n9AG3NtUVV6zN5MOg0FJ9qKTPuHRlW0N9imQywZoVzWrVypKjPt8FYn1HKwND4U5gx04Oc/jYwKRX\n7M0HAU/t6ab35HBxlkLhrJFMjraWeurrkuNWsgRBwImBEQIC9h06yUM/PcTAUAYIF1WsX9XGW193\nTnGz81JTXT/t7VdtoK2lntHs2DnCyWSCFcua2HzZGl2uR5YctXwXiuD0puKZbJ7u40PFlmd5i/Px\n3V109w2THskVW62FaLvs3LO5ZONynnz+KAe7B8jlAzLZPK1NdQwOZxkcztLdN0z/UDj5vzAHta1l\n7AY5s1ltVmjVPvviMQ4dHaCtuY6B4SxrVjRPuPeDyFKg8F0gjhwbilqP4f63mVwYik/t6R43CHe4\nd5C2lnpSqQTDI1mSiQTLWhs4d3Ubt77nUuqSyWLXwYGjp/j5y30MjeRIJROc1dbAif6R6FXDfuET\n/WlSSbj6l6ffYAcmv0TQH3/8zXznx/tnPWhYK1dTFplPCt8FYqJL3QwMZRgYyjCYzrLv8EmCICCR\nSHC4t59jJ4epq0vSUJfifdecx9uuWA+M35YxPZrj5EDYD5vLB2RzeZoaUgyls2TS+WJ3Re/JNHd9\n42d8/pY3AKe7FArB+M2HXigG42SXCGpvb5rTANpSWGwhS4/Cd4EonzLW3FRH57FBBoYzxXN27u1h\nMB1u15gezdEQQF1j8vQEW8YHWf/g6QGwZAIa61OsWdFCejTHYNk107r7hsfVNVEwTnZByrnuL7wU\nFlvI0qMBtxpSaJV+Y/v+SQfTEokEl2xczq+//WI2X7amuJChVCabJ5VM0FCfpK2lniMlYXWoZ2DM\nwF16tGQFGrBiWSMDQxmaGupoqBv747FmRfO415ooGMsXdhT+nuv+wlpsIYuRWr41ZKpfryd6rLzf\nNQAeefZIsWuisX78VXgLG5kX5vg2N6VoakjR2JBiRXsjmWy+2HouDIwlEgk2rm7lv/3GleNqnmjg\nrdhKL+vzneu+BVpsIYuRwreGTPXrdentIJpKVr7AIZvPs//QSQ5297O8rY3RbJ7+wVH84Ak2v2YN\nT/28m4M9A9TXJQkyOUgmirMgMpkwdOvrwt3MIJzX+/63XDjlAFd50B7qGeDx3V0Tfs1c9y3QYgtZ\njBS+NWT9qhae2dfLaDYXbkqTzvAP2/cxnM5y/NQwx04Ok0wmyOeDMQNt+XyeF46cYs/BE4yM5lix\nrJGjfUOMZHIkEwl+6r30nBzmRP8og+kMmWw4qBZOWwvIZHMkkwmGRrKkkgkaG1I01KXYetWGSUOv\nOKMh6oPuO5XmYM8AyWSCZ/cfIwiC4iCfiIyn8K0lJS3FkUyOvQfDX+cLc24LV/oNTz197oM/Pczx\nV9Nko81vINx5LCjsAQkc6h0sPgbQ2JDiiotWsevFY+TyAalkorjkuGiCPueC8hkN+SAglwtIJhOM\nJHPs3Nuj8BWZgsK3hhyJ5udCPX2n0sXZAhCGaSqZoL4u3GHs1NBocWrZQDpDLh+QIFzJVtgAZyQT\nLrLI5oLiXrx1qSiNA7hk43IAnt7bUzhEc0OKFcuawnqODU1aa/mMhuD0rpEiMgOa7VBDSgfGGupS\nxT/h/fCfKh+tSGuoSzI0kiU9mqUu6oqAsF911VlN3HTdL3HtletJJBLFQAyAbO50PD787BEuWr+M\nqy9dzdoVLVy4rp0VyxonrGeyWgv1NTWE+zQ01CVpb2ngjZetOeP3Q2QxU8u3hpSO6q9f1QKJBIej\ngaymhhTp0RwHewYYSmdpba6j71Sa9GgOgvDX/ebGFK/7pVXccoPx1M+7aW1pYKJxsuXR7mUAnceH\n+fh7LwfG9+Me6g33j5ho8Kx83nFDQ5JfvNTH0EiWs9sbeNNrFL4iU1H41pB8EG5oc6hngKF0prgU\nuNSOXZ3FvtahdJZ8dHWIgKC4afnd/3c33X3DLG9vIJvLj/n6Quie6B8Jp5o1JPnb7z5fvLZb+Y5l\n+w+/WuyyKF/eWzoY97fffb64Uu5g9wD3f9+LoS4i4yl8a8h9/7q32P96tC/sby0PsC2b1rHv0El2\nv3Q8bNUGp/tZB4YzPL23pzgjIpVKkEolyUVXHE4AbU31ZHJBcYOeQz2DvNTZX7y2G4yf8rZzT3dx\ntdtky3vLr8FWfl9ExlL41pCJAiybz3Pfv+7lUM8AG1e38ZF3X8xzLx+nf+j0suJCh0D5arf+odFw\nBkIivJRPPoATA6PUpRLkohZx4YrEhSXAT+3pBmBgKBMN/o030fLejavbih8YhfsiMjmFbw2ZKMDK\nW8O7XzxG//DYPReSqQRN9Ska6pMMDmdpbaojPZojH4TXWxvJ5EtmIQQMpXOFvopiV0XpRj2liyw2\nX7aGIAh4JLpCMUw8EHfrey4FKH5IFO6LyMQUvjWkPMAu2nAW//zoi2G/bSIBQVDckLzUOStbecdV\nGzjY3V+8ZFBDfYrmxnAhRd+pEUYzeRobksWrBwPU1SXpOKuR89edNWajnkQiMWbz9Hxxt7TJl/fW\nJZPq4xWZBYVvFZTvR/vm167lieeOjtmhrDDY9dFfuYzHd3fx5PPdfGP7/tNhO8WCh0w013bjmnZe\n6DxFMplgNJ0PF01E/bsN9cloZgTkS1q85687qxiahcG8gkILV8t7ReafwrcKyjfF2XfoJIePDY7d\n7zYa7IJw/m1n72Bxw/Tp9A9lePjZI7Q2nf7nbGup5+z2JvpeLWwDGUSXFApIJMLFFm3NDTSXfI02\nsBGpHoVvFZQPUB3s7mdoJMepoVHy+aA42FV63kTbSUKx92GMoXSWgaEMLY2n/zkHhjKcvayRjavb\nOHxskL5T6WiFXCpqCadoa6lnY8fpgbG5tHB1lQmRual4+JrZa4BvA3e5+1fMrA64D7gIOAXc5O6v\nmtnNwKeBHHCPu99b6dqqpXzbxcaGFN0nhslH2zoWVqcVfs3fd/gkTQ2pcZuZQxiQubL0DQhnNtjG\ns3jT5Wt4ak932KoeHKXv1TQbO9pobaqju2+YlqYUQ+ncvF0/TVeZEJmbioavmbUAdwPbSw7/NtDj\n7jeb2W8BbzWzh4HbgTcAWeBpM9vm7ifHPekCVP7r/MGeAU70jzKSyRIEsGJZExtWtXKoZ4D1Ha1s\nWNVKPh9e1LKvf4RcyZLgIAiiqWJBcQZDXSrc8+FQ7yB27tmcs6qVwXS4D28ikaClqb7YlzzfLVRd\nZUJkbird8k0DNwC3lRx7H/BHAO7+VQAzux7Y6e4D0f3HgC3A9ypcX1WU/zq/Y1cnLxx5lTbCebQb\nVrVy+NggQRDwxC+OkssHtDbV09pcxwrC66eVtnVL92cInz+cwTCUzvLws0fYsKp6V36YzVWMReS0\nioavu+eBETMrPXw+8B4z+19AF/BJYC3QW3JOL7BoR3vKW8KFxRWDw1mGR8OZC/mhcKluc0OSlqY6\nhtJZAiae9JAgQVtLQ3F+bnNjHVuvXM/xwVFWtjaMWzI8n90DGqQTmZs4BtwSwB53v9PM/gD4HPDs\nBOdMq6Ojfb5rOyOzqeeD7zh9PbMfPXWAV7r7yebzlF7sYWgkQ0tTM63N9eTyASOZsds3JhJhq3r1\nihYaG1IEQUD/UIZjp0a4/JdW8YGtlxSvHnF88EBxL+Dw/ui8vX+l38t0FvK/WbXUWk2qpzLiCN+j\nwL9Ft38I3AF8l7A7omA98MR0T9Tb2z/ftc1ZR0f7nOvZdMHZ9PeneWpPd3SF4IDBdJZkIkEyAblc\nwPK2BvqHMqdDOAiXEzc1pLj+9etJJRI8taebk7lRTvSn+Zd/exGgeLXgla0NxcUVhfvVfv/O5D2q\nhFqrB2qvJtUztTP5IIgjfL9P2A/8d8BVgAM7ga+a2TIgD1xDOPNhSSj0Cb/5tWuL+zg01KcYGc1y\n/FSaIICWpiauuOis4taSJwZGAHjjpat5y+vOIRmtQCudIVF6qXZ1D4jUlkrPdng98EXgPCBjZjcB\nvwHcbWYfB/qBW909bWa3AQ8Shu8d7l47H2+TKJ3jetmFK9l0wdnjZhDkg4DHdnexM9qwpjQsyz3x\n3FEOHxskkUxw6tRouCINoqloaZLJBFuvXD9pX+26VS385OddZLJ56uuSbFx9+lNZq9REakulB9ye\nAa6f4KFfm+DcbcC2StYz30oHsV4+eor+/vS4gHt8dxcPPP4K/dEAWnffMIlJgrB0mlYyumRQruQK\nFeXnlHvh0MlwFVsQXottz8vHueqilWf2TYpIRegyQmdgJnNcD/cOll2LLTdpgG7oaCUIwisTZ6Jr\ntrU315NKJoqX65lqKtfh3kFSyQR1qSSpZIJXjp6ay7clIlWg8D0D5UE4UTBu6GgtBieEWzdOFqBb\nNq1jY0cbo9kcrU31NDWkWLOihasvXc2Vl6xi65Xrp+yrLd9D94J1Z83m2xGRKtLeDmegdBCr0Oc7\n0TkBjOnznSxAk9FqtMLVgwE2dLTxkXdcXLyfDwJ27OqccKVa+ZaUn/rQFZw4oRVnIrVI4XsGSgex\nyqfAlG8489kPXzGj5bzTrRibarFE+Z66dXX6xUakVil8K2SuK8qmmxKmvRREFgeFb4VMFZJTbcM4\n3ZQw7aUgsjgofCtkqpA8k30WtFhCZHFQ+FbIVCE5k66DyVrHWiwhsjgofCtkqpCcSdeBNikXWdwU\nvjGYSdeBBtZEFjfNRapRM1nAISILl1q+MZhJl4IG1kQWN4VvDGbSpaCBNZHFTd0OMVCXgoio5RuD\nyboUplp8ISKLi8I3BpN1KWh6mcjSoW6HGqLpZSJLh8K3hqgvWGTpULdDDdH0MpGlQ+FbQzS9TGTp\nULeDiEgMFL4iIjFQ+IqIxEDhKyISg4oPuJnZa4BvA3e5+1dKjr8b+L67J6P7NwOfBnLAPe5+b6Vr\nExGJS0VbvmbWAtwNbC873gjcBnSWnHc7sBW4HviMmS2vZG0iInGqdLdDGrgB6Co7/nngy8BodH8z\nsNPdB9w9DTwGbKlwbSIisalo+Lp73t1HSo+Z2SXAJnf/55LDa4Hekvu9gFYYiMiiFccii7uA/xLd\nnmzLLm3lJSKLWlXD18zOAQz4upklgHVm9gjwx8D7Sk5dDzwx3fN1dLRXpM65qrV6oPZqUj3Tq7Wa\nVE9lVDN8E+7eCVxcOGBmL7v79WbWBHzVzJYBeeAawpkPU+rt7a9YsbPV0dFeU/VA7dWkeqZXazWp\nnqmdyQdBRcPXzF4PfBE4D8iY2Y3AB929cN30AMDd02Z2G/AgYfje4e618w6LiMyzGYWvmS0vCczC\nsQvc/eWpvs7dnyGcOjbZ4xeW3N4GbJtJPSIiC9204WtmSeD/mdlWTg+E1QPfAV5bwdpERBatKaea\nmdlHgL3AtYQrz7LRn0HgYMWrExFZpKZs+br7N4BvmNkd7n5HdUoSEVn8Zjrg9gUz+1VgOSVzcLX/\ngojI3Mw0fL9H2N1wuORYACh8RUTmYKbh2+Tub65oJSIiS8hM93Z4xsxWVbQSEZElZMqWr5ntIOxe\nqAP2mdlewu4HANz9bZUtT0RkcZqu2+EPq1KFiMgSM2W3g7s/6u6PAqkJ/gTRRjkiIjJLMx1w+wPC\nzc33ES62MODfgQvM7H+4+19VqD4RkUVppgNuB4Gr3H2Tu18JvAH4OXAR8B8rVZyIyGI10/C9yN1/\nUbjj7s8Dl0eX/MlVpDIRkUVspt0OQ2b2v4Efc3q/3YboCsQDFapNRGTRmmnL9yPAMPAJ4JNAM3AT\n8DJwS2VKExFZvKab55tw9wA4SXipnzHcPV+pwkREFrPpWr4PRX9ngUzJn8J9ERGZg+m2lNwa/V3R\nS8yLiCw1M72M0NnA54G17n6Lmb0PeNLdeytanYjIIjXTFu1XgUNA4ZprjcB9FalIRGQJmGn4drj7\n3cAogLv/E9BSsapERBa5Gfflmlk90aXezWwN0FqpokREFruZLrL4MvA0sNbMvgO8Efh0xaoSEVnk\nZhq+3wWagLVAmjCI1fIVEZmj2VzDLQMcKTl2AbqGm4jInMzmGm7Xz+UFzOw1wLeBu9z9K2a2kTC0\n6wkH8H7T3XvM7GbCrowccI+ujCwii1lFr+FmZi3A3cD2ksN/Cvy1u19HGMqfjc67HdgKXA98xsyW\nz/b1REQWikpfwy0N3ADcVnLsd6PjAL3AlcBmYKe7D0Sv+xjh5u3fm/F3IiKygFT0Gm7RxjsjZlZ6\nbBjAzJKEO6T9CeFAXulquV5g3Zm8tohILZtub4dHK/GiUfDeD2x390fM7CNlpyQq8boiIrVipgNu\n8+1rgLv7n0X3Oxnb0l0PPDHdk3R0tFegtLmrtXqg9mpSPdOrtZpUT2VUPXyjWQ0j7n5nyeGngHvM\nbBmnr5RtfdbMAAANk0lEQVQx7SKO3t7+yhQ5Bx0d7TVVD9ReTapnerVWk+qZ2pl8EFQ0fM3s9cAX\ngfOAjJndBKwG0mb2COFg3vPu/ikzuw14kDB873D32nmHRUTmWUXD192fIZw6NpNztwHbKlmPiEit\n0CbpIiIxUPiKiMRA4SsiEgOFr4hIDBS+IiIxUPiKiMRA4SsiEgOFr4hIDBS+IiIxUPiKiMRA4Ssi\nEgOFr4hIDBS+IiIxUPiKiMRA4SsiEgOFr4hIDBS+IiIxUPiKiMRA4SsiEgOFr4hIDBS+IiIxUPiK\niMRA4SsiEgOFr4hIDBS+IiIxqKv0C5jZa4BvA3e5+1fMbANwP2HwdwG3uHvGzG4GPg3kgHvc/d5K\n1yYiEpeKtnzNrAW4G9hecvhO4Evufi3wIvCx6Lzbga3A9cBnzGx5JWsTEYlTpbsd0sANhC3cguuA\nB6LbDwDvBDYDO919wN3TwGPAlgrXJiISm4qGr7vn3X2k7HCru2ei2z3AOmAN0FtyTm90XERkUap4\nn+80ErM8PkZHR/s8lnLmaq0eqL2aVM/0aq0m1VMZcYRvv5k1Ri3i9cARoJOxLd31wBPTPVFvb39l\nKpyDjo72mqoHaq8m1TO9WqtJ9UztTD4I4phqth24Mbp9I/ADYCfwBjNbZmZtwDXAjhhqExGpioq2\nfM3s9cAXgfOAjJndBNwM3GdmnwAOAPe5e87MbgMeBPLAHe5eOx9vIiLzrKLh6+7PEE4dK/euCc7d\nBmyrZD0iIrVCK9xERGKg8BURiYHCV0QkBgpfEZEYKHxFRGKg8BURiYHCV0QkBgpfEZEYKHxFRGKg\n8BURiYHCV0QkBgpfEZEYKHxFRGKg8BURiYHCV0QkBgpfEZEYKHxFRGKg8BURiYHCV0QkBgpfEZEY\nKHxFRGKg8BURiYHCV0QkBgpfEZEY1FX7Bc2sFfh74GygAbgTeB64n/DDoAu4xd0z1a5NRKRa4mj5\n/idgr7tvBT4E/CVhAH/Z3a8FXgQ+FkNdIiJVE0f4HgNWRrdXAL3AtcB3omMPAO+IoS4Rkaqpevi6\n+zeB88xsP/Bj4PeB1pJuhh5gXbXrEhGppqqHr5ndDBxw94uBrcBflZ2SqHZNIiLVVvUBN2AL8EMA\nd3/OzNYBg2bW6O4jwHqgcyZP1NHRXrkq56DW6oHaq0n1TK/WalI9lRFH+L4AvAn4f2Z2HtBP2P1w\nE/B14EbgBzN5ot7e/gqVOHsdHe01VQ/UXk2qZ3q1VpPqmdqZfBDEEb5/A9xrZj8GUsAnAAf+3sx+\nBzgA3BdDXSIiVVP18HX3QeDDEzz0rmrXIiISF61wExGJgcJXRCQGCl8RkRgofEVEYqDwFRGJgcJX\nRCQGCl8RkRgofEVEYqDwFRGJgcJXRCQGCl8RkRgofEVEYqDwFRGJgcJXRCQGCl8RkRgofEVEYqDw\nFRGJgcJXRCQGCl8RkRgofEVEYqDwFRGJgcJXRCQGCl8RkRgofEVEYqDwFRGJQV0cL2pmNwO/D2SA\nPwKeA+4n/DDoAm5x90wctYmIVEPVW75mtoIwcK8B3gt8ALgT+JK7Xwu8CHys2nWJiFRTHN0O7wB+\n5O5D7t7t7p8ArgMeiB5/IDpHRGTRiqPb4Xyg1cz+BVgO/AnQUtLN0AOsi6EuEZGqiSN8E8AK4FcJ\ng/iR6Fjp4zPS0dE+r4WdqVqrB2qvJtUzvVqrSfVURhzh2w38xN3zwEtm1g9kzKzR3UeA9UDnTJ6o\nt7e/gmXOTkdHe03VA7VXk+qZXq3VpHqmdiYfBHH0+T4IbDWzhJmtBNqA7cBN0eM3Aj+IoS4Rkaqp\nesvX3TvN7J+AJ4EA+CTwU+B+M/sd4ABwX7XrEpHFIR8EPL67i8O9g2zoaGXLpnUkEzPuzayaWOb5\nuvs9wD1lh98VRy0isrg8vruLh589AsC+wycBeOvrzomzpAlphZuILCqHewenvF8rFL4isqhs6Gid\n8n6tiKXbQUSkUrZsCpcJlPb51iKFr4gsKslEoib7eMup20FEJAYKXxGRGCh8RURioPAVEYmBwldE\nJAYKXxGRGCh8RURioPAVEYmBwldEJAYKXxGRGCh8RURioPAVEYmBwldEJAYKXxGRGCh8RURioPAV\nEYmBwldEJAYKXxGRGCh8RURioPAVEYlBbBfQNLMm4OfAncDDwP2EHwZdwC3unomrNhGRSouz5Xs7\ncDy6fSfwJXe/FngR+FhsVYmIVEEs4WtmBlwKfA9IANcCD0QPPwC8I466RESqJa6W7xeBzxIGL0Br\nSTdDD7AulqpERKqk6uFrZrcAP3H3A5OckpjkuIjIopEIgqCqL2hm/whcAOSB9cBo9NAvu/uImb0N\n+JS7/1pVCxMRqaKqz3Zw918v3DazPwJeAa4BbgK+DtwI/KDadYmIVFPc83wLXQx/DNxqZo8CZwP3\nxVeSiEjlVb3bQURE4m/5iogsSQpfEZEYKHxFRGIQ294Os2VmXwDeAqSAvwCeJub9IGptfwozuxn4\nfSAD/BHwXFw1mVkr8PeEA6gNhO/R83HUY2avAb4N3OXuXzGzDRPVEb1/nwZywD3ufm+V6tkI3AvU\nE069/E1374mrnpLj7wa+7+7J6H5c708d4SD8RcAp4CZ3f7Va9UxS09uAPyf8vzZA+DM0q5oWRMvX\nzK4DLnf3a4AbgP9D+J/5yzHvB1Ez+1OY2QrCwL0GeC/wgZhr+k/AXnffCnwI+Eti+DczsxbgbmB7\nyeFx70t03u3AVuB64DNmtrxK9fwp8Nfufh3hf/DPxlwPZtYI3AZ0lpwXVz2/DfS4+2bgm8Bbq1XP\nFDV9Efho9PP9BPCJ2da0IMIXeJTwPzDASaCVcD+I70THqr4fRA3uT/EO4EfuPuTu3e7+CeC6GGs6\nBqyMbq8Aeonn3yxN+IHdVXLsOsa+L+8ENgM73X3A3dPAY8CWKtXzu8C26HYv4fsWZz0Anwe+zOlF\nUHHW8z7CNQC4+1fd/btVrGeymnqBjuj22YQ/77OqaUGEr7sH7j4c3f04YeDFvR9Ere1PcT7Qamb/\nYmaPmtlWoCWumtz9m8B5ZrYf+DFhd0jV3yN3z7v7SNnhiepYQ/gfqqC3EvVNVI+7D7t7YGZJ4JPA\nPwBr46rHzC4BNrn7P5ccjq0ewp/t95jZI2b2D2Z2drXqmaKmzwLfNrM9hN2hfzfbmhZE+BaY2fsJ\nf1X9FGP3gKjqfhA1uj9FgrCF+avAR4GvEe97dDNwwN0vJvw17K/KTqmVPTwmq6Pa71eSsB96u7s/\nEnM9dxGGy1SvW816EsAed78e+AXwuZjrAfgS8H53v4ywhfvJCc6ZsqYFE75R5//ngP/g7v1Af9Qv\nBeEeEZ1VLOdXgPeb2ROELfHbgYEY6wHoJvxAyLv7S0Dc79EW4IcA7v4cYQtgMOb3qKD8fTkS1VLa\nSql2fV8D3N3/LLofSz1mdg5gwNejn+91ZvYI4XsU1/tzFPi36PYPgctjrgfC3wyejG5vB66abU0L\nInzNbBnwBeC97v5qdHg74T4QUOX9INz91919s7u/Gfgq4QDOdsL9KapeT+RBYKuZJcxsJdAWc00v\nAG8CMLPzCD8MfhRjPaUm+tnZCbzBzJaZWRvhwOWOahQT/ZYw4u53lhx+KoZ6Eu7e6e4Xu/s10c93\nV9TijO39Ab5P2OcKYch5zPUAdJnZpdHtq4H9s61pQSwvNrPfJtz/YR9hUz4AbgX+FmgEDhCOPOZi\nqO2PgZcJP5Hvj7Oe6H36LcL350+Bn8ZVUzTV7F7CvtQU8IeE/2n+vpr1mNnrCfvnzyOcFnQEuJlw\n6tKYOszsg8B/J9xx7253/8cq1bOacFCnn/Df7nl3/1SM9XzQ3U9Gj7/k7hdGt+Oq5zcIZxusI3yP\nbnX33mrUM0VNnwf+N+GAZB/wMXc/NZuaFkT4iogsNgui20FEZLFR+IqIxEDhKyISA4WviEgMFL4i\nIjFQ+IqIxEDhK0uOma0xs29Oc86tZnb/JI/dXJnKZClZMPv5iswXd+8GPjyDU8dNgjezFOHWnV+f\n77pkadEiC1nQzOwl4IpoddE3gQF3/7iZrSFcRvw3wK8RrrLbC/we4e5Tj7n7RjO7gHAVYJ5wg/73\nEO7d8Vbgg4Sbd18OvOLuN5rZ3wG/DvzY3f9DFb9VWWTU7SAL3XbCLf0gXMp8YXT7esLtIj/g7m9z\n9y3Aq4TLr+F0q/ZO4B/d/W2E+2NcXPLclwO/5e5XAa81sysJl7n3KHjlTKnbQRa67cC1ZnaIsGV7\nVnSZoOsJN8z5pJk9TLgnSAunNwcvuAL4nwDu/kMzGyh57OmSfVyPAMsJ1/GLnDGFryx024H/Chwm\n3LR9BeEVM95EeH2977j7fy39gmiXtYIkYZdDQWk/XLbstWplD2JZBNTtIAuau/cR/hy/hzB8dxAO\npnUCjwM3RDusYWa/a2aby55iD+HWf5jZOwm34pxKnvCCoCJnROEri8GPgfPd/Wi0cfubgB+6+zOE\nV9D4sZn9G2GLeFfZ194BfMrMHooeP8z4Fi+cbhF3AkfN7Gkza57370SWDM12kCXNzK4CGt39J9EM\nieeB1XHsDS1Li/p8ZakbAP4yvBg19cDvKHilGtTyFRGJgfp8RURioPAVEYmBwldEJAYKXxGRGCh8\nRURioPAVEYnB/wf3bNbbFaLZtwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 身長(height)と体重(weight)の関係をみる\n", "sns.lmplot('weight', 'height', data=davis, fit_reg=False )\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

Sageとpandasの機能を使って共分散行列を求める

\n", "\t

\n", "\t\tSageとpandas(numpy)の機能を使って共分散を求めみましょう。\n", "\t

\n", "\t

\n", "\t\t平均mxを求めてデータの平均からの差Xcを求めます。pandasではこれらをmeanとマイナス記号で処理することができます。\n", "\t

\n", "\t

\n", "\t\t数式では共分散行列$\\hat{\\Sigma}$は、以下の様に定義されていますが、\n", "$$\n", "\t\\hat{\\mu} = \\frac{1}{N}\\sum_{n=1}^N x^{(n)}, \n", "\t\\hat{\\Sigma} = \\frac{1}{N} \\sum_{n=1}^N (x^{(n)} - \\hat{\\mu})(x^{(n)} - \\hat{\\mu})^T\n", "$$\t\t\n", "\tXcが列ベクトルではなく、行ベクトルなので、以下の様に計算しています。\n", "$$\n", "\t\\hat{\\Sigma} = \\frac{1}{N} Xc^T Xc\n", "$$\t\t\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# pandaデータフレームで、weightとheightの平均からのずれXcを計算する\n", "X = davis[['weight', 'height']]\n", "mx = X.mean()\n", "Xc = X - mx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\tnumpyの固有の理由だと思いますが、マトリックスの積に.dot()関数を使います。この辺はちょっと違和感があります。\n", "\t

\n", "\t

\n", "\t\tPnadasの機能を使うと共分散行列は一発で計算できるのですが欠損値の影響なのか値が一致しません。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " weight height\n", "weight 226.720 34.2040\n", "height 34.204 143.4696\n", " weight height\n", "weight 227.859296 34.375879\n", "height 34.375879 144.190553\n" ] } ], "source": [ "# Pandaで標本共分散行列を求めるには\n", "print Xc.T.dot(Xc) /len(Xc)\n", "# もっと簡単なのはcov関数を使う(でも値が微妙に異なるのは、欠損値があるからか?)\n", "print Xc.cov()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

Sageのマトリックスを使った異常度の計算

\n", "\t

\n", "\t\tSageのmatrixを使った異常度の計算方法を以下に示します。\n", "\t

\n", "\t

\n", "\t\t理論をベースとしているSageと実験値を処理するためのPandas, \n", "\t\tnumpyをどうやって融合するのかがSage使いの腕の見せ所です。\n", "\t

\n", "\t

\n", "\t\tと言ってもやっていることは簡単です。マトリックスの要素単位の処理をするには、\n", "\t\tSageのmatrixにnumpy()メソッドで、numpyのマトリックスに変換するだけです。\t\t\n", "\t

\n", "\t

\n", "\t\tRやnumpyでは多くのデータを一括して処理することが多く、要素単位の処理に四則演算子が\n", "\t\t使われますが、Sageではマトリックスやベクトルの定義に沿って計算が行われます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Sageのmatrixを使った計算\n", "# 中心化したデータをSageのマトリックスにする\n", "Xc =matrix((X - mx).values)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[ 226.7199999999998 34.203999999999986]\n", "[34.203999999999986 143.4696000000001]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 標本共分散行列を求める\n", "Sx = Xc.T*Xc/ Xc.nrows(); Sx" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEBCAYAAAB4wNK4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt8FPW9//HXzG4um5CEDSQEFAjEAwhCIUQpVunxyikg\nqLUqWiq0WhGw9qj8fFhqW6otgvahoo8eDqVAqeCxR4ueaBQvICIWiQkqiAqIXOQiIYiBXHd3vr8/\nZpOwJIEkuxCQ9/PxyGOTmd3Z73x2dt7z/c7sxjLGGERE5Ixnt3UDRETk1KBAEBERQIEgIiJhCgQR\nEQEUCCIiEqZAEBERQIEgIiJhp00gGGMoKytDH5sQETkxvG3dgMaUlBxqMK28/BA9epzFF1/sIjk5\npQ1adfqzbYv09GQOHCjHcRSsraEaRk81jF5TNczIiG7feNr0ECzLiriVlrNtC8uysG3VsLVUw+ip\nhtE7UTU8bQIhFp58Mp7vfCeZSy9NYv36M2rVRUSO64zZK65Z4+HBBxPYs8dmwwYPP/+5r62bJCJy\nSjljAmH3buuYf4uInOnOmEAYNixEly5O3d833hhow9aIiJx6TsmrjE6Ejh0Ny5ZV8MorXjp0MIwa\nFWzrJomInFLOmEAA6NTJMH68egYiIo05Y4aMRETk2BQIIiICKBBERCRMgSAiIoACQUREwhQIIiIC\nKBBERCRMgSAiIkArAmHZsmVkZWVx0003RUxfuXIltm2TlJREUlISPp+PpKQknn/++br7zJ49mz59\n+tC+fXuGDRtGcXFx9GsgIiIx0aJPKj/yyCPMnz+fXr16NTo/OzubrVu3NjovPz+f6dOns2zZMvr3\n788TTzzBqFGj+Pzzz/H59M2jIiJtrUU9BJ/Px9q1a8nJyWnxE82dO5cJEyaQl5dHQkICU6dOxbIs\n8vPzW7wsERGJvRYFwpQpU0hJafpftJWVlXHttdeSkZFB165deeyxx+rmFRUVkZubW/e3ZVkMHDiQ\nwsLCVjRbRERiLWZfbpeamsqAAQO4++67+cc//sGKFSv40Y9+hN/vZ/z48ZSWluL3+yMek56ezv79\n+xssy7Yb/ms4j8euu/V6dS68NY6sobSOahg91TB6J6qGMQuEQYMGsXz58rq/r7jiCiZOnMiCBQsY\nP348AMY07x9qp6cnN/jfyR5PCIDUVB+pqcmxafQZKjVV52yipRpGTzWMXqxreEK//jo7O7vuKqOM\njAxKS0sj5peWltK/f/8GjztwoLxBD6G8vBKAsrJKQiHPCWrxt5vHY5Oa6gvX0Dn+A6QB1TB6qmH0\nmqqh3x/dwXLMAuG5555j//79TJw4sW7axo0b6dmzJwB5eXkUFRUxbtw4ABzHobi4mFtvvbXBshzH\n4DiRvYnalQ6FHIJBbUTRUA2jpxpGTzWMXqxrGLMBqPj4eO69917eeOMNgsEgr7/+OgsXLmTSpEkA\n3HHHHSxatIj33nuPyspKHnroIRITExk5cmSsmiAiIlFoUQ/B5/NhWRaBgPtfx5YuXYplWVRUVDB6\n9Ggef/xxpkyZws6dO8nKymL27NmMGTMGgOHDhzNjxgyuv/56SkpKOP/88ykoKCAhISH2ayUiIi1m\nmeae6T2JSkoONZhWUXGY7OwubNu2m6Skdm3QqtOf12vj9yfz9dfl6qq3kmoYPdUwek3VMCOj6Y8F\nNIeu+xIREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoE\nEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCB\nICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKm\nQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREB\nFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREaAVgbBs\n2TKysrK46aabGsxbvnw5Q4YMIS0tjf79+7NkyZKI+bNnz6ZPnz60b9+eYcOGUVxc3PqWi4hITLUo\nEB555BF++ctf0qtXrwbz9u7dy5gxY5g0aRIlJSU8/vjj3HbbbXU7/fz8fKZPn87TTz/NV199xahR\noxg1ahSVlZWxWRMREYlKiwLB5/Oxdu1acnJyGsxbvHgxvXv35pZbbiE+Pp7LLruM0aNHM2/ePADm\nzp3LhAkTyMvLIyEhgalTp2JZFvn5+bFZExERiUqLAmHKlCmkpKQ0Oq+oqIjc3NyIabm5uRQWFjY6\n37IsBg4cWDdfRETaljdWCyotLaVr164R09LT09m/f3/dfL/f3+T8I9m2hW1bEdM8Hrvu1uvVufDW\nOLKG0jqqYfRUw+idqBrGLBAAjDFRza+Vnp6MZR0dCCEAUlN9pKYmt66BArg1lOiohtFTDaMX6xrG\nLBAyMjIoLS2NmFZaWkpmZuYx5/fv37/Bsg4cKG/QQygvd08+l5VVEgp5YtXsM4rHY5Oa6gvX0Gnr\n5pyWVMPoqYbRa6qGfn90B8sxC4S8vDwWLlwYMa2wsJAhQ4bUzS8qKmLcuHEAOI5DcXExt956a4Nl\nOY7BcSJ7E7UrHQo5BIPaiKKhGkZPNYyeahi9WNcwZgNQN998M9u2bWP+/PlUV1dTUFDAK6+8wu23\n3w7AHXfcwaJFi3jvvfeorKzkoYceIjExkZEjR8aqCSIiEoUW9RB8Ph+WZREIBABYunQplmVRUVFB\nRkYGL730EnfeeSeTJ08mOzubxYsX069fPwCGDx/OjBkzuP766ykpKeH888+noKCAhISE2K+ViIi0\nmGWae6b3JCopOdRgWkXFYbKzu7Bt226Sktq1QatOf16vjd+fzNdfl6ur3kqqYfRUw+g1VcOMjMY/\nFtBcuu5LREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQp\nEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQA\nBYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhI\nmAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEERE\nBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKI\niITFNBBs28bn85GUlFR3e9dddwGwfPlyhgwZQlpaGv3792fJkiWxfGoREYmSN5YLsyyLTZs20bVr\n14jpe/fuZcyYMTz11FOMHTuWVatWMXr0aPr06UNubm4smyAiIq0U0x6CMQZjTIPpixcvpnfv3txy\nyy3Ex8dz2WWXMXr0aObNmxfLpxcRkSjE/BzCfffdR/fu3fH7/UycOJHy8nKKiooa9ARyc3MpLCyM\n9dOLiEgrxXTIaOjQoVx55ZUsWrSIrVu3csMNNzBp0iRKS0sbDCOlp6ezf//+Rpdj2xa2bUVM83js\nuluvV+fCW+PIGkrrqIbRUw2jd6JqGNNAWL16dd3vvXv35uGHH+aqq65i2LBhjQ4lNSU9PRnLOjoQ\nQgCkpvpITU2OTYPPUKmpvrZuwmlPNYyeahi9WNcwpoFwtOzsbEKhELZtU1paGjGvtLSUzMzMRh93\n4EB5gx5CeXklAGVllYRCnhPT4G85j8cmNdUXrqHT1s05LamG0VMNo9dUDf3+6A6WYxYIH3zwAU8/\n/TSPPvpo3bSNGzeSmJjIiBEjWLhwYcT9CwsLGTJkSKPLchyD40T2KGpXOhRyCAa1EUVDNYyeahg9\n1TB6sa5hzAagMjMzmTt3LrNmzaKmpoZNmzbxm9/8httvv50f//jHbN++nfnz51NdXU1BQQGvvPIK\nt99+e6yeXkREohSzQOjSpQsFBQW8+OKLdOzYkYsuuogRI0Ywc+ZMMjIyeOmll3jyySdp374999xz\nD4sXL6Zfv36xenoREYlSTM8hXHTRRREnlo+et27dulg+nYiIxJCu+xIREUCBICIiYQoEEREBFAgi\nIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoE\nEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCB\nICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKmQBAREUCBICIiYQoEEREBFAgiIhKm\nQBAREUCBICIiYQoEEREBFAgt4jht3QIRkRNHgdAM27dbDBuWROfO7bjxRh8VFW3dIhGR2FMgNMNv\nf5vAp596MMZi+XIvc+fGt3WTRERiToHQDN98Y0X8ffCg1cQ9RUROXwqEZrjttgBerwHA7zfcdFOg\njVskIhJ73rZuwOlgxIggb71VwebNNoMHh8jKMm3dJBGRmFMgNFOvXg69eukyIxH59tKQkYiIAOoh\nnPG2bLFITISzz9YwWK0VKzy88YaXPn0cfvzjAJauIZAzhALhBDt0CJ55Jg7LgrFjA7Rr19YtchkD\nd9yRyD//GYdlGR54oJopU1p3svzgQfjXv7x06eLwne+c3sNqK1d6uPFGH8a4KbB3r8XUqTVt3CqJ\nhuPAggVxbN5sM3x4kEsuCbV1k05ZGjJqQk0N/PWvcfzpT/Hs3Nm6Q8RAAK65Jolf/zqRadMSue66\nJILBGDe0ld5/3+af/4wDwBiLhx5KoLKy5cvZv9/i8suTueUWH1dckcy8eXExbunJtWKFty4Mav+W\n09sf/hDP/fcnMn9+PGPH+li92tPWTTplKRCOUFEBS5d6efVVDz//eSL335/IzJkJ/OAHSZSUtDwU\nPv/c5qOP6je+4mIP27efGuMPRw+DWFbDac2Rn+9lx476zejPfz69P7TXr1/k0eO557bN0aQx8Otf\nJ/Cd7yRzzTU+du9uu+3GcWDnTovDh90DpdPNkaHuOBarVikQmqJACKuudo/mb7/dx09+kkRBQf1G\ntG+fTWFhyzaiykq3h2Hb9WPzycmGjIxTY6w+L8/hxhvdISLbNvz+99UkJrZ8Oe3bR65PauqpsX4t\nsXOnxZVXJtG9eztee83LtGlVDB0aZPz4Gn7/++o2adOzz7qfiN+zx2b1ai9Tp7bixWmGYBDuuSeB\nAQOSueEGX4MDn8pKuPZaH4MHtyMnpx1nn53CuHE+qmNQllAInnwynokTE/nHP05cT+y88yKHMfv1\ni82w5uefWxQUeNmz59Q4yIsF9YfD1q3zsG7dkTv9+hfZtg09erRsI3rggQQWLao/Wu7WzeFPf6oi\nNTXalsbO7NlV3HuvGwSZma3bkY8ZE+TNNwM895yXTp0Mjz1WFeNWnnjTpiXwwQfua//ii3Hk5YV4\n8cVWjJ+FHT4MK1d68fsNF17Yuh7Grl2Rx2pfftn8nU5JicUDDySwa5fFd78b4uqrg/Tq5fA//xPH\n4cPu52rmzo1n0yablBTDSy+5w3x799o88EACc+bUv4b/+79xvPuuu5uoHUpbtszL00/H8bOfNf+c\nkzHw+OPxvPOOh5EjHe67DxYtiuPBB933yD//GUdiYiWjR8d+TPWPf6zC5zNs2WIzYkSQq66K/jmW\nL/fwk5/4qKmxSEszvPhiBX37Nm8fkZ/v5f77E3AcmD69mh/96Pjt2bDB5k9/iseyYOrUGvr3j3YN\nGvetDITSUouvvrI45xyH+GaOYHToYLAsU7fR+3yGwYNDlJVZTJxYw7nntiwQjhwqArj66gDf/37z\ndw6BAPy//5fA8uVeEhMhLy/IhAkB8vIatmP9epvt222GDg3RoUPTO/atW20GD46c1q1bdEf0tg1P\nPVXFE0+A5zTqib/7rodf/CKRw4chKSly3v79DXe+X3/t1q9vX/D73WnGwO7dboh4vdR9gn3kyCQ+\n+cQtxi9+Uc2vf+2Os5SVwfTpCWzbZjNmTJCf/KTxHer69TbLl3vweAyhkNuW2t5cc0yalMjKle5b\n+733vDzxRAI5OSE+/9xt06xZhvLyxgNm167I6YEmnvbQoZYdFc+fH8eMGQmAG5j33QdFRZEbzJo1\nnhMSCO3awcyZse3pzZ0bT02NW4NvvrFYsCCORx45/nPs3OlezFH72MmTE5kzx2H27Eref99Lu3aG\nq68O1r2XDh50t7Prr/exf797kPDeex4++KCybjuMpdMiEHbtsli61N2zf/21xdatNiUlFhdcEAoX\nx8PQoSGGDg2xYoWHCRN8VFRY9OsX4oUXKkhLi1zes8962bDBw/e/H+Sii0Js3myTlWX44x+refjh\nBHw+90j3ssuavwPfts1i/XoP550XokcPw8UXB+uOOi3L8L3vNVxWIOAezWVmGrxHvRL//d9xLF5c\nn2ZffBHPSy/FMW9eJYcOWQwaFGLRongWL/aGv1vJwu93mDWrmpwch3vuSeTAATfMfvjDAPffn8jW\nrXEUF8PMmfHcc8+xj+Tz8718+qnNJZcEGw2hozUWBps22dx1VyJ791oMGRJi7NjGQzEQcC/13LPH\n4quvbF591UsgAD/9aYDrrguwbZvNxx/bWBZcdVWwwQ68VkGBl7//PY6MDPeqqaaG54yBn/40kQMH\n3DfYgQNuL9BxLFJTDddd5+6UiottVqzwsmKFzQcfeKmpsbj0Uoc334QtW2yuusrHV1/VH8kvXRrH\npEk1dWEAMGdOPNOm1WBZcO+9ibzwgns0vmqVl/374T//M/Ky1qoquOGG+jd/fLzhqaequPrqYF3b\nS0st0tMNdhMDvp980nBGbRgADcLgyOApLbWYMyeOa64J0qmT4frrAyxZEsf69R7AABadOztcdVWA\n9ettzj7boV07+OADG7/fcM45hs8/t6iutujb1+H5591zTB9+2LBNR9YJYMcOiyFDkjnnHIfHHqtq\nste6YYPNli02Z53lkJ7utr1zZ4eUFHe+47j7ifR0E1HbtWttZs1KwOOBadOqGTCgfrt+8013hODC\nC0ON9ur+9S8P77zjoX//EFdeGWpQ+7/9LY7OnQ133934SZbXX/fw1FPxvP++h0DgyPq7+41Ro5Lr\nXpdXXw3wX/9Vxa23JlJQEIfPZ6isrH9MSYnN3r0WXbo0+lRRsYwxp9ygb0nJobrfy8vh+99Ppn37\nw3z0URqDBh1k3Tp3D+/3O3z9tfvKWJZh6NAg69d7I45eLr88wNy5VVRWWpSVud3d3/2ufjz2rLMc\ndu2y8fkMc+dWMnx45MZw6BC8846XzEyHwYMdHMcdd92716JjR0NhoYc//CGejz7y4DgWPp9h/vxK\nAgF3Izp82OKKK4J897shVq3y8OyzcWzaZHPwoEV5uUUgYJGT4/D3v1ewe7fNJ5/YlJdbfPyxTX5+\nY1fsuG9Kr9cQDDZ1lObe5+i/Bw2C4mLIzYXzzqvh8suDDB0aYvdum+JiD36/Q6dOhrVrPXU1sm3D\nzTcHyMw0dOnicO65Dhs22Lz8chw7d1oEgxbZ2Q5TptTw5Zc2mzfbDB0apF+/ENddl8QXX0S+6UeN\nCjBrVjXvv2/z3HNxbNtmsXOnXfc6Hu3o9ezdO8QPfxhg3ToPmze7O6E776zhww9tHnssAcdx75ua\nWt877NHDobra3dmPGBGgoCCOd97xRNSoX78gY8YE2bfPJi3N3dk+8khCg/bU1vCSSxzeeqthm4/c\nJmvb8eqrFWRmGv7935P58svIx7RrZzjrLIcuXRwcxyIry+HZZyO7tQ89VEVlpXtgVFAQx759Nt27\nOzz6aCUffeThww89VFVBZaVby+3brYgrpY7cBhr+7v7t95uIdtu2YdasasaNC1BR4faODh+Gw4fd\nbfzOO33s2mVjWQaPB4JBC8syDBsWquudpKSYuvdiXJyp2xEeuR2uW+e+xvHxUFFR36asLIeHH67m\nkkuCfPmlTWWlG9zLlsWxYEFc3etcKyXFcPPNNWzY4B4gHj5s0atXiOeeq6R9e0NFBVxwQTvKytzH\nxccb0tMNNTWQlxfitdfc95plGR58sIrSUpv4eBgwIERlJfz857665+zVK8SmTbUBSUQtV6w4zJYt\nHmbPjic+HiZPriYQgEmTfHWh2xxjx9bwzDNND2+4ByY2X39dTjBYH2wZGSnNfo5GmVPQvn1ldT9v\nvnnYgDEDBnxjgPCtadGPx+MYy3IMGJOeHjrmff/jP2rMRRcFjNfrmIEDgyY7u/7+I0fWmISE+mW5\nt06DZdTOB2P++MdKU1R0yHTufOznPfIx7k/jy47mZ9Agt76DBtVPy8wMGa/XOU5bjv9z9DKO3fbY\nr1vzf5p+Xp+vfp7H0/j9Gqvh8Z+rpetbf9+0tGNvN8f6SUlxl2PbjvH5HJOaGjLJye60+PjmvF6O\n8Xgc4/U65r77qsy+fWVmypSqVrfHth3TtWuokRo2VRvH2HZ020nt49u3b0kdo3k/NP++Dd8zzf+p\nreGBA4cj9pfROuWGjIwxlJcfwgr39Tp1ggsvDJKRcYiPPoIePcqifo6zz2563u7d7m3fvm7XMzUV\nBgxwp+3caejdu2Vjp2+84RAMBsnIcMjIaE1rDd27O2zfHv0AfU6OO46dk+Ne4VErKyvqRZ8xGqth\n7XBT8xkij9Ab17dvkJISm5KS1l0M2Lt3iHPOcXj55eg/G/LKK5CbW8nbbxsGDIhmPN6Qk2PV1dCY\nltau9Xr2PFavunWO7Pm0xNlnh7j77hp+85sEDh9u+etbux2Wl1cSCtX3EMrKDCkpKXX7z5Y65YaM\nysrKSDt60F9ERJrlm2++IbWVlzOecoFgjGH79j0NEq6iopy+ff+Nxx//gg4dkti82R3HTk423Hhj\ngJUrPRw44I7r9+vn0Lt3iL/9LZ6dO23c8VT3qCwvL8i0adXMmJHA2rVuBykjw6G01MJxLNLTHcrL\n3ZNicXGGgQODrFvnJSXF0K2bw4cfNuxUJSYahg8P8vXX8PbbtUdj7rjsr35VTdeuDr/7XSKffuoh\nOdkhNzdEZqYhM9Ph3Xe9Ecv0et0TzFVV7voPHx5g8uQaSkosVq70kppqSEoy/PnPCQSDMG5cTd1l\ndI4D27fbfPqpTU6Ow1lnOXz4oYekJEPnzoZOnSz8fh/TpoVYu9ampgaGDg0SDMKHH3qproYjx54H\nDw6yZYtNVZVFQoJ7ZUoo5J6/SE019O8fwhhIS4NLL3XH5gsLPRw8WH/E4/c7/O531SQlGX7720R2\n766f166dQ7duDl26GL73vSCDB9cf6Xz5pc1bb3lISDBccEEIv9/w4IOJfPaZB6/XkJXlkJ7usGWL\nl4oKC4/HHRO+5poAhw9bbN1qkZzsLuvgQXd8fe9eu+5E8uDB7lVbgQC8846HDh0MPXs6PPOM+/qN\nGxegVy8HY+o/sOeeT4rnppviePDBEMY4jB9fg23Db3+byKZNHtq3d/jhDwPk5YWYPTuezz7zRBwB\np6QYpk2r4txzHT75xP18y+ef2xHnAGrPRc2YURW+wCGBL7+08HqhUycHv9+wZo3bzrg4Q/v2Dh07\nOgSDNqmpDued5zBwoENOjkNxsc1DDyUSDLo16tDBsG+fjcdjmDKlmpwc98Rr9+5O+P3nngCdMych\n4mh6wIAg06dX8+ST8Sxf7j73lVcGGDAgRHm5RU2NYcEC9xxOXJzhqqsCfPWVRadOhq1b3ZPyAEOG\nBAgGPTz+uM2SJQF+8IMaQiFYvdpDMOiO5y9dGsdrr3k5dOjII2f3/evxuJ/lue++alJT3XMAH3/s\nYd68eKqqLDIzHfLyguzYYbNhQ/37yrIMkyfXcMUVQZ55Jo516zwkJjrs3u2p64FlZ4fYts3ticfH\nG4wh4ug/KyvE3r3u/G7dQvzmN9VkZhoOHYJf/Sox3It325mU5GBZVvhEsTvt8ssD3HlnTcSJ7u3b\nbR59NJ79+22+970gGzfafPmlu40nJxtuuCHAyJFBKivdq+EqK2HnTpuzzoKzz/ZRVhbZQ/D7k79d\nPQSIPKlcq6LiMNnZXdi2bTdJSe4XAjkOTV5pUWvDBpv33vPQpYv79dU9ergnDINBeOklLzU17pUr\n27fb7Nplcf75IQ4csPjwQw/9+oU45xxT9zzGuDuPykrIyjLs2eO+4BdfHCIjw92ACgq8HDoEw4cH\nG1wWdvCgOwR1ZJsdB15+2Ut5OVxySbBuWGntWnejaM4VPs3l9dr4/cl1J6JCocirgyoq4IsvLN59\n1/1MwdGXAH7zjXv1RrduTV/hArBxozuzfXt3B5QQPjdrjFuDwkIPVVXuyXafr/ntD4Xc6/E7dDB1\n3wl16JD7ifDu3Z3jXoYXCLiXmyYnt76uR9fwSI29vuBeJlhd7X6Fut9fX4/GfPqpzWef2eTmhuja\ntem35vr1NmVl7vZ6vEurP/rIvWggNzfEv/2be5CQmenQs2fTy//4Y5v16208HveS7IsvDhEXVz8v\nLo4GXwdfVGSzYYOHCy4IRVymHQy6dU9MNFxwgXPMGtaqqIDXX/fi8bhXLqWnO3i9kJ7eeHurqtwT\n3h07uutU+74qLYWMDENOjqFPn4bPVVnpti0tzTB4sMOKFe62edllQRISYP9+yM+PIz4errsuwGef\n2ZSWup/vOHLbraqCzZvdE++VlRZ5eSG8Xtixw6ZDBwePp+m2H8lx4MCBY19FBk1vh9GeVD6tA0Fa\npjlvRDk21TB6qmH0zqhAaEztuYVoxsdERKRpp00gGGM4dOhQVONjIiLStNMmEERE5MTSt52KiAig\nQBARkTAFgoiIAAoEEREJUyCIiAigQPjWsm0bn89HUlJS3e1dd90FwPLlyxkyZAhpaWn079+fJUuW\ntHFrTw3Lli0jKyuLm266qcG849Vs9uzZ9OnTh/bt2zNs2DCKi4tPVrNPOU3VceXKldi2TVJSUsR2\n+fzzz9coVYZRAAAEZUlEQVTdR3WEHTt2cO2119KxY0c6d+7MhAkTKCtzv9TzhG+HUX9fqpySbNs2\nO3bsaDB9z549pl27dmbhwoWmurravPHGGyYpKckUFRW1QStPHbNmzTJ9+vQxF198sRk7dmzEvOPV\n7P/+7/9Menq6KSwsNFVVVWbmzJmmc+fOpqKioi1WpU0dq45vvfWW6dGjR5OPVR1dAwYMMD/72c9M\nRUWF2bVrlzn//PPNbbfddlK2Q/UQvqWMMZhGPmKyePFievfuzS233EJ8fDyXXXYZo0ePZt68eW3Q\nylOHz+dj7dq15OTkNJh3vJrNnTuXCRMmkJeXR0JCAlOnTsWyLPLz80/2arS5Y9XxeFRH95tKzz//\nfGbMmIHP56NLly7ccsstvP322ydlO1QgfIvdd999dO/eHb/fz8SJEykvL6eoqIjc3NyI++Xm5lJY\nWNhGrTw1TJkyhZSUxr8H5ng1O3q+ZVkMHDjwjKzpseoI7lfQXHvttWRkZNC1a1cee+yxunmqI6Sl\npTFv3jwyjvjnKTt37uSss846KduhAuFbaujQoVx55ZVs2bKFNWvWsGbNGiZNmkRpaSn+o74WND09\nnf3797dRS099x6uZato8qampDBgwgLvvvps9e/Ywf/58pk+fzsKFCwHVsTHvv/8+Tz31FNOmTTsp\n26EC4Vtq9erVTJgwgbi4OHr37s3DDz/MkiVLCAaDjQ4lybEdr2aq6fENGjSI5cuXc9FFF+H1erni\niiuYOHEiCxYsqLuP6lhv9erVDB8+nJkzZ3LppZcCJ347VCCcIbKzswmFQti2TWlpacS80tJSMjMz\n26hlp76MjIxj1ux486Vp2dnZ7A7/31rVsV5+fj4jR45k9uzZTJ48GTg526EC4Vvogw8+4N57742Y\ntnHjRhITExkxYgTvv/9+xLzCwkKGDBlyMpt4WsnLy6OoqChi2pE1O3q+4zgUFxerpkd57rnnmDNn\nTsS0jRs30rNnT0B1rPXuu+8yfvx4nn/+eW6++ea66SdlO4z1JVPS9nbt2mVSUlLMzJkzTXV1tfns\ns89Mv379zC9/+Uuzb98+k5aWZv7617+aqqoq8/LLL5vk5GSzYcOGtm72KWH8+PENLpc8Xs1effVV\n4/f7zZo1a0xFRYWZPn266d69u6mqqmqLVTglNFbHF1980SQnJ5vXX3/dBAIB89prr5mUlBTzwgsv\nGGNUR2OMCQaDpm/fvuYvf/lLg3knYztUIHxLrVq1ylx44YUmJSXFZGRkmKlTp5rq6uq6eQMHDjSJ\niYmmT58+dW/IM1liYqLx+XzG6/Uar9db93et49Vszpw5plu3bsbn85lhw4aZjz/++GSvwinheHX8\ny1/+Ynr37m2SkpJMz549zYIFCyIef6bXcdWqVca2bePz+epqV3u7Y8eOE74d6v8hiIgIoHMIIiIS\npkAQERFAgSAiImEKBBERARQIIiISpkAQERFAgSAiImEKBBERARQIIiISpkAQERFAgSAiImH/Hwa7\nH1+6+DOkAAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# (XcとSxの逆行列の積)の各要素とXcの各要素を掛け合わせる\n", "# ここでのトリックは、numpyの*が要素毎の積であることを利用\n", "# numpy()関数を使ってマトリックスをnumpyのマトリックスに変換する\n", "a_prev = (Xc*Sx.inverse()).numpy() * Xc.numpy()\n", "# 行単位の和を求める\n", "a = [sum(x) for x in a_prev]\n", "# リストプロット\n", "list_plot(a, figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

Pandaを使った異常度の計算

\n", "\t

\n", "\t\tpandaとnumpyの機能を使って異常値を計算します。ここでは、$\\Sigma^{-1} Xc$がsolveの解であることを\n", "\t\t使っています。また、次数の和には、pythonのsumを使っています。\n", "\t

\n", "\t

\n", "\t\t出力結果では、1つだけ突出した異常データが見つかります。\n", "\t

\n", "\t

\n", "\t\t全体的な処理から以降ではPandaを使った分散行列の方式を使います。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEBCAYAAAB4wNK4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt8FPW9//HXzG4um5CEDSTcBALxAIJQCEGKVXq8cgoI\naq0VLQVarQhYPSo/HpbalmqLqH2o6KOHQylQK3js0aIHRfECImKBkACCeAGVO0gIYCDX3Z3v74/Z\nJCxJIMkuBOT9fDzy2GRmd/Y7n52d93y/M7uxjDEGERE579nN3QARETk7KBBERARQIIiISJgCQURE\nAAWCiIiEKRBERARQIIiISNg5EwjGGIqLi9HHJkRETg9vczegLoWFR2tNKyk5SpcuHfjqqz0kJ6c0\nQ6vOfbZtkZ6ezKFDJTiOgrUpVMPoqYbRq6+GGRnR7RvPmR6CZVkRt9J4tm1hWRa2rRo2lWoYPdUw\neqerhudMIMTCM8/E853vJHPllUls2nRerbqIyCmdN3vF1as9PPxwAvv22Wze7OEXv/A1d5NERM4q\n500g7N1rnfRvEZHz3XkTCIMHh2jf3qn++5ZbAs3YGhGRs89ZeZXR6dC6tWHp0lLeeMNLq1aG4cOD\nzd0kEZGzynkTCABt2hjGjlXPQESkLufNkJGIiJycAkFERAAFgoiIhCkQREQEUCCIiEiYAkFERAAF\ngoiIhDU6EJYuXUrbtm259dZbI6avWLEC27ZJSkoiKSkJn89HUlISL7/8cvV9Zs6cSY8ePWjZsiWD\nBw+moKAg+jUQEZGYaNQH0x5//HHmzp1Lt27d6pyflZXFl19+Wee8xYsXM23aNJYuXUrv3r15+umn\nGT58OF988QU+n75oTkSkuTWqh+Dz+Vi7di3Z2dmNfqLZs2czbtw4cnNzSUhIYPLkyViWxeLFixu9\nLBERib1GBcKkSZNISan/P/IUFxdz4403kpGRQceOHXnyySer5+Xn55OTk1P9t2VZ9O3bl7y8vCY0\nW0REYi1m32WUmppKnz59uO+++/jHP/7B8uXL+dGPfoTf72fs2LEUFRXh9/sjHpOens7BgwdrLcu2\na/8nII/Hrr71enUuvCmOr6E0jWoYPdUweqerhjELhH79+rFs2bLqv6+55hrGjx/PvHnzGDt2LADG\nNOz/p6anJ9f6V5keTwiA1FQfqanJsWn0eSo1VedsoqUaRk81jF6sa3hav+00Kyur+iqjjIwMioqK\nIuYXFRXRu3fvWo87dKikVg+hpKQMgOLiMkIhz2lq8bebx2OTmuoL19A59QOkFtUweqph9Oqrod8f\n3cFyzALhpZde4uDBg4wfP7562pYtW+jatSsAubm55OfnM3r0aAAcx6GgoIDbb7+91rIcx+A4kb2J\nqpUOhRyCQW1E0VANo6caRk81jF6saxizAaj4+HgeeOAB3nnnHYLBIG+//Tbz589nwoQJANx11108\n99xzrFmzhrKyMh555BESExMZNmxYrJogIiJRaFQPwefzYVkWgYD7T2YWLVqEZVmUlpYyYsQInnrq\nKSZNmsSuXbto27YtM2fOZOTIkQAMGTKE6dOnc/PNN1NYWMiAAQNYsmQJCQkJsV8rERFpNMs09Ezv\nGVRYeLTWtNLSY2RltWf79r0kJbVohlad+7xeG78/mcOHS9RVbyLVMHqqYfTqq2FGRv0fC2gIXfcl\nIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQ\nREQkTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIU\nCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKA\nAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQREQk\nTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIUCCIiAigQREQkTIEgIiKAAkFERMIaHQhL\nly6lbdu23HrrrbXmLVu2jIEDB5KWlkbv3r1ZuHBhxPyZM2fSo0cPWrZsyeDBgykoKGh6y0VEJKYa\nFQiPP/449957L926das1b//+/YwcOZIJEyZQWFjIU089xR133FG901+8eDHTpk3j+eef5+uvv2b4\n8OEMHz6csrKy2KyJiIhEpVGB4PP5WLt2LdnZ2bXmLViwgO7duzNmzBji4+O56qqrGDFiBHPmzAFg\n9uzZjBs3jtzcXBISEpg8eTKWZbF48eLYrImIiESlUYEwadIkUlJS6pyXn59PTk5OxLScnBzy8vLq\nnG9ZFn379q2eLyIizcsbqwUVFRXRsWPHiGnp6ekcPHiwer7f7693/vFs28K2rYhpHo9dfev16lx4\nUxxfQ2ka1TB6qmH0TlcNYxYIAMaYqOZXSU9PxrJODIQQAKmpPlJTk5vWQAHcGkp0VMPoqYbRi3UN\nYxYIGRkZFBUVRUwrKioiMzPzpPN79+5da1mHDpXU6iGUlLgnn4uLywiFPLFq9nnF47FJTfWFa+g0\nd3POSaph9FTD6NVXQ78/uoPlmAVCbm4u8+fPj5iWl5fHwIEDq+fn5+czevRoABzHoaCggNtvv73W\nshzH4DiRvYmqlQ6FHIJBbUTRUA2jpxpGTzWMXqxrGLMBqNtuu43t27czd+5cKioqWLJkCW+88QZ3\n3nknAHfddRfPPfcca9asoaysjEceeYTExESGDRsWqyaIiEgUGtVD8Pl8WJZFIBAAYNGiRViWRWlp\nKRkZGbz22mvcfffdTJw4kaysLBYsWECvXr0AGDJkCNOnT+fmm2+msLCQAQMGsGTJEhISEmK/ViIi\n0miWaeiZ3jOosPBorWmlpcfIymrP9u17SUpq0QytOvd5vTZ+fzKHD5eoq95EqmH0VMPo1VfDjIy6\nPxbQULruS0REAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiE\nKRBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCQURE\nAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiI\nSJgCQUREAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiEKRBE\nRARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCQUREgBgH\ngm3b+Hw+kpKSqm/vueceAJYtW8bAgQNJS0ujd+/eLFy4MJZPLSIiUfLGcmGWZfH555/TsWPHiOn7\n9+9n5MiRPPvss4waNYqVK1cyYsQIevToQU5OTiybICIiTRTTHoIxBmNMrekLFiyge/fujBkzhvj4\neK666ipGjBjBnDlzYvn0IiIShZifQ5gyZQqdO3fG7/czfvx4SkpKyM/Pr9UTyMnJIS8vL9ZPLyIi\nTRTTQBg0aBDXXnst27ZtY/Xq1axevZoJEyZQVFSE3++PuG96ejoHDx6M5dOLiEgUYnoOYdWqVdW/\nd+/enUcffZTrrruOwYMH1zmUVB/btrBtK2Kax2NX33q9ujiqKY6voTSNahg91TB6p6uGMQ2EE2Vl\nZREKhbBtm6Kiooh5RUVFZGZm1vm49PRkLOvEQAgBkJrqIzU1+fQ0+DyRmupr7iac81TD6KmG0Yt1\nDWMWCBs2bOD555/niSeeqJ62ZcsWEhMTGTp0KPPnz4+4f15eHgMHDqxzWYcOldTqIZSUlAFQXFxG\nKOSJVbPPKx6PTWqqL1xDp7mbc05SDaOnGkavvhr6/dEdLMcsEDIzM5k9ezaZmZnce++9bN++nd/8\n5jfceeed/OQnP2HatGnMnTuX2267jXfffZc33niDNWvW1LksxzE4TuQQU9VKh0IOwaA2omiohtFT\nDaOnGkYv1jWM2QBU+/btWbJkCa+++iqtW7fmsssuY+jQocyYMYOMjAxee+01nnnmGVq2bMn999/P\nggUL6NWrV6yeXkREohTTcwiXXXZZxInlE+etX78+lk8nIiIxpNP8IiICKBBERCRMgSAiIoACQURE\nwhQIIiICKBBERCRMgSAiIoACQUREwhQIIiICKBBERCRMgSAiIoACQUREwhQIIiICKBBERCRMgSAi\nIoACQUREwhQIIiICKBBERCRMgSAiIoACQUREwhQIIiICKBBERCRMgSAiIoACQUREwhQIIiICKBBE\nRCRMgSAiIoACQUREwhQIIiICKBBERCRMgSAiIoACQUREwhQIIiICKBBERCRMgSAiIoACQUREwhQI\nIiICKBBERCRMgSAiIoACQUREwhQIjeA4zd0CEZHTR4HQADt2WAwenES7di245RYfpaXN3SIRkdhT\nIDTAb3+bwKefejDGYtkyL7Nnxzd3k0REYk6B0ADffGNF/H3kiFXPPUVEzl0KhAa4444AXq8BwO83\n3HproJlbJCISe97mbsC5YOjQIO+9V8rWrTb9+4do29Y0d5NERGJOgdBA3bo5dOumy4xE5NtLQ0Yi\nIgKoh3De27bNIjERLrhAw2BVli/38M47Xnr0cPjJTwJYuoZAzhMKhNPs6FF44YU4LAtGjQrQokVz\nt8hlDNx1VyL//GcclmV46KEKJk1q2snyI0fgX//y0r69w3e+c24Pq61Y4eGWW3wY46bA/v0WkydX\nNnOrJBqOA/PmxbF1q82QIUGuuCLU3E06a2nIqB6VlfDXv8bxpz/Fs2tX0w4RAwG44YYkfv3rRKZO\nTeSmm5IIBmPc0CZat87mn/+MA8AYi0ceSaCsrPHLOXjQ4uqrkxkzxsc11yQzZ05cjFt6Zi1f7q0O\ng6q/5dz2hz/E8+CDicydG8+oUT5WrfI0d5POWgqE45SWwqJFXt5808MvfpHIgw8mMmNGAj/4QRKF\nhY0PhS++sPnoo5qNr6DAw44dZ8f4w4nDIJZVe1pDLF7sZefOms3oz38+tz+016tX5NHjRRc1z9Gk\nMfDrXyfwne8kc8MNPvbubb7txnFg1y6LY8fcA6VzzfGh7jgWK1cqEOqjQAirqHCP5u+808dPf5rE\nkiU1G9GBAzZ5eY3biMrK3B6GbdeMzScnGzIyzo6x+txch1tucYeIbNvw+99XkJjY+OW0bBm5Pqmp\nZ8f6NcauXRbXXptE584teOstL1OnljNoUJCxYyv5/e8rmqVNL77ofiJ+3z6bVau8TJ7chBenAYJB\nuP/+BPr0SebHP/bVOvApK4Mbb/TRv38LsrNbcMEFKYwe7aMiBmUJheCZZ+IZPz6Rf/zj9PXELr44\nchizV6/YDGt+8YXFkiVe9u07Ow7yYkH94bD16z2sX3/8Tr/mRbZtQ5cujduIHnoogeeeqzla7tTJ\n4U9/Kic1NdqWxs7MmeU88IAbBJmZTduRjxwZ5N13A7z0kpc2bQxPPlke41aeflOnJrBhg/vav/pq\nHLm5IV59tQnjZ2HHjsGKFV78fsOllzath7FnT+Sx2u7dDd/pFBZaPPRQAnv2WHz3uyGuvz5It24O\n//M/cRw75n6uZvbseD7/3CYlxfDaa+4w3/79Ng89lMCsWTWv4f/+bxwffujuJqqG0pYu9fL883H8\n/OcNP+dkDDz1VDwffOBh2DCHKVPguefiePhh9z3yz3/GkZhYxogRsR9T/eMfy/H5DNu22QwdGuS6\n66J/jmXLPPz0pz4qKy3S0gyvvlpKz54N20csXuzlwQcTcByYNq2CH/3o1O3ZvNnmT3+Kx7Jg8uRK\neveOdg3q9q0MhKIii6+/trjwQof4Bo5gtGplsCxTvdH7fIb+/UMUF1uMH1/JRRc1LhCOHyoCuP76\nAN//fsN3DoEA/L//l8CyZV4SEyE3N8i4cQFyc2u3Y9Mmmx07bAYNCtGqVf079i+/tOnfP3Jap07R\nHdHbNjz7bDlPPw2ec6gn/uGHHn75y0SOHYOkpMh5Bw/W3vkePuzWr2dP8PvdacbA3r1uiHi9VH+C\nfdiwJD75xC3GL39Zwa9/7Y6zFBfDtGkJbN9uM3JkkJ/+tO4d6qZNNsuWefB4DKGQ25aq3lxDTJiQ\nyIoV7lt7zRovTz+dQHZ2iC++cNv02GOGkpK6A2bPnsjpgXqe9ujRxh0Vz50bx/TpCYAbmFOmQH5+\n5AazerXntARCixYwY0Zse3qzZ8dTWenW4JtvLObNi+Pxx0/9HLt2uRdzVD124sREZs1ymDmzjHXr\nvLRoYbj++mD1e+nIEXc7u/lmHwcPugcJa9Z42LChrHo7jKVzIhD27LFYtMjdsx8+bPHllzaFhRaX\nXBIKF8fDoEEhBg0KsXy5h3HjfJSWWvTqFeKVV0pJS4tc3osvetm82cP3vx/ksstCbN1q07at4Y9/\nrODRRxPw+dwj3auuavgOfPt2i02bPFx8cYguXQyXXx6sPuq0LMP3vld7WYGAezSXmWnwnvBK/Pd/\nx7FgQU2affVVPK+9FsecOWUcPWrRr1+I556LZ8ECb/i7lSz8fofHHqsgO9vh/vsTOXTIDbMf/jDA\ngw8m8uWXcRQUwIwZ8dx//8mP5Bcv9vLppzZXXBGsM4ROVFcYfP65zT33JLJ/v8XAgSFGjao7FAMB\n91LPffssvv7a5s03vQQC8LOfBbjppgDbt9t8/LGNZcF11wVr7cCrLFni5e9/jyMjw71qqr7hOWPg\nZz9L5NAh9w126JDbC3Qci9RUw003uTulggKb5cu9LF9us2GDl8pKiyuvdHj3Xdi2zea663x8/XXN\nkfyiRXFMmFBZHQYAs2bFM3VqJZYFDzyQyCuvuEfjK1d6OXgQ/vM/Iy9rLS+HH/+45s0fH2949tly\nrr8+WN32oiKL9HSDXc+A7yef1J5RFQZArTA4PniKiixmzYrjhhuCtGljuPnmAAsXxrFpkwcwgEW7\ndg7XXRdg0yabCy5waNECNmyw8fsNF15o+OILi4oKi549HV5+2T3HtHFj7TYdXyeAnTstBg5M5sIL\nHZ58srzeXuvmzTbbttl06OCQnu62vV07h5QUd77juPuJ9HQTUdu1a20eeywBjwemTq2gT5+a7frd\nd90RgksvDdXZq/vXvzx88IGH3r1DXHttqFbt//a3ONq1M9x3X90nWd5+28Ozz8azbp2HQOD4+rv7\njeHDk6tflzffDPBf/1XO7bcnsmRJHD6foays5jGFhTb791u0b1/nU0XFMsacdYO+hYVHq38vKYHv\nfz+Zli2P8dFHafTrd4T16909vN/vcPiw+8pYlmHQoCCbNnkjjl6uvjrA7NnllJVZFBe73d3f/a5m\nPLZDB4c9e2x8PsPs2WUMGRK5MRw9Ch984CUz06F/fwfHccdd9++3aN3akJfn4Q9/iOejjzw4joXP\nZ5g7t4xAwN2Ijh2zuOaaIN/9boiVKz28+GIcn39uc+SIRUmJRSBgkZ3t8Pe/l7J3r80nn9iUlFh8\n/LHN4sV1XbHjvim9XkMwWN9RmnufE//u1w8KCiAnBy6+uJKrrw4yaFCIvXttCgo8+P0ObdoY1q71\nVNfItg233RYgM9PQvr3DRRc5bN5s8/rrcezaZREMWmRlOUyaVMnu3TZbt9oMGhSkV68QN92UxFdf\nRb7phw8P8NhjFaxbZ/PSS3Fs326xa5dd/Tqe6MT17N49xA9/GGD9eg9bt7o7obvvrmTjRpsnn0zA\ncdz7pqbW9A67dHGoqHB39kOHBliyJI4PPvBE1KhXryAjRwY5cMAmLc3d2T7+eEKt9lTV8IorHN57\nr3abj98mq9rx5pulZGYa/v3fk9m9O/IxLVoYOnRwaN/ewXEs2rZ1ePHFyG7tI4+UU1bmHhgtWRLH\ngQM2nTs7PPFEGR995GHjRg/l5VBW5tZyxw4r4kqp47eB2r+7f/v9JqLdtm147LEKRo8OUFrq9o6O\nHYNjx9xt/O67fezZY2NZBo8HgkELyzIMHhyq7p2kpJjq92JcnKneER6/Ha5f777G8fFQWlrTprZt\nHR59tIIrrgiye7dNWZkb3EuXxjFvXlz161wlJcVw222VbN7sHiAeO2bRrVuIl14qo2VLQ2kpXHJJ\nC4qL3cfFxxvS0w2VlZCbG+Ktt9z3mmUZHn64nKIim/h46NMnRFkZ/OIXvurn7NYtxOefVwUkEbVc\nvvwY27Z5mDkznvh4mDixgkAAJkzwVYduQ4waVckLL9Q/vOEemNgcPlxCMFgTbBkZKQ1+jjqZs9CB\nA8XVP+++e8yAMX36fGOA8K1p1I/H4xjLcgwYk54eOul9/+M/Ks1llwWM1+uYvn2DJiur5v7DhlWa\nhISaZbm3Tq1lVM0HY/74xzKTn3/UtGt38uc9/jHuT93LjuanXz+3vv361UzLzAwZr9c5RVtO/XPi\nMk7e9tivW8N/6n9en69mnsdT9/3qquGpn6ux61tz37S0k283J/tJSXGXY9uO8fkck5oaMsnJ7rT4\n+Ia8Xo7xeBzj9TpmypRyc+BAsZk0qbzJ7bFtx3TsGKqjhvXVxjG2Hd12UvX4li0bU8do3g8Nv2/t\n90zDf6pqeOjQsYj9ZbTOuiEjYwwlJUexwn29Nm3g0kuDZGQc5aOPoEuX4qif44IL6p+3d69727On\n2/VMTYU+fdxpu3YZundv3NjpO+84BINBMjIcMjKa0lpD584OO3ZEP0Cfne2OY2dnu1d4VGnbNupF\nnzfqqmHVcFPDGSKP0OvWs2eQwkKbwsKmXQzYvXuICy90eP316D8b8sYbkJNTxvvvG/r0iWY83pCd\nbVXX0JjG1q7punY9Wa+6aY7v+TTGBReEuO++Sn7zmwSOHWv861u1HZaUlBEK1fQQiosNKSkp1fvP\nxjrrhoyKi4tJO3HQX0REGuSbb74htYmXM551gWCMYceOfbUSrrS0hJ49/42nnvqKVq2S2LrVHcdO\nTjbcckuAFSs8HDrkjuv36uXQvXuIv/0tnl27bNzxVPeoLDc3yNSpFUyfnsDatW4HKSPDoajIwnEs\n0tMdSkrck2JxcYa+fYOsX+8lJcXQqZPDxo21O1WJiYYhQ4IcPgzvv191NOaOy/7qVxV07Ojwu98l\n8umnHpKTHXJyQmRmGjIzHT780BuxTK/XPcFcXu6u/5AhASZOrKSw0GLFCi+pqYakJMOf/5xAMAij\nR1dWX0bnOLBjh82nn9pkZzt06OCwcaOHpCRDu3aGNm0s/H4fU6eGWLvWprISBg0KEgzCxo1eKirg\n+LHn/v2DbNtmU15ukZDgXpkSCrnnL1JTDb17hzAG0tLgyivdsfm8PA9HjtQc8fj9Dr/7XQVJSYbf\n/jaRvXtr5rVo4dCpk0P79obvfS9I//41Rzq7d9u8956HhATDJZeE8PsNDz+cyGefefB6DW3bOqSn\nO2zb5qW01MLjcceEb7ghwLFjFl9+aZGc7C7ryBF3fH3/frv6RHL//u5VW4EAfPCBh1atDF27Orzw\ngvv6jR4doFs3B2NqPrDnnk+K59Zb43j44RDGOIwdW4ltw29/m8jnn3to2dLhhz8MkJsbYubMeD77\nzBNxBJySYpg6tZyLLnL45BP38y1ffGFHnAOoOhc1fXp5+AKHBHbvtvB6oU0bB7/fsHq12864OEPL\nlg6tWzsEgzapqQ4XX+zQt69DdrZDQYHNI48kEgy6NWrVynDggI3HY5g0qYLsbPfEa+fOTvj9554A\nnTUrIeJouk+fINOmVfDMM/EsW+Y+97XXBujTJ0RJiUVlpWHePPccTlyc4brrAnz9tUWbNoYvv3RP\nygMMHBggGPTw1FM2CxcG+MEPKgmFYNUqD8GgO56/aFEcb73l5ejR44+c3fevx+N+lmfKlApSU91z\nAB9/7GHOnHjKyy0yMx1yc4Ps3GmzeXPN+8qyDBMnVnLNNUFeeCGO9es9JCY67N3rqe6BZWWF2L7d\n7YnHxxuMIeLov23bEPv3u/M7dQrxm99UkJlpOHoUfvWrxHAv3m1nUpKDZVnhE8XutKuvDnD33ZUR\nJ7p37LB54ol4Dh60+d73gmzZYrN7t7uNJycbfvzjAMOGBSkrc6+GKyuDXbtsOnSACy7wUVwc2UPw\n+5O/XT0EiDypXKW09BhZWe3Zvn0vSUnuFwI5DvVeaVFl82abNWs8tG/vfn11ly7uCcNgEF57zUtl\npXvlyo4dNnv2WAwYEOLQIYuNGz306hXiwgtN9fMY4+48ysqgbVvDvn3uC3755SEyMtwNaMkSL0eP\nwpAhwVqXhR054g5BHd9mx4HXX/dSUgJXXBGsHlZau9bdKBpyhU9Deb02fn9y9YmoUCjy6qDSUvjq\nK4sPP3Q/U3DiJYDffONevdGpU/1XuABs2eLObNnS3QElhM/NGuPWIC/PQ3m5e7Ld52t4+0Mh93r8\nVq1M9XdCHT3qfiK8c2fnlJfhBQLu5abJyU2v64k1PF5dry+4lwlWVLhfoe7319SjLp9+avPZZzY5\nOSE6dqz/rblpk01xsbu9nurS6o8+ci8ayMkJ8W//5h4kZGY6dO1a//I//thm0yYbj8e9JPvyy0PE\nxdXMi4uj1tfB5+fbbN7s4ZJLQhGXaQeDbt0TEw2XXOKctIZVSkvh7be9eDzulUvp6Q5eL6Sn193e\n8nL3hHfr1u46Vb2vioogI8OQnW3o0aP2c5WVuW1LSzP07++wfLm7bV51VZCEBDh4EBYvjiM+Hm66\nKcBnn9kUFbmf7zh+2y0vh61b3RPvZWUWubkhvF7YudOmVSsHj6f+th/PceDQoZNfRQb1b4fRnlQ+\npwNBGqchb0Q5OdUweqph9M6rQKhL1bmFaMbHRESkfudMIBhjOHr0aFTjYyIiUr9zJhBEROT00red\niogIoEAQEZEwBYKIiAAKBBERCVMgiIgIoED41rJtG5/PR1JSUvXtPffcA8CyZcsYOHAgaWlp9O7d\nm4ULFzZza88OS5cupW3bttx666215p2qZjNnzqRHjx60bNmSwYMHU1BQcKaafdapr44rVqzAtm2S\nkpIitsvFNjeDAAAEZ0lEQVSXX365+j6qI+zcuZMbb7yR1q1b065dO8aNG0dxsfulnqd9O4z6+1Ll\nrGTbttm5c2et6fv27TMtWrQw8+fPNxUVFeadd94xSUlJJj8/vxlaefZ47LHHTI8ePczll19uRo0a\nFTHvVDX7v//7P5Oenm7y8vJMeXm5mTFjhmnXrp0pLS1tjlVpVier43vvvWe6dOlS72NVR1efPn3M\nz3/+c1NaWmr27NljBgwYYO64444zsh2qh/AtZYzB1PERkwULFtC9e3fGjBlDfHw8V111FSNGjGDO\nnDnN0Mqzh8/nY+3atWRnZ9ead6qazZ49m3HjxpGbm0tCQgKTJ0/GsiwWL158plej2Z2sjqeiOrrf\nVDpgwACmT5+Oz+ejffv2jBkzhvfff/+MbIcKhG+xKVOm0LlzZ/x+P+PHj6ekpIT8/HxycnIi7peT\nk0NeXl4ztfLsMGnSJFJS6v4emFPV7MT5lmXRt2/f87KmJ6sjuF9Bc+ONN5KRkUHHjh158sknq+ep\njpCWlsacOXPIOO6fp+zatYsOHTqcke1QgfAtNWjQIK699lq2bdvG6tWrWb16NRMmTKCoqAj/CV8L\nmp6ezsGDB5uppWe/U9VMNW2Y1NRU+vTpw3333ce+ffuYO3cu06ZNY/78+YDqWJd169bx7LPPMnXq\n1DOyHSoQvqVWrVrFuHHjiIuLo3v37jz66KMsXLiQYDBY51CSnNypaqaanlq/fv1YtmwZl112GV6v\nl2uuuYbx48czb9686vuojjVWrVrFkCFDmDFjBldeeSVw+rdDBcJ5Iisri1AohG3bFBUVRcwrKioi\nMzOzmVp29svIyDhpzU41X+qXlZXF3vD/rVUdayxevJhhw4Yxc+ZMJk6cCJyZ7VCB8C20YcMGHnjg\ngYhpW7ZsITExkaFDh7Ju3bqIeXl5eQwcOPBMNvGckpubS35+fsS042t24nzHcSgoKFBNT/DSSy8x\na9asiGlbtmyha9eugOpY5cMPP2Ts2LG8/PLL3HbbbdXTz8h2GOtLpqT57dmzx6SkpJgZM2aYiooK\n89lnn5levXqZe++91xw4cMCkpaWZv/71r6a8vNy8/vrrJjk52WzevLm5m31WGDt2bK3LJU9Vszff\nfNP4/X6zevVqU1paaqZNm2Y6d+5sysvLm2MVzgp11fHVV181ycnJ5u233zaBQMC89dZbJiUlxbzy\nyivGGNXRGGOCwaDp2bOn+ctf/lJr3pnYDhUI31IrV640l156qUlJSTEZGRlm8uTJpqKionpe3759\nTWJiounRo0f1G/J8lpiYaHw+n/F6vcbr9Vb/XeVUNZs1a5bp1KmT8fl8ZvDgwebjjz8+06twVjhV\nHf/yl7+Y7t27m6SkJNO1a1czb968iMef73VcuXKlsW3b+Hy+6tpV3e7cufO0b4f6fwgiIgLoHIKI\niIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhImAJBREQABYKIiIQpEEREBFAgiIhI2P8H\nJFocU+ZnmD4AAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Pandaのみを使った計算\n", "# 中心化したデータの値をXcにセット\n", "Xc =(X - mx).values\n", "# 標本共分散行列を求める\n", "Sx = (X - mx).cov().values\n", "# 異常値aを求めるΣ^-1 * Xcは、solveの解を利用\n", "a_prev = np.linalg.solve(Sx, Xc.T).T * Xc\n", "# 行単位の和を求める\n", "a = [sum(x) for x in a_prev]\n", "# リストプロット\n", "list_plot(a, figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

マハラノビス=タグチ法

\n", "\t

\n", "\t\t異常値の検出では、全体に対する個々の観測値の異常を検出するのに使用しましたが、\n", "\t\t個々の観測値の変数値の異常を検出するには、マハラノビス=タグチ法を使用します。\n", "\t

\n", "\t

\n", "\t\tM変数のなかからいくつかの変数を選び、1変数当たりの異常度を以下の様に検出します。\n", "\t\t\n", "$$\n", "\tSN_q = -10 log_{10} \\left \\{ \\frac{1}{N'} \\sum_{n=1}^{N'} \\frac{1}{a_q(x'^{(n)})/M_q} \\right \\}\n", "$$\t\t\n", "\t

\n", "\t

MASSパッケージのroadデータで実験

\n", "\t

\n", "\t\tマハラノビス=タグチ法の実行例2.6をSageを使って試してみましょう。\n", "\t

\n", "\t

\n", "\t\tMASSパッケージのroadデータには、アメリカ26州の交通事故死亡者数deaths、運転者数drivers、\n", "\t\t人口密度popden、郊外地区の道路長rural、1月における1日の最高気温の平均値temp、1年度との燃料消費量\n", "\t\tfuelの6変数を記録しています。\n", "\t

\n", "\t

\n", "\t\tデータの前段階でroadをdriversで割って対数化する必要があるのですが、pandaで処理する方法が見つからず、\n", "\t\tRで前処理することにしました(残念)。\n", "\t

\n", "\t

\n", "\t\t取り込んだデータは、infoで個数と欠損値の有無を調べ、内容を調べます。\n", "\t\t$rowという変な名前が存在するので、stateに置き替えます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Int64Index: 26 entries, 0 to 25\n", "Data columns (total 6 columns):\n", "_row 26 non-null object\n", "deaths 26 non-null float64\n", "fuel 26 non-null float64\n", "popden 26 non-null float64\n", "rural 26 non-null float64\n", "temp 26 non-null float64\n", "dtypes: float64(5), object(1)\n", "memory usage: 1.4+ KB\n" ] } ], "source": [ "# マハラノビス=タグチ法\n", "# MASSパッケージのroadデータを使って計算\n", "r('library(MASS)')\n", "# roadをPandaのデータフレームに変換\n", "# road / road$driversとlogを使った対数化がPandaで処理する方法がみつからない\n", "r('X <- road / road$drivers')\n", "r('X <-log(X[, -2] + 1)')\n", "road = RDf2PandaDf(\"X\")\n", "road.info()" ] }, { "cell_type": "code", "execution_count": 24, "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", "
_rowdeathsfuelpopdenruraltemp
0Alabama1.96380.56140.34010.34910.3310
1Alaska1.59110.44700.03570.42941.3157
2Arizona2.00980.53900.12390.30940.5326
3Arkanas2.07400.59020.31450.58420.4411
4Calif1.78880.10460.09990.11680.0660
\n", "
" ], "text/plain": [ " _row deaths fuel popden rural temp\n", "0 Alabama 1.9638 0.5614 0.3401 0.3491 0.3310\n", "1 Alaska 1.5911 0.4470 0.0357 0.4294 1.3157\n", "2 Arizona 2.0098 0.5390 0.1239 0.3094 0.5326\n", "3 Arkanas 2.0740 0.5902 0.3145 0.5842 0.4411\n", "4 Calif 1.7888 0.1046 0.0999 0.1168 0.0660" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# _rowという変な名前が存在するので、headで内容をチェック、_rowが州の名前だと判明\n", "road.head()" ] }, { "cell_type": "code", "execution_count": 25, "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", "
statedeathsfuelpopdenruraltemp
4Calif1.78880.10460.09990.11680.066
\n", "
" ], "text/plain": [ " state deaths fuel popden rural temp\n", "4 Calif 1.7888 0.1046 0.0999 0.1168 0.066" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 州のカラムにstateをセットする\n", "road.columns = ['state', 'deaths', 'fuel', 'popden', 'rural', 'temp']\n", "# Calif州のデータをピックアップする\n", "road[road.state == 'Calif']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

1変数当たりの異常度とCalifのSN比を求める

\n", "\t

\n", "\t\t異常度の計算と同様に異常度を求め、次元数で割って1変数当たりの異常度を求めます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# pandaデータフレームで、平均からのずれXcを計算する\n", "X = road[['deaths', 'fuel', 'popden', 'rural', 'temp']]\n", "mx = X.mean()\n", "Xc = X - mx\n", "Xc.shape[1]" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAECCAYAAAD+VKAWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAGPhJREFUeJzt3XuUk/Wdx/FPLpOZZGYSZpjRgamD3ITKStVKUY7lLKho\ndUE7x6MW6AF0FUG0KrR6VnYPeO1QW88q9LDluiJseyquLHp2FQSEAsql6ooUK3bXYZFbh3aCc0/y\n2z/yMBjCJTOT5Ekm79c/4cnk8uWbJ/k8v+fqMMYYAQByntPuAgAAmYFAAABIIhAAABYCAQAgiUAA\nAFgIBACAJAIBAGAhEAAAktIYCMYYBYNBcRwcAGQmdype9NixE3H3NTScUN++lfqf/zmowsLiVLxt\n1nE6HSotLdTx4w2KRAhKiZ6cjn7EoyexzteP8vLEf2/TNkJwOBwxt91JW5tUV9fx/5fT6ZDD4ZDT\n2f160ln0JBb9iEdPYiWzH2xD6KLf/96poUML9c1vFun73/eqsdHuigCgcwiELpo9u0B1ddE2bt3q\n1iuv5NlcEQB0DoHQRS0tsdPNzQxjAWQnAqGLZs1qlccT3ZDTt29E48e32VwRAHROSvYyyiXf+15I\n773XoIMHnRoyJKyiIrsrAoDOIRCS4BvfMPrGN8J2lwEAXcIqIwCAJAIBAGAhEAAAkggEnMef/uTQ\nrl1OtbHzFNDtEQg4qyVL8nTNNYW6+eZCVVd74465ANC9EAg4q2eeyZcx0QPt3n/frXXr2CkN6M4I\nBJxVQUHsmRPz8zmzJNCdEQg4q+efb5HXGw2B6uo2XX89x1oA3RnrAHBWN98c0r59X6mx0aGePRkd\nAN1dp0cIjzzyiJxOBhjdndcrwgDIEZ36Rf/www+1YsWKbnmxGwDIVR0OBGOMpk2bppkzZ6aiHgCA\nTTocCAsXLpTX69X48eNTUQ8AwCYd2qh85MgRzZkzR5s3bz7n45zO+Ot7ulzO9lu3m20PUmxPEEVP\nYtGPePQkVjL74TDGJLzFcOLEiaqqqtKzzz6rL774Qv369VM4HL8rojEmbvtCMBhUIBBQfX29/H5/\nlwsHACRXwiOEd955R9u2bdOiRYskRX/0z+b48Ya4EUJDQ5MkKRhsUjjs6kyt3Y7L5ZTf77V6ErG7\nnIxAT2LRj3j0JNb5+lFSUpjwayUcCCtXrtTRo0dVVVUlSYpEIjLG6IILLtD8+fN1xx13tD82EjGK\nRGID42Sh4XBEoRAf4tfRk3j0JBb9iEdPYiWjHwkHwgsvvKCnn366ffrAgQO65ppr9NFHH6mkpKRL\nRQAA7JdwIAQCAQUCgfbptrY2ORwO9erVKyWFAQDSq9Obpfv06XPGDcoAgOzEflsAAEkEAgDAQiAA\nACQRCAAAC4EAAJBEIAAALAQCAEASgQAAsBAIAABJBAIAwEIgAAAkEQgAAAuBAACQRCAAACwEAgBA\nEoEAALAQCAAASQQCAMCS8YFgjNTaancVAND9ZXQgbNrk0qBBRaqqKtLjj+fbXQ6A04RC0ksveTRz\nZr7eecdldznooowOhAcfLNBf/+pQJOLQ0qUebdrEDAdkkn/6p3w99VS+VqzwaMIEr7Zt4zuazTI6\nEE6ccJxzGoC9fve7UwEQiTgIhCyX0YHw0EOnNh78zd+ENXp0yMZqAJxu6NBIzPS3vhW2qRIkg9vu\nAs7l0UdbNXp0SMePO3T11WH5fHZXBODramqaFQgY/elPTv3d34V0ww0EQjbL6ECQpMsvj5z/QQBs\nUVgoPfNMi91lIEkyepURACB9CAQAgCQCAQBgIRAAAJIIBACAhUAAAEgiEAAAFgIBACCJQAAAWAgE\nAIAkAgEAYCEQAACSOhEIH330ka6//nr16NFDvXr10l133aUjR46kojYAQBp1KBBaW1t14403avTo\n0Tp27Jj27NmjI0eOaPr06amqDwCQJh0KhMbGRj377LN6/PHHlZeXp549e6q6ulp79uxJVX0AgDTp\n0PUQevToobvvvrt9+tNPP9Xy5ct11113Jb0wAEB6deoCObW1tRo4cKDC4bDuu+8+zZkzJ+bvTqdD\nTmfs9Y9dLmf7rdvNtmwptieIoiex6Ec8ehIrmf1wGGNMZ5/8+eef67777lNFRYVWrlzZfr8xRg5H\nbCAEg0EFAgHV19fL7/d3vmIAQEp0KRAk6b333tOIESN07Ngx9ezZU5JUV/dV3AihoeErXXRRhQ4c\nOKzCwqKuvGW34XI55fd7FQw2KRzmUqESPTkd/YhHT2Kdrx8lJYUJv1aHVhlt3LhR06ZN0759+9rv\nczgccjgc8ng87fdFIkaRSGzOnCw0HI4oFOJD/Dp6Eo+exKIf8ehJrGT0o0Mrnb797W8rGAzqscce\nU1NTk44dO6a5c+dq5MiRKi4u7lIhAAB7dSgQ/H6/1q1bpx07dqi8vFyXXXaZSkpKtGrVqlTVBwBI\nkw7vZTRkyBBt3LgxFbUAAGzEflsAAEkEAgDAQiAAACQRCAAAC4EAAJBEIAAALAQCAEASgQAAsBAI\nAABJBAIAwEIgAAAkEQgAAAuBAACQRCAAACwEAgBAEoEAALAQCAAASQQCAMBCIAAAJBEIAAALgQAA\nkEQgAAAsBAIAQBKBAACwEAgAAEkEAgDAQiAAACQRCAAAC4EAAJBEIAAALAQCAGSxP/4x+jNeV+fo\n8msRCACQpVatcmv69AJJ0tSpBfrii66FAoEAAFlq0SKPjImGwF/+4tRrr+V16fUIBADIUqWl5pzT\nHUUgAECWmjevWX37hiVJo0e3acKEti69HoEAAFmqf3+jxYubJUlPPNEqt7trr9fhQKitrVV1dbXK\nysrUq1cvTZkyRcFgsGtVAABs1+FAGDt2rEpLS3XgwAHt3r1bn3zyiWbNmpWK2gAAadShQKivr9ew\nYcP03HPPyev1qnfv3po0aZI2b96cqvoAAGnSoTVOgUBAixcvjrmvtrZWlZWVSS0KAJB+XdoEsWvX\nLs2fP19vvPFGsuoBANik04GwdetWjRs3TvPmzdOoUaNi/uZ0OuR0xh4x53I522/dbnZukmJ7gih6\nEot+xKMnsZLZD4cxpsNHMqxdu1Y//OEPtWDBAk2YMCHu78YYORyxgRAMBhUIBFRfXy+/39/5igEA\nKdHhEcK2bds0efJkrV69Wtddd90ZH3P8eEPcCKGhoUmSFAw2KRx2daLU7sflcsrv91o9idhdTkag\nJ7HoRzx6Eut8/SgpKUz4tToUCOFwWPfee69qamrOGgaSFIkYRSKxA4+ThYbDEYVCfIhfR0/i0ZNY\n9CMePYmVjH50aKXT9u3btW/fPj300EPyer3y+XzttwcOHOhSIQAAe3VohHDttdcqHA6nqhYAgI3Y\nTA8AkEQgAAAsBAIAQBKBAACwEAgAAEkEAgDAQiAAACQRCAAAC4EAAJBEIAAALAQCAEASgQAAsBAI\nAABJBAIAwEIgAAAkEQgAAAuBAACQRCAAACwEAgBAEoEAALAQCAAASQQCAMBCIAAAJBEIAAALgQAA\nkEQgAAAsBAIAQBKBAACwEAgAAEkEAgDAQiAAACQRCAAAC4EAAJBEIAAALAQCAEASgQAAsBAIAABJ\nnQiEt956SxUVFRo/fnwq6gEA2MTdkQf/7Gc/09KlS3XJJZekqh4AgE06NELwer3asWOH+vfvn6p6\nAAA26dAIYcaMGamqAwBgMzYqAwAkdXCEkCin0yGn0xFzn8vlbL91u8khKbYnqRYOS/X1DvXoYeTM\n4PansyfZgH7EoyexktmPlARCaWmhHI7TAyEsSfL7vfL7C1PxtlnL7/em5X3KytLyNkmRrp5kC/oR\nj57ESkY/UhIIx483xI0QGhqaJEnBYJPCYVcq3jbruFxO+f1eqyeRhJ7z1786tHOnU2VlRldckdhz\namo8evvtvPbp225r04MPtnaq5lTrTE+6M/oRj57EOl8/SkoSXwBPSSBEIkaRiIm572Sh4XBEoRAf\n4tcl2pNjxxwaM8argwejQ8PHHmvRzJnn/2Hft0/64INT05deajL+M2A+iUU/4tGTWMnoR4d3O/X5\nfHrllVf029/+tn0a6fHWW+72MJCkpUvzzvHoU2bMaFUgEA3osrKI7r+/LSX1AchuHRohNDU1paoO\nJKCsLHLatDnLI2N961sRbd/eoM8/d+qSS8IqKUlFdQCyHZvps8hNN4U1dWqrioqMBgwIa/785oSf\nW1ZmNHw4YQDg7FKyDQGp89RTLXrqqRa7ywDQDTFCAABIIhAAABYCAQAgiUAAAFgIBACAJAIBAGAh\nEAAAkggEAICFQAAASCIQAAAWAgEAIIlAADLa7t3Rr+hXX9lcCHICgQBkqCef9OgnP4leFnHGDK/q\n620uCN0egYCsZ0z3XIJessTT/u8DB5zauJGTE7e1Sfv387OVKnQWWe3jj50aOrRQ/foV6847vWpO\n/BIRGa9nz9gLICV6QaTuqqlJuu02n6ZOjY6aNmzg2uzJRiAgq/3DP+TryJHobLxxo1srVyZ2WdFs\nsHBhk3r1il4lb8KEVl17bdjmiuz1H//h1s6dp0LgV7/ynOPR6AzGoMhqDQ2Oc05ns+98J6IRI5ok\nFeruu9sUCtldkb1crnNPo+sYISCrPfJIq9zu6KqUqqqI7ryzzeaKkCrjxoX0t397KhUfeIArByYb\nIwRktbFjQxoypEH/939OXX55WH6/3RUhVTwe6Te/aVJdnUuSTyNGRHJ+1JRsBAKyXr9+Rv365fb6\n9VzhcEgVFbm9cT2VWGUEAJDECAE5qrVV+vd/dysclm69NaTCQrsrQio1NkqHDjlUWWlUUJDa9woG\npTVr8uTzGd16a0juLPqVzaJSgeQwRpo40atNm6Kz/7/+a1hr1zbKw16M3dL+/Q5VV/t0+LBTVVUR\nvfZao6qqUrPaqbFRGjvWpz/8IboL1BtvtGnZsuw5OKZbrjI6eNChmhqPXnrJo4YGu6vJPY2N0X3G\n1693yWTg6t7Dhx3tYSBJH3zg0h/+0C2/CpD0/PP5Onw4+vnW1jr14ouJJX9Li7RwYZ6eecajzz5L\nbP6Izkun9od98808BYMdr9ku3W6EUF8v3XKLT19+Gf0A33nHpddfb7K5qtzR3Cx9//s+ffBB9Etx\n551teumlzFpCKi428vmMGhujxyzk5RmVl2dgciElEl1ImT69QGvXRg90XL7co02bGlRZee4nX3hh\nRE6nUSQSnbd69DDy+bpUblp1u8Wi//5vV3sYSNK2be6sSuhs9/vfu9rDQJJ+85u8jDvPUFGRtHRp\nk/r2jeiiiyKaP79ZvXsTCN3Vww+3qrw8esR3794RzZjRmtDz3n771PJyfb1D779//iPhBgww+sUv\nmlVZGdHAgWEtX97ENgQ79ekTUX6+UUtLNKF79YqoqMjmonJIaWnsD2thYeo34nXG6NFhvf8+6xNP\nam6Wpk4t0Lp1bg0aFNHy5U3q06d7hOTgwRG9/36DDhyIbkNIdAeCSy6J6OOPoyHgdBoNHBhJ6Hnj\nx4c0fnx2HiDR7UYIVVVGS5Y06aqrwho5MqR/+7cmObvd/zJzDR4c0dy5zSouNrrggogWLcquJaRU\naG2VZs/O1/e+59OTT3oUzsBDJpYsydN//meeQiGHPvnEpX/8x/yUvl9Li/TrX7u1YkV6RpBFRdI3\nv5l4GEjSsmVNuummNl11VVgLFjTrsssSC4Rs1i2/qmPGhDVmTGOHn/frX7vbl5AeeaRVed3nPGlp\nNW1am6ZN4xQSJ/385572E7Ht3u1SSYn04IOJrbZIl7/8xXHO6WQ6uZfXu++e3MsrT2++2aj81GZQ\nh1VVGb38cmZt/0o1lp0ta9e69dBDXq1dm6fnn8/XU09l2NyJjPDxx07NmZOvf/mXPLUlmHmn78G0\nb1/mfe3uvDOkHj2iq4hcLqO///vUBfqhQ472MJCi2/0ysSe5qFuOEDpj167YDUZfP80uIEmff+7Q\n2LG+9r2T9uxxJbQH1ZgxYf3Xf50abt5wQ+atXx44MKJNmxq0c6dLAwZENGRI6laPBAJGhYWm/cy0\n7OWVOYhly/DhsSt2r746A1f0dmORiDR3br7uvju6BfrECZsLOoPf/c7dHgaStH59YgsNEye2afHi\nJk2f3qqXX27UbbdlXiBIUu/e0SNrUxkGklRYGN3Lq3//iKqqIlqwgL28jJGee86jUaN8uv/+Atsu\nl8oIwXLzzSH98pdNWr/ercGDE981rbP273fo1Vc9+vnPpc8+c6pv3+6/wepcVqzI04IFHl1xRXT6\nl7/0aObMzFp/O2hQRA6HkTHRUBg8OPHPbNy4kMaNy8wgsMOoUWFt385eXietWpWnF16Irqb+5BOX\n8vJky/E7BMLX3H57SLffnvovbWurdPvtPl14YXSA9pOfFOhXv/pKJSUpf+uM9b//GztYPXgw8wav\nV18d1j//c7NWrsxTRYXR009zPn4kx+nXibbrutGZ963LAUePOmIOngsGHXE/iLnmpptCyss7tdpg\n5MjMXJq+666Q1q5t0qJFzbrwwtxezYHkufHGUPuFniTpllvs2UuPEYINKiqMBgwIS4qugy4ri2jA\ngNxeZTR8ePQEc3/8Y54kj26/PcTFT5Azrr46rDVrGrVhQ3SV9a232jPzpyUQjhxxaM6c6MbCxx/P\n15NPSj16pOOdM5PbLb32WpPWrMmXlKdf/KJZxcV2V2W/K6+M6DvfaZPEaUe76vBhhzweo9JSuytB\nooYNi2jYMHuPT0nLeoqf/tSjzz6LLg3v3etu33iSyyoqjB54IPrhn++EWUBHzJqVr6FDizRkSJGW\nLuXoSiQuLYFQX+84bTod7wrkng8+cOrll6MjrHDYodmz89WcWTtrIYOlJRDuuadNHk90KbigwGjS\nJE5rAKTC6UdPh8PRYzyARDiMSe4lTIwx+uKLQ3I4YkcFtbUNGjlyoLZu/UyVlVyvUJJcLqf8fq+C\nwSaFw3xrJXpyuo72wxhp3jyPtm6NriqaOLFVd9zRvRbAmEdina8fJSWFKi4ujvtNPpOkB0IwGFQg\nEEjmSwIAuqC+vl5+v/+8j0vbCKGxsUGXXjpQe/d+Jp+PEYLEks6Z0JNY9CMePYmV0SOEszk5ckg0\nqQAA6ZW2QDDG6MSJEwknFQAgvdIWCACAzJbbJ9ABALQjEAAAkggEAICFQAAASCIQAAAWAgEAIIlA\nsI3T6ZTX65XP52u//dGPfmR3WWn11ltvqaKiQuPHj4/724YNGzR8+HAFAgFddtllWrVqlQ0VptfZ\n+vHuu+/K6XTK5/PFzC+rV6+2qdL0qa2tVXV1tcrKytSrVy9NmTJFwWBQUm7OI2fqR319ffLmEQNb\nOJ1OU1tba3cZtpk3b54ZPHiw+e53v2t+8IMfxPzt0KFDpqioyCxfvty0tLSY9evXG5/PZ3bv3m1T\ntal3rn5s2rTJ9O3b16bK7DV06FBzzz33mMbGRnPw4EEzbNgwc++99+bkPGLM2fuRrHmEEYJNjDEy\nOXxMoNfr1Y4dO9S/f/+4v61cuVKDBg3SpEmT5PF4dN1112ncuHFavHixDZWmx7n6kavq6+s1bNgw\nPffcc/J6verdu7cmTZqkzZs35+Q8cq5+JAuBYKPHHntMffr0UWlpqaZOnaqGhga7S0qbGTNmqPgs\n1w3dvXu3rrzyypj7rrzySu3cuTMdpdniXP2QoucCq66uVnl5uS666CK98MILaazOHoFAQIsXL1Z5\neXn7fQcOHFBlZWVOziNn6kdtba0qKyslJWceIRBscs0112jMmDHav3+/tm/frvfee08PPPCA3WVl\nhLq6OpWUlMTcV1paqj//+c82VWQvv9+voUOH6tFHH9WhQ4e0dOlSzZ07V8uXL7e7tLTatWuX5s+f\nryeeeIJ5RNF+LFiwQLNnz07aPEIg2GTr1q2aMmWK8vLyNGjQINXU1GjVqlVqO/2SVzkql1enne6K\nK67Qhg0bdO2118rtduuGG27Q/fffr2XLltldWtps3bpVN954o2pqajR69GhJuT2PfL0fo0aNSto8\nQiBkiIsvvljhcFhHjx61uxTblZeXq66uLua+uro6XXDBBTZVlHkuvvhiffnll3aXkRZr167VLbfc\nohdffLF9FJ3L88iZ+nEmnZlHCAQbfPjhh5o1a1bMfXv37lV+fr569+5tU1WZ46qrrtLu3btj7tu5\nc6eGDx9uU0X2evXVV7Vw4cKY+/bu3at+/frZVFH6bNu2TZMnT9bq1as1YcKE9vtzdR45Wz+SNo90\neT8ldNjBgwdNcXGxqampMS0tLebTTz81Q4YMMQ8//LDdpaXd5MmT43azPHr0qAkEAmbJkiWmubnZ\nvPnmm6awsNDs2bPHpirT50z9WLNmjSksLDTr1q0zbW1t5u233zbFxcXm9ddft6nK9AiFQubSSy81\nixYtivtbLs4j5+pHsuYRAsEmW7ZsMSNGjDDFxcWmvLzc/PjHPzYtLS12l5U2BQUFxuv1Grfbbdxu\nd/v0SVu2bDGXX365KSgoMIMHD+72P37n68eiRYvMoEGDjM/nM/369TPLli2zr9g02bJli3E6ncbr\n9bb34+RtbW1tzs0j5+tHMuYRLpADAJDENgQAgIVAAABIIhAAABYCAQAgiUAAAFgIBACAJAIBAGAh\nEAAAkggEAICFQAAASCIQAAAWAgEAIEn6f7pp9Kvxj4G7AAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Pandaのみを使った計算\n", "# 中心化したデータをnumpyのマトリックスにする\n", "Xc =(X - mx).values\n", "# 標本共分散行列を求める\n", "Sx = (X - mx).cov().values\n", "# 異常値aを求めるΣ^-1 * Xcは、solveの解を利用\n", "a_prev = np.linalg.solve(Sx, Xc.T).T * Xc\n", "# 行単位の和を求め、1変数当たりの異常度を計算する\n", "a = [sum(x)/Xc.shape[1] for x in a_prev]\n", "# リストプロット\n", "list_plot(a, figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\t最も異常度の大きい、CalifのデータをXcから取り出し、SN比を計算し、棒グラフで表示します。\t\t\n", "\t

\n", "\t

\n", "\t\t表示順はdeaths fuel popden rural tempで、図2.7(b)と並びが異なりますが、\n", "\t\t値は同じ結果になっています。\n", "\t

\n", "\t

\n", "\t\tこの結果から、異常に寄与している変数は、通事故死亡者数deaths、1年度との燃料消費量\n", "\t\tfuel、人口密度popdenであることが分かります。\n", "\t\tこの結果は、車が多く使用される過密な大都市で交通事故死亡者が増加するという一般的な\n", "\t\t常識と一致します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-14.80328085 12.82597259 -6.26565335 -0.20576912 0.21784965]\n" ] } ], "source": [ "# XcからCalifのデータを取得\n", "Xc_prime = Xc[4]\n", "# 要素毎の計算をしたいので、numpyの形式に変換\n", "SN1 = 10*np.log10(Xc_prime^2/Sx.diagonal())\n", "print SN1" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEBCAYAAACXArmGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAFAJJREFUeJzt3X1wVPW9x/HPOZsHdjcczGIiYikhPCnIIEEu4o1DO1o7\nQ2+pzr0MMxBbGJoyllKpYUoFEe5U2ivSqXZsO9rS2jFY5lqL6bRQvAUdpz6QFPWKUh/ipVIvF8kT\nbLK7PGR37x9bfnGbVAibk99i3q9/Fn7n7P6++e7JfvY8bNZJp9NpAQAgybVdAAAgfxAKAACDUAAA\nGIQCAMAgFAAABqEAADAIBQCAQSgAAIy8DoV0Oq1oNCo+XwcAg6PAdgGS1NLS2ed4LNapceOu0KFD\n/6twePggV5U/XNdRJBJWe3tMqdTQDkh6kUEfetCLjHP1oazs/F5D83pPwXGcrNuhynUdOY4j1x3a\nfZDoxVn0oQe9yBioPuR1KAAABldeHD7CRzt9+rSamg4qGk0omUxZrWXq1GkqKiqyWgMA/xAKF4HX\nXz+gm276H0lTLVfyhnbvlmbMmGm5DgB+IRQuGlMlzbJdhKSY7QIA+IhzCgAAg1AAABiEAgDAIBQA\nAAahAAAwCAUAgEEoAAAMQgEAYBAKAACDUAAAGIQCAMAgFAAABqEAADAIBQCAQSgAAIwBCYXdu3dr\n1KhRWrRoUa9le/fu1ezZszVixAhNmzZNjz/++EBMCQDwQc6hcP/992vVqlWaNGlSr2VHjx7VF77w\nBX31q19VS0uLHnjgAdXW1urll1/OdVoAgA9yDoVgMKjGxkaNHz++17Jt27Zp8uTJ+tKXvqSioiLd\neOONmj9/vn7605/mOi0AwAc5h8LXvvY1DR8+vM9l+/fvV1VVVdZYVVWVmpqacp0WAOADX080t7W1\nqbS0NGssEomotbXVz2kBABeowO8J0un0OddxXUeu6/QaDwRcc1tQMHQvlOqrN7bYfi4+vE0MZfSh\nB73IGKg++BoKZWVlamtryxpra2tTeXl51lgkEpbj9BUKSUmS5wXleWH/Cs1zJSXDbJdgeF5QpaX2\nnwvPC9ouIS/Qhx70IiPXPvgaCtdee60effTRrLGmpibNnj07a6y9Pdbnu+FYLCFJikYTSiYDvtWZ\n77q6TkrKj2CIRhPq6IhZmz8QcOV5wb9tEylrddhGH3rQi4xz9eF838z5GgqLFy/Wxo0b9bOf/UyL\nFy/Wnj17tGvXLu3bty9rvVQqrVSq92Gmsz9YMplSd/fQfbL76o0t+fJc5EsdttGHHvQiI9c+5BwK\nwWBQjuPozJkzkqQdO3bIcRzF43GVlZXpt7/9rVauXKkVK1aooqJC27Zt09SpU3OdFgDgg5xDIZFI\nfOTy6upqvfLKK7lOAwAYBEP7dD0AIAuhAAAwCAUAgEEoAAAMQgEAYBAKAACDUAAAGIQCAMAgFAAA\nBqEAADAIBQCAQSgAAAxCAQBgEAoAAINQAAAYhAIAwCAUAAAGoQAAMAgFAIBBKAAADEIBAGAQCgAA\ng1AAABiEAgDAIBQAAAahAAAwCAUAgEEoAAAMQgEAYBAKAADD91BwXVfBYFChUMjc3nHHHX5PCwC4\nAAV+T+A4jt5++22NGTPG76kAADnyfU8hnU4rnU77PQ0AYAAMyjmFNWvWaOzYsYpEIlq+fLlisdhg\nTAsA6CffQ2HOnDm6+eab1dzcrBdffFEvvfSSVqxY4fe0AIAL4Ps5heeff978e/Lkybrvvvs0f/58\n/eQnP1FhYaEkyXUdua7T676BgGtuCwqG7oVSffXGFtvPxYe3iaGMPvSgFxkD1QffQ+HvVVRUKJlM\n6tixY7riiiskSZFIWI7TVygkJUmeF5TnhQe1znxSUjLMdgmG5wVVWmr/ufC8oO0S8gJ96EEvMnLt\ng6+h8Oqrr6q+vl5btmwxYwcPHlRxcbFGjx5txtrbY32+G47FEpKkaDShZDLgZ6l5ravrpKT8CIZo\nNKGODnvnhAIBV54X/Ns2kbJWh230oQe9yDhXH873zZyvoVBeXq5HHnlE5eXlWrVqlf7yl7/onnvu\n0fLly7P2DFKptFKp3lconf3BksmUuruH7pPdV29syZfnIl/qsI0+9KAXGbn2wdeDcKNHj9bOnTvV\n0NCgSy+9VNXV1Zo3b57uu+8+P6cFAFwg388pVFdXZ51sBgDkr6F9uh4AkIVQAAAYhAIAwCAUAAAG\noQAAMAgFAIBBKAAADEIBAGAQCgAAg1AAABiEAgDAIBQAAAahAAAwCAUAgEEoAAAMQgEAYBAKAACD\nUAAAGIQCAMAgFAAABqEAADAIBQCAQSgAAAxCAQBgEAoAAINQAAAYhAIAwCAUAAAGoQAAMAgFAIBB\nKAAADEIBAGAQCgAAo8B2AQCQi9OnT6up6aCi0YSSyZS1OqZOnaaioiJr8w8U66GQTqcVi3XKcZxe\ny06c6JAk7dv3goYNCw12acaUKVNUWGjvyU4k4pLekNRlrYaMN3XyZIXicXt1pFLdevXVRsVip5RK\npa3VYXubyJc+SPZ7ceDAf+vzn98l6ZPWapAOq6EhoenTr7FWQSDgKhBIKhbrOxyj0bSGDx/e52vt\nhznpdNrqFhWNRjVixAibJQDAkHDixAl5nveR61gPhXQ6rffe+78+02vfvhe0cOG/SfqrpI/+Qfyz\nXw0NJ62/A/C8oPXd43yQeVdYLGmmxSrsbxOZPhySdKW1GjLeVENDBb8feeBcfSgtDZ/XnoL1w0eO\n4ygcHt7nsp5DRp7shUKJhg1zFAqVWJpfKihw5XlhJZMBdXcP3Y1ekoLBkKSg7G0PUj5sE7NmzVZj\n48g8eCG8yvqxdH4/Ms7VB8/r+3W21+MMdGEA/FdUVKRZs2apoyM2pF8IMfC4JBUAYBAKAACDUAAA\nGIQCAMAgFAAABqEAADAIBQCAQSjkub17/6Arr6zUokWLbJdi3fvv/1UbN95tuwzrXn/9gG699V90\nySWX6KqrxusrX1miY8eO2S7Lqm984xsaOfL8Ppz1cXTZZSM0evSlCoVCGj36Uo0de5nWrfvmBT0W\noZDHHnroQa1f/y2NHz/Bdil5oaZmoUpKbH6S2b7Tp09r4cJbdcMNc9XS0qLnn29US0uL1qy503Zp\n1hw48Joee+yxc/75ho8zx3HU2Piq4vG4jhxp1XvvfaBNmzZf0GMRCnksGBym3bufUUVFpe1SrItG\nT2jGjCp9+cu1tkuxKpGIa926DVq1qk6FhYWKREbqc5/7vN5886Dt0qxIp9Oqq1uluro626VYlfkT\ndgPzZ+wIhTy2bNlylZQM3V3iD/O8Efr+9x/SJZeU2i7FqhEjLtGiRbfJdTO/uu+887a2b39ct9zy\nr5Yrs+PRR7cqGBzG4VVJGzeu19ixY1VZOUZ1dXcoFotd0OMQCsBF6P33/6ri4mJdf/0sVVXN1De/\nudZ2SYPu2LFjuv/+72rLlgdsl2Ldtdf+kz796RvV3Nys3bv3aP/+Jn3rWxe290QoABehT3xijE6d\nOqXGxlf07rvNuv32L9suadBt2LBWixd/URMnTrJdinW/+91/afHi21RYWKiJEyfpnnv+Xb/+9RM6\nc+ZMvx+LUMBF5cUXn7ddQl4ZN65Sd921Xjt2/Ert7W22yxk0zz33rJqaGnXnnZkrbCx/LUzeGTNm\nrJLJpFpbW/p9X0IBF43Gxn3avPk/bJdh1R//+Jyuvz77C4Ycx5HjOB+L7wc+X08++Z9qbW1RVdUU\nTZw4VjNnzlQ6ndaUKZVqaPi17fIG1YEDr2nDhnVZY2+99aaKi4s1atTl/X48QgEXhWQyqbq6laqt\nXW67FKumT79GnZ2d2rhxvRKJhFpbW7Rly39ozpx/HlIXJXz729/VSy+9rGeeeUHPPfeSdu7cKUl6\n5pkX9NnPzrNc3eAqKyvTY489qh/84Ps6ffq0mpvf0ebNm/TFLy69oMt0+ZKdPPbJT5bLcRxzXHDH\njh1yHEfvvfeB5coGX1NTo95552398Ic/kLTAdjnWDB/u6YknGrR27WqVlZUpHA6runquHnhgk+3S\nBpXnjZDnZb7bvaDAVThcKMdxdNlloyxXNvhGjbpcjz/+K9177wZ973ubVVxcrIULF+uuu9Zf0OMR\nCnns8OHMp1QLClyVloaH9LdsXXfdHB09elyvvfaKbrrJdjV2XXnlVfrNb3YN+W3iw8aOHavW1uiQ\n7cV1183R73//hwHZJjh8BAAwCAUAgEEoAAAMQgEAYBAKAACDUAAAGIQCAMAgFAAABqEAADAIBQCA\nQSgAAAxCAQBgEAoAAMO3UPjFL36hQCCgUCikUCikYDCoUCikP/3pT35NCQDIka9/Onvu3Lnau3ev\nn1MAAAYQh48AAIavoXD48GHdfPPNikQimjBhgrZt2+bndACAHPkWCmVlZZo8ebK2bNmiDz74QJs2\nbdLSpUv17LPP+jUlACBHvp1TmDdvnubN6/kC7YULF2rHjh36+c9/rk996lNZ67quI9ft/QXTfY3Z\nEAi4Kiiwd6QtEHCzbocytome+T98O5TRi4yB6sOAhUJ9fb1qa2slSY7jKB6P91qnoqJC+/fv7zUe\niYTlOL1/2cPh4oEqLyeeF1Rpadh2GfK8oO0SrCspGWa7BElsE/mIXmTk2ocBC4WamhrV1NSY/z/8\n8MOKRCJasGCBGfvzn/+sysrKXvdtb4/1+Q4wFjs1UOXlJBpNqKMjZm3+QMCV5wUVjSaUTA7NLyY/\nq6vrpCT7wcA2kT/oRca5+nC+b2J8O3x06tQprVy5UpWVlZo+fbqeeOIJ7dq1S42Njb3WTaXSSqXS\nfY7ng2Qype5u+xtbvtRhE9tEftaRD+hFRq598C0Uvv71r6urq0sLFizQ0aNHNW7cODU0NOiaa67x\na0oAQI58/fDa2rVrtXbtWj+nAAAMoKF9uh4AkIVQAAAYhAIAwCAUAACGryeaAX+8kQfzV1iuAfAH\noYCLytVXT1Njo+0PKlVo6tRpluYG/EUo4KJSVFSkWbNmqaMjxgeVAB9wTgEAYBAKAACDUAAAGIQC\nAMAgFAAABqEAADAIBQCAQSgAAAxCAQBgEAoAAINQAAAYhAIAwCAUAAAGoQAAMAgFAIBBKAAADEIB\nAGAQCgAAg1AAABiEAgDAIBQAAAahAAAwCAUAgJFzKHR3d2v16tUKBAJ6+umns5al02mtW7dO48eP\n18iRIzVv3jwdOnQo1ykBAD7JKRTi8biqq6vV0dHR5/KHHnpI27dv165du3T48GFNmDBBt956ay5T\nAgB8lFModHV1admyZdq6davS6XSv5Y888ojuvPNOTZo0SeFwWN/5znd08OBBNTY25jItAMAnOYVC\neXm5amtr+1x28uRJHTx4UDNmzDBjJSUlmjhxopqamnKZFgDgkwK/Hrijo0PpdFqlpaVZ45FIRK2t\nrVljruvIdZ1ej9HXmA2BgKuCAnvn5AMBN+t2KKMXGfShB73IGKg++BYKZ/V1WOnvRSJhOU7vAAiH\ni/0oqd88L6jS0rDtMuR5Qdsl5A16kUEfetCLjFz70K9QqK+vN4eLHMdRPB7/h+tGIhG5rqu2tras\n8ba2NpWXl2eNtbfH+twriMVO9ac830SjCXV0xKzNHwi48rygotGEksmUtTryAb3IoA896EXGufpw\nvm9s+xUKNTU1qqmpOa91i4uLdfXVV2v//v264YYbJEnHjx9Xc3OzZs+enbVuKpVWKtV7j6KvMRuS\nyZS6u+1vbPlSRz6gFxn0oQe9yMi1D74ehLv99tv14IMP6q233lJnZ6fWrFmjmTNnqqqqys9pAQAX\nKKdQqK+vVzAYVCgUkuM4mj9/vkKhkJYvXy5JWr58uZYsWaK5c+fq8ssv15EjR/Tkk08OSOEAgIGX\n04nm8zmctGHDBm3YsCGXaQAAg2RoX8MFAMhCKAAADEIBAGAQCgAAg1AAABiEAgDAIBQAAAahAAAw\nCAUAgEEoAAAMQgEAYBAKAACDUAAAGIQCAMAgFAAABqEAADAIBQCAQSgAAAxCAQBgEAoAAINQAAAY\nhAIAwCAUAAAGoQAAMAgFAIBBKAAADEIBAGAQCgAAg1AAABg5h0J3d7dWr16tQCCgp59+OmvZ0qVL\nVVhYqFAopFAopGAwqEgkkuuUAACf5BQK8Xhc1dXV6ujo+IfrrF+/XvF4XPF4XIlEQu3t7blMCQDw\nUU6h0NXVpWXLlmnr1q1Kp9MDVRMAwJKCXO5cXl6u2traj1xnz549euqpp9Tc3KwpU6boRz/6kaqq\nqvo5035JJRdcZ27ekFRhaW4AGFw5hcK5jB8/XgUFBbr33nsVDoe1ceNGfeYzn1Fzc7NKS0vPef8p\nU6ZIkhoaTmrYMMfPUj9ChaZOnWZpbgAYXL6Gwt133531/82bN+uXv/ylnnrqKS1dutSMu64j1+39\noj9s2DBJUlVVlcJhW3sK9gUCbtbtUEYvMuhDD3qRMVB96Fco1NfXm8NFjuMoHo/3azLXdTVmzBgd\nOXIkazwSCctxeodCIJCUJHleUJ4X7tdcH0eeF7RdQt6gFxn0oQe9yMi1D/0KhZqaGtXU1Jz3+nV1\ndVqyZImmTcscfjlz5ozeffddVVZWZq3X3h7rc08hFktIkqLRhJLJQH9K/VgJBFx5XvBvfUjZLscq\nepFBH3rQi4xz9aG09PzeWPt6+OjQoUNasWKFtm/fLs/ztH79ehUVFemWW27JWi+VSiuV6n310tkf\nLJlMqbt76D7ZZ9GHHvQigz70oBcZufbBSedwLenZw0mO4+jUqVMqLCyU67q67bbb9PDDD+v48eOq\nq6vTzp071dnZqdmzZ+vHP/6xJk2adF6PH41GNWLECJ04cUKe511omQCA85RTKPgtnU6rs7NTw4cP\n7/OcAwBgYOV1KAAABtfQvoYLAJCFUAAAGIQCAMAgFAAABqEAADAIBQCAQSgAAAxCAQBgEAoAAINQ\nAAAYhAIAwPh/Sx365IQS/gIAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 左からdeaths fuel popden rural tempで、図2.7(b)と並びが異なる\n", "bar_chart(SN1, 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 }