{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\t
\n", "\t\t井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。\n", "\t
\n", "\n", "\n", "\n", "\t\n", "\t\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\t最初に必要なライブラリーやパッケージをロードしておきます。jupyter用に新しくなったRUtil.pyも使います。\n", "\t
\n", "" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "\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\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", " | height | \n", "repht | \n", "repwt | \n", "weight | \n", "
---|---|---|---|---|
count | \n", "200.000000 | \n", "183.000000 | \n", "183.000000 | \n", "200.000000 | \n", "
mean | \n", "170.020000 | \n", "168.497268 | \n", "65.622951 | \n", "65.800000 | \n", "
std | \n", "12.007937 | \n", "9.467048 | \n", "13.776669 | \n", "15.095009 | \n", "
min | \n", "57.000000 | \n", "148.000000 | \n", "41.000000 | \n", "39.000000 | \n", "
25% | \n", "164.000000 | \n", "NaN | \n", "NaN | \n", "55.000000 | \n", "
50% | \n", "169.500000 | \n", "NaN | \n", "NaN | \n", "63.000000 | \n", "
75% | \n", "177.250000 | \n", "NaN | \n", "NaN | \n", "74.000000 | \n", "
max | \n", "197.000000 | \n", "200.000000 | \n", "124.000000 | \n", "166.000000 | \n", "
\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": [ "\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\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\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\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\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", " | _row | \n", "deaths | \n", "fuel | \n", "popden | \n", "rural | \n", "temp | \n", "
---|---|---|---|---|---|---|
0 | \n", "Alabama | \n", "1.9638 | \n", "0.5614 | \n", "0.3401 | \n", "0.3491 | \n", "0.3310 | \n", "
1 | \n", "Alaska | \n", "1.5911 | \n", "0.4470 | \n", "0.0357 | \n", "0.4294 | \n", "1.3157 | \n", "
2 | \n", "Arizona | \n", "2.0098 | \n", "0.5390 | \n", "0.1239 | \n", "0.3094 | \n", "0.5326 | \n", "
3 | \n", "Arkanas | \n", "2.0740 | \n", "0.5902 | \n", "0.3145 | \n", "0.5842 | \n", "0.4411 | \n", "
4 | \n", "Calif | \n", "1.7888 | \n", "0.1046 | \n", "0.0999 | \n", "0.1168 | \n", "0.0660 | \n", "
\n", " | state | \n", "deaths | \n", "fuel | \n", "popden | \n", "rural | \n", "temp | \n", "
---|---|---|---|---|---|---|
4 | \n", "Calif | \n", "1.7888 | \n", "0.1046 | \n", "0.0999 | \n", "0.1168 | \n", "0.066 | \n", "
\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 }