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

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

\n", "\t

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

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

5章 不要な次元を含むデータからの異常検知

\n", "\t

\n", "\t\tこの章でのポイントは、主成分分析で求めた正常な空間からのずれの距離を異常度とする\n", "\t\t再構成誤差と思われます。\n", "\t

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

分散最大化による定式化

\n", "\t

\n", "\t\t井出本のラグランジュ乗数の説明が私が習った方法と違うので、ここではPRMLの12章の定式化を元に\n", "\t\t式を整理します。\n", "\t

\n", "\t

\n", "\t\tPRMLの図12.2から1次元に投影されたD次元ベクトル$u_1$を使って、分散最大化による定式化の\n", "\t\t様子をまとめます。\n", "\t\t

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

\n", "\t\t$u_1$は、単位ベクトルであり、$u_1^T u_1 = 1$であり、各データ点は$x_n$は、$u_1^T x_n$として投影されます。($\\tilde{x}_n$の点)\n", "\t\t投影されたデータの平均値は、$u_1^T \\bar{x}$であり、$\\bar{x}$はサンプルデータの平均です。\n", "$$\n", "\t\t\\bar{x} = \\frac{1}{N} \\sum_{n=1}^N x_n\n", "$$\t\t\n", "\t\tまた、投影されたデータの分散は、以下の様に求まります。\n", "$$\n", "\t\t\\frac{1}{N} \\sum_{n=1}^N \\left \\{ u_1^T x_n - u_1^T \\bar{x} \\right \\} = u_1^T S u_1\n", "$$\t\t\n", "\t\tここで、Sはデータの共分散行列であり、以下の様に定義されます。\n", "$$\n", "\t\tS = \\frac{1}{N} \\sum_{n=1}^N (x_n - \\bar{x})(x_n - \\bar{x})^T\n", "$$\t\t\n", "\t\t投影された分散Sの最大値を求めます。ここで$u_1^T u_1 = 1$の制約が加わり、これをラグランジュ乗数を導入して、\n", "\t\t以下の様な式を$u_1$で微分し、0となる値を求めます。\n", "$$\n", "\t\tu_1^T S u_1 + (\\lambda_1 - u_1^T u_1)\n", "$$\n", "\t\t以下の公式にSが対象行列であることを使うと、第1項の偏微分は以下のようになります。\n", "$$\n", "\t\t\\frac{\\partial u_1^T S u_1}{\\partial u_1} = (S + S^T)u_1 = 2 S u_1\n", "$$\t\t\n", "\t\t右辺をまとめると、\n", "$$\n", "\\begin{eqnarray}\n", "\n", "\t\t0 & = & \\frac{\\partial}{\\partial u_1} \\left [ u_1^T S u_1 - \\lambda_1(1 - u_1^T u_1) \\right] \\\\\n", "\t\t & = & ( 2S u_1 - 2 \\lambda_1 u_1 )\n", "\\end{eqnarray}\t\t\n", "$$\n", "\t\t固有値の定義から、u_1が固有値$\\lambda_1$の固有ベクトルであることが分かります。\n", "$$\n", "\t\tS u_1 = \\lambda_1 u_1\n", "$$\t\n", "\t\t上記の式に左から$u_1^T$を掛け、$u_1^T u_1 = 1 $から、分散$\\lambda_1$(固有値)は以下の様に求まります。\n", "$$\n", "\t\tu_1^T S u_1 = \\lambda_1\n", "$$\t\t\n", "\t

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

準備

\n", "\t

\n", "\t\tいつものように必要なライブラリを読み込み、テストデータとしてRのMASSパッケージに含まれているCars93を使用します。\n", "\t

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

データ選択

\n", "\t

\n", "\t\tCars93は、93年のアメリカの車の指標と集めたもので、価格(Min.Price, Price, Max.Price)と燃費(MPG.city, MPG.highway)、\n", "\t\tエンジンサイズ、馬力など、15項目を使用します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Index: 93 entries, Acura Integra to Volvo 850\n", "Data columns (total 15 columns):\n", "Min.Price 93 non-null float64\n", "Price 93 non-null float64\n", "Max.Price 93 non-null float64\n", "MPG.city 93 non-null int64\n", "MPG.highway 93 non-null int64\n", "EngineSize 93 non-null float64\n", "Horsepower 93 non-null int64\n", "RPM 93 non-null int64\n", "Rev.per.mile 93 non-null int64\n", "Fuel.tank.capacity 93 non-null float64\n", "Length 93 non-null int64\n", "Wheelbase 93 non-null int64\n", "Width 93 non-null int64\n", "Turn.circle 93 non-null int64\n", "Weight 93 non-null int64\n", "dtypes: float64(5), int64(10)\n", "memory usage: 11.6+ KB\n" ] } ], "source": [ "# 以下の15変数を選択し、Make(製品名)をインデックスとする\n", "cc = ['Min.Price', 'Price', 'Max.Price', 'MPG.city', 'MPG.highway', 'EngineSize', 'Horsepower', 'RPM', 'Rev.per.mile', 'Fuel.tank.capacity', 'Length', 'Wheelbase', 'Width', 'Turn.circle', 'Weight', 'Make']\n", "X = cars93[cc]\n", "X = X.set_index(['Make'])\n", "X.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

データのスケーリング

\n", "\t

\n", "\t\tデータを平均と分散でスケーリングし、データ個数Nをセットします。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 中心にスケーリングしたデータ\n", "Xc = (X - X.mean())/X.std()\n", "N = Xc.shape[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

\n", "\t\tPandasの関数で共分散行列を求めます。これとSageで定義に従って計算した共分散行列Sigの値を比較し、\n", "\t\t正しく計算できていることを確認しました。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Pandasで共分散行列Sigを求める\n", "Sig = matrix(Xc.cov().values)\n", "# show(Sig)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 定義に基づきSageで散布行列Sを計算し、正規化して共分散行列Sigを表示して上記と一致することを確認\n", "Xcm = matrix(Xc.values).T\n", "S = Xcm*Xcm.T\n", "Sig = S/(N-1)\n", "# show(Sig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

固有値と固有ベクトル

\n", "\t

\n", "\t\tSageのeigenmatrix_leftを使って散布行列Sの固有値と固有ベクトルを求めます。\n", "\t

\n", "\t

\n", "\t\t固有値の値を大きい順にプロットしてみます。2〜3成分程度で十分表現できることが分かります。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEBCAYAAAB4wNK4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAHgZJREFUeJzt3X9QlPeBx/H3s4vA8mNhiZCImnhNWhw9eoTD0LGeTYLR\npqL2bKuTGEpopBY0mabEZnqTqXLxx9n2lBhtM4Z6nhEud6NzObF46VTbmtNoCWgyNU3aWPFHRH6s\nPzD8lN3n/nAl2YDK6uLDYz6vmcxmn198noV9Pvt9nt3VME3TREREPvMcVgcQEZGhQYUgIiKACkFE\nRAJUCCIiAqgQREQkQIUgIiKACkFERAJsUwimadLa2oo+NiEiMjgirA7wac3NF/qd3tZ2gb/5m5Ec\nPfohsbHxNznVjXE4DJKSYjlzpg2/316FZtfsds0Nym6VWzV7cvLAj5e2GSEYhhF0aycOh4FhGDgc\nyn6z2DU3KLtVlN0mhbB9ewT/8A8uAH77W6fFaUREbk1DvhAaGw2Ki6M5ceJS1OLiKC70f1ZJRERu\nwJAvhOZmg+7uj4dBHR0GZ8/ab0gnIjLUDflCSEvz8/d/7+u9P3Gij1Gj7HXBR0TEDoZ8IQwbBlu3\ntrNqVRcAmzZ14hjyqUVE7McWh9bYWJg7tweAyEiLw4iI3KJsUQgiIjL4VAgiIgKoEEREJECFICIi\ngApBREQCVAgiIgKoEEREJECFICIigApBREQCVAgiIgKoEEREJECFICIigApBREQCVAgiIgKoEERE\nJECFICIigApBREQCQi6EQ4cOkZOTg8fjITU1lby8PLxeLwC7d+8mOzubhIQE0tPTqaysDFp37dq1\njB07lsTERCZPnkxdXV149kJERG5YSIXg8/mYPn06EydOpLm5mcOHD9PU1ERxcTGnT59m1qxZFBcX\n09zcTFlZGYWFhb0H/aqqKkpLS9myZQuNjY3k5uaSm5tLR0fHoOyYiIiEJqRCaGhooKGhgccee4yI\niAg8Hg+zZ8/m4MGDVFRUkJaWRn5+PpGRkeTk5DBz5kzKy8sB2LBhAwUFBWRlZREVFcXixYsxDIOq\nqqpB2TEREQlNSIUwcuRI7r33XjZs2EBbWxtNTU1s27aN3NxcamtryczMDFo+MzOTmpoagD7zDcMg\nIyOjd76IiFgrIpSFDcNg69atTJkyhbKyMgDuv/9+VqxYwaxZsxg9enTQ8klJSbS0tADg9XrxeDxX\nnH+Zw2HgcBh9frbT6ei9jYiw17XwT2a3G7tmt2tuUHarKHuIhdDd3c2MGTOYO3cu//RP/8RHH31E\ncXEx8+bNA8A0zauuf635AElJsRhGf4XgA8DtduF2x4YSe8hwu11WR7huds1u19yg7Fb5LGcPqRB2\n7dpFfX09K1asACAuLo6lS5eSkZHBww8/3Ptuo8u8Xi8pKSkAJCcn9zs/PT09aNqZM239jhDa2i5d\nfG5t7cDnc4YS23JOpwO32xXI7rc6Tkjsmt2uuUHZrXKrZvd4Bv4COqRC8Pl8+P1+/H4/DseloUln\nZyeGYTBlyhQ2bdoUtHxNTQ3Z2dkAZGVlUVtbS15eHgB+v5+6ujrmz58ftI7fb+L39x1JXN5Jn89P\nT4+9flmXKfvNZ9fcoOxW+SxnD+mE08SJE4mLi2PJkiV0dHTg9XpZsWIFX/nKV8jLy+PYsWNs3LiR\nrq4uqqur2blzJwsWLACgqKiIzZs3c+DAATo6Oli2bBnR0dFMnz79usOLiEj4hFQISUlJvP766+zd\nu5dRo0aRnp5OTEwMlZWVDB8+nB07dvDiiy+SmJhISUkJFRUVjB8/HoBp06axcuVK5syZw2233cau\nXbuorq4mKipqUHZMRERCY5gDudJ7EzU3X+h3env7R4wZk0p9/SliYuJucqobExHhwOOJ5ezZNtsN\nRe2a3a65QdmtcqtmT06OH/B27Pf+KhERGRQqBBERAVQIIiISoEIQERFAhSAiIgEqBBERAVQIIiIS\noEIQERFAhSAiIgEqBBERAVQIIiISoEIQERFAhSAiIgEqBBERAVQIIiISoEIQERFAhSAiIgEqBBER\nAVQIIiISoEIQERFAhSAiIgEqBBERAVQIIiISoEIQERFAhSAiIgEqBBERAVQIIiISoEIQERFAhSAi\nIgEqBBERAVQIIiISoEIQERFAhSAiIgEqBBERAVQIIiIScF2FsHz5clJTU4mPj2fq1KkcO3YMgN27\nd5OdnU1CQgLp6elUVlYGrbd27VrGjh1LYmIikydPpq6u7sb3QEREwiLkQli/fj2VlZXs2bOHhoYG\nxo0bx5o1azh9+jSzZs2iuLiY5uZmysrKKCws7D3oV1VVUVpaypYtW2hsbCQ3N5fc3Fw6OjrCvlMi\nIhK6kAth9erVrFixgnvuuYe4uDjKysooKyujoqKCtLQ08vPziYyMJCcnh5kzZ1JeXg7Ahg0bKCgo\nICsri6ioKBYvXoxhGFRVVYV9p0REJHQhFcKpU6c4evQoXq+X8ePHM3z4cObMmUNLSwu1tbVkZmYG\nLZ+ZmUlNTQ1An/mGYZCRkdE7X0RErBURysInT54EYOvWrezevRufz8c3vvENCgsLaW9vZ/To0UHL\nJyUl0dLSAoDX68Xj8Vxx/mUOh4HDYfT52U6no/c2IsJe18I/md1u7JrdrrlB2a2i7CEWgmmaADz7\n7LPcfvvtAJSWlvLwww/z0EMP9c6/1vpXk5QUi2H0Vwg+ANxuF253bCixhwy322V1hOtm1+x2zQ3K\nbpXPcvaQCuGOO+4AICEhoXfamDFjME2Tixcv4vV6g5b3er2kpKQAkJyc3O/89PT0oGlnzrT1O0Jo\na7t08bm1tQOfzxlKbMs5nQ7cblcgu9/qOCGxa3a75gZlt8qtmt3jGfgL6JAKYdSoUbjdbg4dOkRG\nRgYAR48eJTIykq997Wts3rw5aPmamhqys7MByMrKora2lry8PAD8fj91dXXMnz8/aB2/38Tv7zuS\nuLyTPp+fnh57/bIuU/abz665Qdmt8lnOHtIJJ6fTyRNPPMHy5cs5cuQITU1NPP/88+Tl5fHtb3+b\nY8eOsXHjRrq6uqiurmbnzp0sWLAAgKKiIjZv3syBAwfo6Ohg2bJlREdHM3369OsOLyIi4RPSCAFg\n5cqVdHd3c99999HT08M3v/lNXnjhBWJiYtixYwdPPvkkCxcuZMyYMVRUVDB+/HgApk2bxsqVK5kz\nZw7Nzc1MmDCB6upqoqKiwr5TIiISOsMcyJXem6i5+UK/09vbP2LMmFTq608RExN3k1PdmIgIBx5P\nLGfPttluKGrX7HbNDcpulVs1e3Jy/IC3Y7/3V4mIyKBQIYiICKBCEBGRABWCiIgAKgQREQlQIYiI\nCKBCEBGRABWCiIgAKgQREQlQIYiICKBCEBGRABWCiIgAKgQREQlQIYiICKBCEBGRABWCiIgAKgQR\nEQlQIYiICKBCEBGRABWCiIgAKgQREQlQIYiICKBCEBGRABWCiIgAKgQREQlQIYiICKBCEBGRABWC\niIgAKgQREQlQIYiICKBCEBGRABWCiIgAKgQREQlQIYiICKBCEBGRgOsuhKeffhqH4+PVd+/eTXZ2\nNgkJCaSnp1NZWRm0/Nq1axk7diyJiYlMnjyZurq6608tIiJhd12FcOjQIV555RUMwwCgoaGBWbNm\nUVxcTHNzM2VlZRQWFvYe9KuqqigtLWXLli00NjaSm5tLbm4uHR0d4dsTERG5ISEXgmmaFBUVUVJS\n0jutoqKCtLQ08vPziYyMJCcnh5kzZ1JeXg7Ahg0bKCgoICsri6ioKBYvXoxhGFRVVYVvT0RE5IaE\nXAgvvfQSLpeLRx99tHdaXV0dmZmZQctlZmZSU1MDQG1tbdB8wzDIyMjonS8iItaLCGXhxsZGli5d\nyp49e4Kme71eRo8eHTQtKSmJlpaW3vkej+eK8z/J4TBwOIw+051OR+9tRIS9roV/Mrvd2DW7XXOD\nsltF2UMshJKSEp544gnS0tI4duxY0DzTNK+67rXmX5aUFNt7beKTnE4fAG63C7c7doCJhxa322V1\nhOtm1+x2zQ3KbpXPcvYBF8KuXbvYt28fL7/8MhB8gE9OTsbr9QYt7/V6SUlJuer89PT0Pj/nzJm2\nfkcIbW2XLkC3tnbg8zkHGntIcDoduN2uQHa/1XFCYtfsds0Nym6VWzW7xzPwF9ADLoSKigqampq4\n8847AfD7/ZimSUpKCiUlJX3eZlpTU0N2djYAWVlZ1NbWkpeX17tuXV0d8+fP7/Nz/H4Tv7/vaOLy\nTvp8fnp67PXLukzZbz675gZlt8pnOfuATzitWbOGP//5z7z99tu8/fbbVFdXA/D222/z6KOPcuzY\nMTZu3EhXVxfV1dXs3LmTBQsWAFBUVMTmzZs5cOAAHR0dLFu2jOjoaKZPn37dwUVEJLwGPEJISEgg\nISGh9/7FixcxDIMRI0YAsGPHDp588kkWLlzImDFjqKioYPz48QBMmzaNlStXMmfOHJqbm5kwYQLV\n1dVERUWFeXdEROR6hXRR+ZPuuusufD5f7/1JkyZx8ODBKy6/YMGC3hGDiIgMPfZ7f5WIiAwKFYKI\niAAqBBERCVAhiIgIoEIQEZEAFYKIiAAqBBERCVAhiIgIoEIQEZEAFYKIiAAqBBERCVAhiIgIoEIQ\nEZEAFYKIiAAqBBERCVAhiIgIoEIQEZEAFYKIiAAqBBERCVAhiIgIoEIQEZEAFYKIiAAqBBERCVAh\niIgIoEIQEZEAFYKIiAAqBBERCVAhiIgIoEIQEZEAFYKIiAAqBBERCVAhiIgIoEIQEZEAFYKIiAAq\nBBERCQi5EI4fP87s2bMZPnw4I0aMoKCggNbWVgB2795NdnY2CQkJpKenU1lZGbTu2rVrGTt2LImJ\niUyePJm6urrw7IWIiNywkAthxowZJCUlceLECWprazl8+DDPPPMMp0+fZtasWRQXF9Pc3ExZWRmF\nhYW9B/2qqipKS0vZsmULjY2N5ObmkpubS0dHR9h3SkREQhdSIZw/f54JEyawcuVKXC4Xqamp5Ofn\ns2fPHioqKkhLSyM/P5/IyEhycnKYOXMm5eXlAGzYsIGCggKysrKIiopi8eLFGIZBVVXVoOyYiIiE\nJqRCSEhIoLy8nOTk5N5pJ06cYOTIkdTW1pKZmRm0fGZmJjU1NQB95huGQUZGRu98q7S3w49+FMU/\n/qOLn/98mKVZRESsFHEjK7/11lusW7eO7du3s2rVKkaPHh00PykpiZaWFgC8Xi8ej+eK8y9zOAwc\nDqPPz3I6Hb23ERHhuxb+z/8cycaNl4pg794Ibr/dYO7cnrBtH4Kz241ds9s1Nyi7VZT9Bgph7969\nzJw5k1WrVvHggw+yatUqTNO86jrXmg+QlBSLYfRXCD4A3G4Xbnfs9YXux5/+FHz/gw+i8Hiiwrb9\nT3K7XYOy3ZvBrtntmhuU3Sqf5ezXVQhVVVXk5eWxfv165s2bB0BycjJerzdoOa/XS0pKylXnp6en\nB007c6at3xFCW9uli8+trR34fM7rid2vyZOH8eabkQAYhsmXvtTF2bO+sG0fLrW22+0KZPeHdduD\nza7Z7ZoblN0qt2p2j2fgL6BDLoR9+/bx+OOPs23bNnJycnqnZ2VlsWnTpqBla2pqyM7O7p1fW1tL\nXl4eAH6/n7q6OubPnx+0jt9v4vf3HUlc3kmfz09PT/h+WT/4QRfDh/t5/30HOTk9fOUrPnrCe8ao\nV7iz30x2zW7X3KDsVvksZw/phJPP56OwsJBVq1YFlQHAvHnzqK+vZ+PGjXR1dVFdXc3OnTtZsGAB\nAEVFRWzevJkDBw7Q0dHBsmXLiI6OZvr06dcdPhwMA/LzL7JiRRc5OeEdGYiI2ElIhfDmm2/y3nvv\n8dRTT+FyuYiJiem97ezsZMeOHbz44oskJiZSUlJCRUUF48ePB2DatGmsXLmSOXPmcNttt7Fr1y6q\nq6uJihqc8/UiIhKakE4ZTZo0CZ/vyq+iR48ezcGDB684f8GCBb0jBhERGVrs9/4qEREZFCoEEREB\nVAgiIhKgQhAREUCFICIiASoEEREBVAgiIhKgQhAREUCFICIiASoEEREBVAgiIhKgQhAREUCFICIi\nASoEEREBVAgiIhKgQhAREUCFICIiASoEEREBVAgiIhKgQhAREUCFICIiASoEEREBIMLqALe69993\nsHXrMEaPhnnzwOm0OpGISP9UCIPo5EmD6dNjaG01ANizJ4ry8g6LU4mI9E+njAbR/v3O3jIA+PWv\nNTwQkaFLhTCIvvAFPw6H2Xt/7Fi/hWlERK5OhTCIvvhFPz//eSfZ2T6+/nX493/vCuv2Dx1ycN99\nsYwZE8ePfxwV1m2LyGePCmGQzZ7dw86dnfz3f8Po0ea1VwhBcXE09fUO2tsNXnopkt/8JnynpEwT\n/uVfIpk0ycXcuXD2bNg2LSJDlC4q25jXG9znLS3GFZYM3X/+ZwSrV18adbz7Lvh8UfziF7ogLnIr\n0wjBxp54orv3/++8089DD/nCtu2jR4P/NI4cCV/ZiMjQpBGCjf3wh91MnOijqcng/vt7SEoK37Yf\neqiHdesiuXjxUhHMmBG+shGRoUmFYHOTJg3OgTory09VVTu//e0w7r03kqlTL9LTMyg/SkSGCBWC\nXFFmpp/77ruIxxMZ9ovKf/iDgx//OJrubnj22S6mTQt/sX30Udg3KXJL0zUEuek6OyEvL4a6Oid/\n/KOT+fNdnD4dvmsUp04ZTJzoIj4e7r8/OqwX20VuZSoEuenOnTM4e/bjg3RXl8GpU+E7aK9aFcV7\n7136037nHSdr1kSGbdsA//u/Th5+OIZvftPFn/6kp5DcOvTXLDfd7bebfOlLH1+Q+PznfWH9FPen\nTxV99FH4yqa+3mD+fBe1tU727IngkUdcmOH9eMmg6uyEujoHH36oUZP0pUKQm84w4NVXO1i+vJMl\nSzqpqmonJiZ82//e97qJi7t0lE5IMJk/v/saawzcsWMOurs/PpieOuWgrS1smx9UFy7Aww/H8NWv\nxnLffbH8z/+E9xLi737nZOxYFy4XrFo1LKzbHmx+P+zZ4+D3v2dQCt40Lz3+Q50KQSwREwOFhRdZ\nuPBiWN8uCzBhgp833+zg9ddh//4O0tPDN/r4u7/zkZr68fYmTeohLi5smx9U27YN4/DhS59mv3jR\nYMWK8H7dyYIFLpqaHHR2wqpVkdTWhvfw8u67DjZtGsZbb4V3u6YJBQXRfP3rLu6/H4qKwvu4/OUv\nDrKyYrn77nhmzHCF/c0O//EfEdx/fzTTpsFf/3pjIz8VgtySRo40mTr10umpcEpMhF/9qp1nnuli\nyZJOtmwJ76e3T540mDkzmrvugtLS8L7Kjow0r3r/RvT0wPnzwdPOnQvfaan9+51MmxbDD38YTW5u\nDNu3h2908/77Dnbu/Pix/q//iuDkyfBlX7o0ihMnLh1qDxyIYMOG8F3TOnTIwfe/H8077zj59a8h\nP//GymxIve3UNE3a2i5gGH1/Ge3tbUG3duJ0OnA6fbS1deDz2esbT+2afTBzezywaNHH99vbw7ft\nJ5+MZu/eS3lfeAHuuaczbB8K/NrXYOvWTv7v/5zExZksWdJJe3v4Hpv8/G7+7d8uHVjHjfOTkdER\ntsfm1Vcj6eq6dN3J74ctW3qYMiU8XxYZEWEAPcCl447DYWIY7WHLfuFCD/Dx94ydO9dNe/vFsGz7\nvfecmObHfx8ffGDS/qngra0m8fHx/R5XP80wzaFzSay1tZWEhASrY4iI3FLOnz+P2+2+5nJDqhBM\n0+TYsYYrjhDGjfs87777F2JiYi1Id/2cTgdut4vWVnu9ygb7Zrdr7hdfHMa//uulUwrx8SavvdbB\n3XcPmafoNQ3W497VBSUlUezb5+Rv/9bHCy904fGEbfMAmKaDhATXoIwqT540qK83GD/eH/bc9fUG\nVVXDGDFiGN/4RgeGEZzd44m15wgBoLm5/0vx7e0fMWZMKvX1p4iJsclVvICICAceTyxnz7bR02Of\ngxPYN7tdcwP85jfDaGqK5stfbueuu+z1HVJ2ftxv1ezJyfED3064g4nIjfnqV314PHD2rKnvj5Kb\nasiNEK7k8vWFgZ4LExGR0NimEEzT5MKFCwM+FyYiIqGxTSGIiMjg0gfTREQEUCGIiEiACkFERAAV\ngoiIBKgQREQEUCEMuuPHjzN79myGDx/OiBEjKCgooLW11epYIXn66adxOOz1p7J8+XJSU1OJj49n\n6tSpHDt2zOpIA3Lo0CFycnLweDykpqaSl5dHS0uL1bH69frrr3PHHXfw6KOP9pm3e/dusrOzSUhI\nID09ncrKSgsSXtnVsv/+979n4sSJJCQkcPfdd7N8+XILEl7Z1bJfZpomWVlZPPjggyFt217Pchua\nMWMGSUlJnDhxgtraWg4fPswzzzxjdawBO3ToEK+88oqtPvuxfv16Kisr2bNnDw0NDYwbN441a9ZY\nHeuafD4f06dPZ+LEiTQ3N3P48GGamppYuHCh1dH6+OlPf8r3v/99vvCFL/SZd/r0aWbNmkVxcTHN\nzc2UlZVRWFhIXV2dBUn7ulr2EydOkJubS0FBAWfOnOHVV1/lZz/72ZAptKtl/6R169Zx5MiRkLev\nQhhE58+fZ8KECaxcuRKXy0Vqair5+fns2bPH6mgDYpomRUVFlJSUWB0lJKtXr2bFihXcc889xMXF\nUVZWRllZmdWxrqmhoYGGhgYee+wxIiIi8Hg8zJ49m4MHD1odrQ+Xy8Uf/vAH7r777j7zKioqSEtL\nIz8/n8jISHJycpg5cybl5eUWJO3ratkbGxspLCyksLAQp9PJhAkTmDJlypB5zl4t+2UNDQ0sX76c\np556KuTtqxAGUUJCAuXl5SQnJ/dOO378OCNHjrQw1cC99NJLuFyuqw5Nh5pTp05x9OhRvF4v48eP\nZ/jw4XzrW98asqddPmnkyJHce++9bNiwgba2Npqamti2bRszZsywOlofixYtIj6+/y9Nq62tJTMz\nM2haZmYmNTU1NyPaNV0te1ZWFqtXrw6aduLEiSHznL1a9suefvppioqK+NznPhfy9lUIN9Fbb73F\nunXreO6556yOck2NjY0sXbqUX/ziF1ZHCcnJkycB2Lp1K7t37+add97h5MmTfPe737U42bUZhsHW\nrVt57bXXcLvdjBgxAp/Px4oVK6yOFhKv14vnU9/xnJSUZItS/rQXX3yRv/71r3zve9+zOsqAvP76\n69TV1fGjH/3outZXIdwke/fuZdq0afzkJz/hgQcesDrONZWUlPDEE0+QlpZmdZSQXP4mlmeffZbb\nb7+d1NRUSktL2b59O93d3Ranu7ru7m5mzJjB3LlzOX/+PB9++CFut9tWI7TLboVvxFm3bh1Llixh\n+/btQaP8oaqrq4tFixaxbt06IiOv75/p1Ndf3wRVVVXk5eWxfv165s2bZ3Wca9q1axf79u3j5Zdf\nBuz15L7jjjsAgv7lvTFjxmCaJk1NTYwaNcqqaNe0a9cu6uvre0cEcXFxlJaWkpGRwblz50hMTLQ4\n4cAkJyfj9XqDpnm9XlJSUixKFLrnnnuOTZs28bvf/Y4vfvGLVscZkGXLlpGZmcnUqVOB63veqhAG\n2b59+3j88cfZtm0bOTk5VscZkIqKCpqamrjzzjsB8Pv9mKZJSkoK69atY86cORYnvLJRo0bhdrs5\ndOgQGRkZABw9epRhw4aRmppqcbqr8/l8+P1+/H5/79t8Ozs7bfUOL7h0Hn7Tpk1B02pqasjOzrYm\nUIhWr17Nq6++yv79+4f0C4hPq6io4OzZs72jma6uLjo7O0lJSeHgwYMDuw5iyqDp6ekxx40bZ778\n8stWRwnJuXPnzA8//LD3v/3795uGYZinTp0yOzo6rI53TT/4wQ/Me+65x/zggw/MxsZG88tf/rI5\nf/58q2Ndk9frNZOTk83nnnvObG9vN1taWsxZs2aZDzzwgNXRrujxxx83H3nkkaBpTU1NZkJCgvnL\nX/7S7OzsNH/1q1+ZsbGx5h//+EeLUvavv+xHjhwx4+PjzcOHD1uUamD6y97Y2Bj0vF2zZo05ceJE\n89SpU6bf7x/QdlUIg+iNN94wHQ6H6XK5zOjo6KDb48ePWx1vwOrr602Hw2F1jAHr6uoyFy1aZCYl\nJZlut9v8zne+Y7a1tVkda0Dq6urMBx54wExKSjJHjBhhPvLII2ZDQ4PVsfq4/HccERFhRkRE9N6/\n7I033jAzMjLM6Ohoc+zYseZrr71mYdpgV8v+/PPPm06n03S5XL3/Xd6HoeBaj/snbdq0KeQXE/r3\nEEREBNC7jEREJECFICIigApBREQCVAgiIgKoEEREJECFICIigApBREQCVAgiIgKoEEREJECFICIi\ngApBREQC/h+4JMkes8VwDwAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 固有値Lamと固有ベクトルUを求める\n", "(Lam, U) = S.eigenmatrix_left()\n", "# 固有ベクトルのプロット\n", "list_plot(Lam.diagonal(), figsize=4, zorder=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

第1と第2固有ベクトル

\n", "\t

\n", "\t\t各点を第1と第2固有ベクトルに投影します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "93 x 2 dense matrix over Real Double Field (use the '.str()' method to see the entries)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 固有ベクトルの第1成分と第2成分上にプロットする\n", "U12 = U.submatrix(0, 0, 2, 15)\n", "Z = Xcm.T*(U12.T); Z" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

第1と第2固有ベクトル空間での分布

\n", "\t

\n", "\t\t観測値を第1と第2固有ベクトル空間での分布をプロットしてみます。井出本の分布と比べると第1成分が正負逆になって\n", "\t\tいますが、分布は同じように思われます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xd4FNX6wPHvzOxukk0jQAihhCIoYAUFQVEQREGUqqJI\nUVEsKKhclSJFQLlUr1csqKBehB+IiqgUr2ABQYqIXjSIIKIUEQgxZdN2Z+f3x8iGUJJsyZbk/TwP\nT9hkZ/Y9O+Wdc+bMOYphGAZCCCEEoIY6ACGEEOFDkoIQQggPSQpCCCE8JCkIIYTwkKQghBDCQ5KC\nEEIID0kKQgghPCQpCBFEhmGQnZ2NPB4kwpXF1wWPHs1BVRWqV4/l+HEHbndk7+RSlvBVmcrjcOTQ\nqFFdfv31ILGx8aEOxy+VabtUlbIkJ5e9z/lVU1BVBUVRUFXFn9WEBSlL+KpM5VEUpcTPSFaZtouU\n5aTlAxyPEEJUOg4HfPGFxs6dlf+UWflLKIQQfsjJgRtusHPrrXY6drTz+uvWUIdUoSQpCCFEKVas\nsLBzpwaAYSjMnm0LcUQVS5KCEEKUIja25Ou4uNDEESySFIQQohTdu7u4+WYnANWqGcyeXRDiiCqW\nz11ShRCiKlBVeOmlAmbOLCA62nxdmUlSEEKIcrDbQx1BcFTynCeEEMIbkhSEEEJ4SFIQQgjhIUlB\nCCGEhyQFIYQQHpIUhBBCeEhSEEII4SFJQQghhIckBSGEEB7yRLMQXkpJSSQqKgpFUTAMA0VRGDBg\nMM88Mz3UoQnhN0kKQnhJURS+/vpb6tatF+pQhAg4aT4SwkuGYWAYkT2PrxBnI0lBCB9MnjyeVq3O\n59xz0xg5cgQOhyPUIQkRED43H6mqgqaZOeXEz0gmZQlf4Vae1q0vp1Ona3nlldfZt+9X7r57EGPG\n/IMXX5xb5rInl8ViCY/y+Crctos/pCzFFMPHevCJG2xCVHWrV6+mR48eOBwOrNbS5+/Nzs4mMTGR\nrKwsEhISghShEOXnc03h+HEHVqtGQkIM2dn56Lo7kHEFnaapUpYwFe7lSUqqha7r/PzzPurUqVPq\nex2OfIC/y6IFI7wKE+7bxRtVpSxJSbFnWaqYz0nB7TY8H6jrblyuyP4iT5CyhK9wKM+OHf/j3XeX\n8PTTz3h+l56+k6ioKJKTU8qMT46Z8CZlkRvNQnglOTmZBQve5IUX/kVRURG//LKb6dOfYdCgu6Q5\nVVQKkhSE8ELt2qksWvQuq1evoFmzRtx00/V07nwd48ZNCnVoQgSEPLwmhJfatm3HihWfhjoMISqE\n1BSEEEJ4SFIQQgjhIUlBCCGEhyQFIYQQHpIUhBBCeEhSEEII4SFJQQghhIckBSGEEB6SFIQQQnhI\nUhBCCOEhSUEIIYSHJAUhhBAekhSEEEJ4SFIQQgjhIUlBCCGEhyQFIYQQHpIUhPDDuHGjSElJDHUY\nQgSMJAUhfLRjx/9YunSxzM0sKhVJCkL4wDAMnnjiUR544OFQhyJEQElSEMIHb745j5iYGPr0uSXU\noQgRUBZfF1RVBU0zc8qJn5FMyhK+wq08R478ycyZU/n440+wWMyYTvwsy8llKe8y4Srctos/pCzF\nFMMwDF8WNAxD2lJFlTRgwADS0tJ49tln+e2332jcuDG6rpdr2ezsbBITE8nKyiIhIaGCIxXCez7X\nFI4fd2C1aiQkxJCdnY+uuwMZV9BpmiplCVPhVJ4vv/ycr77awMaN/yIz08FffzkAyMx0lGt5hyMf\n4O+yaBUWZzCE03bxV1UpS1JSbJnL+5wU3G7D84G67sbliuwv8gQpS/gKh/K8884Sjh49ykUXNQPA\n7XZjGAbnntuQqVNn0rNnn1KXD9Qx43DA229bKSxU6N/fSc2aPlX4AyIctkugSFn8SApCVEWTJ09l\n9OhxntcHDx7ghhuu5fPPN5KYWC0oMRgG3H57DJs2mYfvokVW1q51EFv2RaAQZZKkIIQXEhISSUgo\nfljN6XSiKAopKbWDFsOxY4onIQDs3auSnq7SunXluMIVoRX5t9qFCKH69dM4fPivoH5mYqJBjRrF\nCSAqyqBu3dA1H4nKRZKCEBHGZoOFC/Np3Vrnoot05s3Lp04dSQpV3bZtKi+/bGXTJv9O69J8JEQE\natXKzYoVeaEOQ4SJzz7TuOOOGHRdQVEMli6FTp18W5fUFIQQIsK9954VXTefGzMMhUWLfF+XJAUh\nhIhw9euX7GTQoIHv65LmIyGEiHAjRhSxf7/Kpk0al17qZtIkC06nb+uSpCCEEBEuJgZefLEAMMfh\niouzkJnp27qk+UgIIYSHJAUhhBAekhSEEEJ4SFIQQgjhIUlBCCGEhyQFIYQQHpIUhBBCeEhSEEII\n4SFJQQghhIckBSGEEB6SFITw0g8/7KBv3x40aVKfCy5oytChd3LkyJFQhyVEQEhSEMILRUVF9OvX\nm6uuupqdO/eybt0mjh49ypNPPhbq0IQICEkKQnghPz+PsWMnMHz4Y1itVqpXr0H37jfx00/poQ5N\niICQUVKF8EJiYjX69x/oeb1nz24WL15Er159QxiVEIHjc1JQVQVNMysaJ35GMilL+ArH8hw4sJ/L\nLrsYXdcZPPguxox5CkVRylzu5LJYLOFTHl+E43bxlZSlmGIYhk8zfhuGUa6DQIjK7JdffmHo0KHU\nrl2bhQsXlvn+7OxsEhMTycrKIiEhIQgRCuEdn5NCRkYuVqtGQkIM2dn56Lq77IXCmKapUpYwFe7l\n2bp1C127dmb37n1Ur16j1Pc6HLnUr1+b/fsPExsbF6QIK0a4bxdvVJWyJCXFlrm8z81Hbrfh+UBd\nd+NyRfYXeYKUJXyFQ3m++modTzzxKBs3bvP8zu02a82qaikzPjlmwpuURXofCeGViy++hJycHCZN\nGk9+fj7Hjh1j5sx/0q7dlcTFxYc6PCH8JklBCC/ExyewdOlytm/fRvPmjenQoS2JidV45ZV5oQ6t\nUsjNDXUEQrqkCuGlZs2as2zZilCHUan8+adCv34xpKdrNG+u8847+aSk+HS7U/hJagpCiJCbMcNG\neroGwM6dGjNm2EIcUdUlSUEIEXIOh1LqaxE8khSEECF3331FxMebzUXx8Qb33VcU4oiqLrmnIIQI\nuJwceP11G4WFMHiwk9TU0u8PXHKJm40bHaSnqzRv7qZ2bbmfECqSFIQQAWUYcOutdrZtM+8RLF1q\n5YsvHMSX0WM3JcUgJUUPQoTh75tvVHbvVrniCp0GDYKbICUpCCEC6uhRxZMQAPbvV0lP17j8cjnh\nl8fixRZGjIjGMBTi4w1WrMijWbPgPVAn9xSEEAGVlGRQq1bxSSwmxqBBg8rxlHAwvPWWDcMwb7Tn\n5Ci8/35wr90lKQghAspqhcWL8+nQwUXbti7eeitf7hF4ITnZfcpraT4SQkS4Cy5ws3RpfqjDiEhT\npxZy9Kh5T6FLFxd33ukM6udLUhAiQL77TmXGjChUFZ58spALLpAmE+G9unUNVq3KC9nnS1IQIgCy\nsqBfPzuZmWZb8LZtKt9848BuD3FgQnhJ7ikIEQAHDqiehABw7JjK4cPyVK6IPJIUhAiARo3cpKUV\nNxc1aaJTr57cXBWRR5qPhAgAux0+/DCPuXNtqCrcf38RNhnTTUQgSQpCBEidOgZPP10Y6jCE8Is0\nHwkhhPCQpCCElw4c2M+dd95Bs2YNueCCpgwf/gA5OdmhDkuIgJCkIISXBgzoR1JSEtu372TNmnXs\n2rWTiROfCnVYQgSEJAUhvJCdnUXLlq0YO3YiMTEx1K6dyq239ufrrzeEOjQhAkJuNAvhhYSERJ57\nbk6J3x08eIDU1DpnXSYjQ0HXoVYt6aIqwp/PSUFVFTTNrGic+BnJpCzhK5zLs337t8yf/yr/939L\nsVhOj2/OHAsTJpijXg4b5mTUKHP4Ak1Tz/j+SHK27fLrrwoLF1pISIB773USExOK6LwTzvuYt/wt\ni2IYhk+XL4ZhoCjyxKaoujZs2MCNN95Bu3ZLad68NfffD02bFv89MxNq1DAnnTlh8+ZsLr88kays\nLBISEoIfdAU7cgQuvND8CdC1K6xaFdqYKlpuLgwbBtu3Q5cuMH06aFrZy4Urn2sKx487sFo1EhJi\nyM7OR9cje/AvTVOlLGEqHMuzevVK7r//XpKTd7FqVQqrVsGCBQYbN+ZRo4b5nowMMIzYEsv99VcB\nwN9lieAzB2feLmvXahw5Eu15z+rVcOSIA6s1VFGWjz/72BNP2PjPf8wC7tgB1asX8uCDrooIs1xK\nK0tSUuxZlirmc1Jwuw3PB+q6G5crPA5Wf0lZwle4lGfLls0MG3YfL764iIEDUzy/P3pUYft2hY4d\nzRnGEhNhxIhCnn8+CoB+/Zycd575t3ApSyCcXJYGDcBiMXC5zFaExo3dKIobV+jOkV7xZbvs3l2y\nxWTPHiUstq2v+1jkN6AJEUS6rjNy5MOMGzeJ66+/qsSMYtHRBuecU/IgHDu2iK+/zuXLLx288EJB\nQGM5dEihTZtY0tLi6No1htzcgK7eJ02bunn11QJat9bp3NnFwoWhGwI6WLp3L854qmrQrVuEZMCz\n8PmewtGjOVgsKklJsWRmOsIiM/pDyhK+wqk8mzZ9Ta9e3YiKisIwDAzjHJzOZ4BYXnyxKX37Vi91\n+by8XBo2rMO+fYew2+P8iqVNm1j27Su+ruvZ08lrrwU28ZwqI0Nh/PgoDh5UuO02nYcfjgrodklP\nV5k1yxw0auTIIlq0CM729ncfW7XKwo4dKlddpdOuXWjnoi6tLMnJ8WUvX1GBCVEZtW3bjsOH/wp1\nGACnDc29e3fFV/wffjiaNWvM08bGjRYuuAAuuSQw687NhVtuieHoUbMcX3+tsWWLgzj/cmdQdOvm\nolu3UEcRGNJ8JESEatPm5CtSg379Kn7axh9/LHnK+OGHwK370CHVkxDAnJPi0CE5RQWbfONCRKhF\ni/K5664i2rRxMX16IQ88UPFJoXPn4vbyqCiDjh0Dt+60tJJzUqSlualfP/KbPyONNB8JEaFsNpg2\nLbhDdU+fXsi557o5eFClb1+dCy+MITMzMOuOjoYPPsjjpZdsGAYMG1YUEQ++VTaSFIQQ5WaxwP33\nO//+f+AbGurVM3j2WZmTIpSk+UgIIYSHJAUhhBAekhSEEEJ4SFIQQgjhIUlBCCGEhyQFIYQQHpIU\nhBBCeEhSEEII4SFJQQghhIckBSGEEB6SFIQQQnhIUhBCCOEhSUGIv+3dq/Df/2ocOaKU+r7PPlvD\n+ec34f777w5SZJEpMxNef93KwoVWiopCHY0oLxklVQjg00817rorhqIihRo13Hz8cR7nnHP6TLVz\n5jzP//3fAs45p0kIoowcubnQvbudPXs0AFassLBoUX6IoxLlITUFIYCXX7ZRVGTWEDIyVBYssJ3x\nfTEx0Xzyyec0bNgomOGFnSNHFI4cOfvfv/9e8yQEgDVrLGRlBSEw4TdJCiJs7NunsHmzRn4ILijj\n4krWChISTq8lAAwZch9xcWVPfl6ZTZtm44IL4mjWLJbx48/8njp13Gha8XdYvbo7IuZaFn40H6mq\ngqaZOeXEz0gmZQmtJUssPPSQDV1XOP98nRUrCkhIMP8WjPI884yTX37R2L1b5eqrdR580FXqJDKK\noqAoitcTzZxcloqYpKaiHTigMGtWlOf15MnQr59GnTol39e0Kbz8ciHTp9uIiTGYPr2IqKjwLW8k\nHjNn429ZFMMwznxJVAbDMFCU0m/ICVFejRrBvn3Fr+fOhaFDg/PZug7PPw+7dsENN0DPnmUvc9dd\nd1FYWMiiRYu8+qzs7GwSExPJysoi4UTWiyD79pnb6mS//AKNG4ckHFEBfK4pHD/uwGrVSEiIITs7\nH12P7Am2NU2VsoSQ1RrDya2Zul5IZqY5SXxFl2fcOBsvvmgF4LXXDJYtK+Dqq0v/nMJCF0VFLjIz\nHV59lsNhto2ZZdHKeHf4SUyEoUNtvPqq+X09/DDUrJlPZmZk7GdnE4nHzNmUVpakpNgyl/c5Kbjd\nhucDdd2NyxXZX+QJUpbQmDatgDvvjCE7W+G661z07FmEy1XyPRVVni+/LE5GhqGwbp3KFVe4SlnC\nrCkbhuF1PJXhmJkypYA77yxEVVUuv9xOZmbkluVUkbxdTuVrWaRLqggL7dvrpKfnkpsL1asH97Mv\nukjnhx+Kr9ovvLBynBQqUpMmBhaLTy3PIsxJUhBhw2YLfkIAeOaZQuLiYM8elW7dXHTvfvZaQlpa\nLRRFwel0ArBy5ccoisJvv/0ZrHCFqFCSFESVFxsLU6YUluu9v/9eSud8ISqByO9/JYQQImAqbVJ4\n4w0rt98ew9NPR1FQEOpohBAiMlS65iOHAwYNimH9erNoa9dCQQFMnVq+5gEhhKjKApoUfv1VYeVK\nC3XqGPTq5SIUz7ZNmRLlSQgnfPdd5PUHF6Gxf7/C559baNDATYcOeqjDESLoApYUfv9d4frrY/nr\nLzMTbN9exKRJwb86//nn01vE2rcvvc95VeR2w2uvWfnpJ5Vrr9VL7XFzNnv2KHzzjcb557srRTfO\nffvMfTgz09yHx40r5OGHfRvzOS8PPvjAPLx693YRExOwMIWoUAG7p7B2rcWTEACWLQtNy1TXrief\n3AwGDChi1CgZzP1U06bZGDcumoULbdx1VwyffupdbWrbNpXOnWMZPjyG666zs3Jl5LdErlhh8SQE\ngIULrT6tx+WCW2+N4ZFHzH/9+sWgS6VDRIiAHckNGpS8UqxfPzQPttx7r5PkZIMffjAHNrv66spx\nNB4+rDBjho38fIUHHijy+8r8q68sp73u0qX839XixVby880TqK4rvPWWlawsKChQ6NvXSQQO60Pt\n2sYpr337jvfuVdmypfj73bTJwr59yhnnZxAi3AQsKXTqpPPUU4UsWmQlNdXNc8+FrstPr14uevUK\n2ccHnGHALbfEsGuXeTX/6acWNm50kJzs+0nm4ot1tm7VSrz2Rq1aJT/7p59URoww20jefNPK6tV5\nEddk0qePi2+/LWLZMgtpaQazZ/u2D9eoYRAdbVBQYCbNmBiD6tUlIYjIENA6//DhRQwfLk01gfbX\nX3gSAkBWlsKuXSrJyb7XgsaPLyQqCnbtUunc2UWfPt7dUxg2rIgff1RZv95C8+Z6iSvjnTs1duxQ\nadMmsu4zKIr5dPMzz/h3L6xGDYNXXilg4sQoFAUmTSogKSlAQQpRwSK/IbgKqFYNzj1X5+efzcSQ\nmGhw7rn+nXCjo2HCBN9PfnY7vPmmeSVdVATnnx9HVpZ5ZWyxGKSkVO0r4xtucHHDDdLBQUSeSvvw\nWmWiKLB0aT4DBxZx881O3n8/77Tmm1Cy2eCtt/Jp1kynYUM3L7xQQIMG4ROfEKL8pKYQIVJTDWbN\nCt8H8K64QmfdurxQhyGE8JPUFAQA8+ZZ6ds3hieeiCI3N9TRCCFCRWoKgo8/tjB6dDQA69dDfr7C\nCy/41vPm++9Vtm/XuPRSPawfaHM64aOPLDid0KOHPFwmxAmSFAQ//KCW+rq81qzRGDgwBl1XsFgM\nFi7M55prwu85EcMwx8dau9bc/d96S2f58jysvj2rJkSlIs1Hgg4ddFS1+MZwx46+ncgXL7ai62YP\nJJdL4Z13wvMse+iQ4kkIAN98o/HTT3IoVEW//abwxx8hGKQtjMmREKa2blVZutTC4cMVv8O2a6ez\neHE+Q4YU8eyzBYwb59sN7VOfCE5NDc/mo8REA7u9OFaLxaBGDektVdWMHBlF69ZxXHJJLM8/H54X\nMKEgSSEMvfGGle7dYxk2LIZOnez89tvZE0NuLgwZEk2rVrE89FB0ibkj/vhDYcyYKB5/PIq9e0tP\nLh076kydWsg99zhRfdwrnnyykK5dndSs6eaGG5w89lh4PsgYFwevvppPWpqb1FQ3//53AXXqSFKo\nSnbsUFmwwAaAYShMmmQlOzvEQYUJuacQhubPL75qOXZMZfly61mfFJ86NYqPPjLf/847KmlpbsaM\nceFyQa9e0ezebZ7hV60yh8aoyDGJ4uPhP/8J/xmN/vlPG6+8YiMpyXzyuF27ks1lBw8q/PijSosW\nburVOz1Z/Pmnwtdfa6SluWnVKjxrQ8J7hlwXAFJTCEs1a5bcO5OTz37i2b+/ZA3gwAFzkx4+jCch\nABw5orJnz+mbe80ajVmzbGzaVDXmnPj6a43Zs6PIy1M4eFDlvvuiS/x92zaV9u1jGTDAzlVXxbJt\nW8nvbP9+hU6d7AwdGkPXrrG8+aY0O0SiCy90c8cdxRdaY8c6SUwMYUBhpMonhT17FDZu1Mir4Oeu\n8vLgyy810tPL/spnziygRQsdu93g9tud3HrrmYdLWL1a+7upx0wimmbQu7cTgJQUqFevOJlUq2bQ\nsGHJ5LJokYX+/e1MmxZFr14xfPZZ5U8MGRklk+jx40qJK8TXX7fhcJjvcTgUXnvNVuL9H3xg5ejR\n4m04b174J4UPP7Tw8stWfv1Vbqie7LnnCtm8OZfvvsvlscecoQ4nbPjUfGQYBg5HDhaLhqbpOBz5\n6HrkVaMXL7YwerQNw3DTrJmbTz/NBgJfltxcuPnmaH76yTz7jB9fyN13n31cnNRUWLkyx/O68Az3\nfd9918I//hEFmOu5+WYnAwa4uOQSNw6HiqbpLFpUwMyZFoqKYNgwJ9HR7hLJb9myKMA8GNxuWL7c\nSdu24XcfQNPUgO1nrVvDOecU8ssv5ol90CAn+fnFZbbbiwDrSa+d5OUV/z0uzgJEeV4nJOjk5Z25\nySwzE9av10hJMbj8cjPuvDxHiZ8VbepUK3Pnmvvd7NkGH3+cH7Bh7QO5XY4fh9xchbS04LbhpKSY\nP08cM5F6LjtZadslO9sgPj4epZRpMRXD8L4lLTs7m0SpawkhRMTJysoioZSbiz4lBcMw+O23P7BY\nNBISYsjOjszses01Mfz668lNAXDddYEvy3//qzF0aHHbdd26bjZsyPdrndOmWXn55eKmjccfL2LY\nMPOqX9PUcm0XhwOeeiqKHTtUrrxSZ/z4IrSTWpAKCqBDhxj+/NP8jmrUMPjiizzi4/0K3WvlLU8g\nOZ14HmZzOmHKFBvffKNy6aVuxo0rKvNBtyVLLDz5ZHGNokYNg23b8sjLc9CiRVPS03djt8dWYAlM\nfftGs21b8UadP7+ATp0C80BhILaLYUDz5nbP3BMAb7+dT/v2wT2fhGIf81dBAVx1VYynOTMpyWDd\nujyqVTt7WZKSYsusKfjUfKQoCrGx8VgsKgkJsei6hssVGV/kyWbN0rjrrhiyshRuusnF4MEWsrMD\nX5aePWHt2miWLLGSmGjw/PP52O3+td+PHg35+VF8+61Gu3Y6jzyiYLWaJ6Hybhe7HebOPfFKA0q2\nnx8+rPDnn3Ge1xkZkJFhISUluNs61PvZjBk23nrL/G5//BFq1iwsc4rX1FQNsHteJyfr2O3FFyB2\neyx2e9wZlgysl19WeOSRaA4eVLn9dic33hi4+wqB2C5uNxhGHFAcl81mwW4P7pPwod7HfHH4sMLR\no8X7UGYmHDtmoW5dzlqWhISyr+iqdJfU9u11du7MJS8PqldX0bSK+ToUBV54oYBp0wqIjsbr5wAK\nC2HfPpXUVLenS2l0NEyfXrGjptata5CW5ub3382A69RxnzbtalVwaq+tM/XiOlXXrjpDhxbx9ttW\nUlIMn8eS8lfDhgYffOBfrdRbf/6p4HRyxu68p1JVePrpQsaMicLtVujWzVlpptCtaPXqGTRu7Gbv\nXvXv1+6/O5P413+oyvc+slgI2nzCdrv3CeHYMbML5FVXxXLppXFs3Rq8TRYVBcuW5XHnnUUMGlTE\nsmV5xP7d4vHHHwo5OaUvX1l07VqyU0C3bi62blV59VUr27effXtMmVLIvn25bN7s4OKLq0Yyffll\nKxddFEurVnGMHBlV9gLA3Xc72b7dwYYNDt58s8DnhyfDnWEQ0GPGZjOPz6FDi7jnniI++KD4+PSH\nT/cUAI4ezcFiUUlKiiUz0xExVa6zCdeyTJ9uY+bM4oPriitcZV75VWRZ3G4YOjSaDz+0EhVlMGdO\nAT17VuwMY+Gwbdau1diyReOyy3SKihSGDInG7VbQNHPgv/K20+fl5dKwYR327TsUlOajimAYZlfc\nzZstdOhgYfBgB7ruJjcXmjSJw+0ubgpas8bBRReFz/F0NhW9j+3apXLbbTEcPKjStq2LRYvyiaug\nzV9aWZKTy24+qqQ5ufLavFnj889D9zzB2rUaH35o3mUtLFR44onoMpaoHDp31hk9uoguXXTefdfi\nOfHpusJ774X/swqBNH++lTFjolm+3MJjj8FLLxU3u556ialLSxAAEyZEcfCgebrdtMnC66/bylgi\ndCQphLkhQ5wlBpbTdYXHHw/didjpVE55HaJAQujUtvKTHxKsCrZuLXlRsnmz+TouDkaPLr4Bf+ut\nTlq29P67cbnMmsjTT0fx3XeV4xR16sOxFf2wrD8qxzdeidWoYfDEEyVvKBeEcHiha691cdVVZnOR\nqhqMHx++U4RWlCefLOTGG53UqeOmTx8njzwSfg/8VaTWrUte/rdtW/z6kUeK2LYtl40bc5kzx7cd\nddSoKEaPjubFF2306GFn587IP00NH15EVJR5MVG7tpuBA8P3aqpK9z6KFD17upg/X2fHDg1VNcrs\nDlmRbDZYsiSfnTtVqlUzAvZ0bEVyueD33xWSk42APGMRF2f296+qhgwxT2hbtpj3FAYOdJVoJvJ3\nn/j00+LTUkGBwldfaTRvHtm1sWuv1dmwwcG+fSoXXaRTrVqoIzo7SQoRIC4OVqzI4/vvNZKT3TRu\nHNoTscVCWE+1ebKcHOjb185332nExxssWJDPFVdIQ7e/hgxxct99OklJFjIzA7vuZs3c/PFHce3g\nvPMiY18rS1qaQVpa+O97ftfLXnzxRdq0aUnjxnW58srLmDdvbtkLhak//viDAQNuo2HDVC64oCnP\nPPN0qEPyiI6Gyy/XvU4IDoeDli1bMGLEgxUUWcX76KPlXHLJJaSl1ebKKy/j7bffKveyCxZY+e47\ns807J0cctoW6AAAVXUlEQVRh0qTydZMUpTtwYD+DBt1OzZo1ad78HIYPf4CcnMBMSDBnTgF9+jhp\n08bFzJkFQXtu4dFHH6VGjSA/rl8BZs2aTp06dUhLq80tt/Rk//7fvVrer6Tw6aef8OSTTzJ37jz2\n7j3ICy+8wuTJE1m79r/+rDZkevfuTYMGDUlP/4Xly1fy1VdfsmHD+lCH5Zdp057B4cgNdRheMQxz\nEEGA7du3cf/99zBlyhT27TvEpEnPMmrUSLZs2Vyudbndpb8WvhkwoB/VqlVn//79fPbZV+zatZOJ\nE58KyLqTk815Lj7+OJ9Bg4LT9r5jx/9YsGBBqcM/RIJ5817l3XffYd26daSn7+G885oxd+6LXq3D\nr6Tw/fffceGFF9KyZSsAWrW6jObNW7Bjx//8WW1IbNiwnl9//ZWnn56C3W7nnHOasmrVZ1x55VWh\nDs1nP/74Ax988B79+t0R6lDKbedOlZYtY2ncOJ7evWP4449sHnvscW688UZUVaVz5+uoX38Qb7xx\nvFwPAg0c6KRFC/NK0243GDu26t0YD7Ts7CxatmzF+PETiYmJITU1lVtv7c/XX28IdWg+MQyDkSMf\nYeTIkaEOxW+vvDKHceMm0qRJE+Li4pgyZRpTpkzzah1+JYVOna4lPT2dDRvW43Q62bp1M3v27KZT\np2v9WW1IbN68iQsvvJDJkyfSvHkj2rS5mJdfnhPqsPwycuQjjBkzvtQREcPNuHFRHDpk7pYbNljY\nu/c6Ro58wvP3KVOs/PLLa7z33q1062YvMzEkJsJ//5vHl186+PbbXDp0CP823XCXkJDIc8/NoWbN\nZM/vDh48QGpqnRBG5bs335xHTEw0/fv3D3Uofjl8+A9+//03jh8/zvnnn0+TJmkMGTKIjIwMr9bj\nV1Jo1epSZs2aRZ8+PahfP5nevbszatRTXHTRJf6sNiQOHTrIxo0bqVUrhe3bd/LPf85k6tRJrF69\nMtSh+WTu3Llomsptt0VOLQEgL69k9f3EhDcnvPRS8S77888aX3xRdl8Jmw2aN3dTvXpgYhQlbd/+\nLfPnv8qjjz4e6lC8duTIEWbMmMrMmf8KdSh+O3ToIAAffvgBn332GevXb+bQoYOMHDncq/X4lRTW\nr/+S0aNH8957y/n99yO8//4KnntuRkSeSA3DoFatWgwb9jDR0dF06tSFG264keXL3w91aF47evQI\nEyZMYNas50MditdGjCjEZjNvpqemuhkwoLhNeeLEcbjdx0q8Pykp/LvEVgavvmqlc2c7gwZFc/hw\ncaLesGEDN9/ci3HjJtG+/dUhjNA3EyaM4Y47BtG06bmhDsVvJ0YsGjHiUVJSUkhNTeWJJ8bwyScr\nKSoqfzd2r5LC22+/TUxMDHa7nQYNUnjjjXn07duX9u2vxmaz0abN5fTufTOLFi3wrjQhsHTpYtLS\natGgQQoNGqSQklKbaqd0Hq5fvwFHjhwJUYTld2pZxo8fy+DBg2nWrHmoQ/NadvZCoDlW6w1kZtaj\nbl0DwzAYPHgwn3yyirlzi0hNdRMdbfDQQ4W0by/NQRXt8881nnoqmh07NFavtjJsmPlE/erVK+ne\nvTv//OcMhgwZ6tU63W549lkb3brZefzxKPKDO5ArAOvWfcHWrVt47DGzedLHYeDCRq1a5jRyCQnF\nE6ClpaVhGAbHjh0t93q8ek5hwIABDBgwADAHxBsyZCD6KYObFBVFxo28W265jVtuuc3z+uOPP2DO\nnOfJy8vDZjN3+v37f6N+/fqhCrHcTi1LSkoiSUlJzJs3D4D8/Hzcbjeffrqa9PS9oQqzXE4uS24u\nvPGGlY8+ep+8vL188sla7PZ4evYMzlSWZ7JmjcaiRVZq1zZ48snCKjHZ+4mpS0/Ys0dly5bNDBt2\nP++99x6tWrX1ehC5efOs/OtfZvfgbds0oqNh8uTgnjvee+8djh07SqtWLQAzKRiGQYsWjZk6dSY9\ne/YJajz+qlOnLvHxCezY8T+uvrodAL/99htWq5XatVPLvR6/Hl7r1q07Y8Y8QZ8+/WjZ8jL+97/v\n+PDDZUyc+Iw/qw2Jrl1vICkpiQkTxjJu3GS2bdvK6tUrWbLkg1CH5rUffviZatXs/PVXHrru5qWX\nXuDw4UNMnvzPUIdWbk4n3HyznW+/1YD+tG59C7GxRRhG6PqU7tihMmhQDC6X2Xyyb5/KokUhuMQN\nso4dXdjthud+T7duTkaOfJgJEybRuXNnMjO9T9I//VQy0ezaFfyhLCZPnsro0eMAc+a1nJwM2rVr\nx+efbyQxMYwfOT4LTdPo338gs2fPoGvXa9F1ldmzp3PLLbehejEeuV9J4bbb+uNyFTBixIMcOvQH\nqampjBgxkn79Iu8ufnR0NKtXr2bIkHto3rwRNWsmM2PGv7j88rahDs0rhgF5eXWoXt1Oaqo5dG58\nfDxZWXZSUmqHOrxy++UX9e+EYNq61Uq9eq1Q1V8AaNv2CpYsWRbUmL7/XvMkBDCvcHfvVhk5MoqM\nDIUhQ5zcfXdox7TZt09h1SoLqakGPXu6CES3+yZNDBYtymPUqGjy8hTc7n3s3v0zo0c/zqhR/0BR\nFAzDQFEUNm7cRt269cpc53XXuViwwFbidbAlJCR6mlosFpXYWCuKokTUcXKqp56aiK47adOmDU6n\ni5tu6ul1l1SZT+FvlaEsug533RXN6tVWNA1mzSqkf//IHKwtI0OhVatY8vPNs1psLPzwg4PY2NBt\nm/R0lS5d7J6RYi+77A8cjlrs3FmcvFascNC69dljnD+/gFGjajFnzp/cemtMQOPbv1/h2mtjycw0\n47vvvqKANckMHhzNqlXFQ4S/9VY+N93k9uuY+ewzjfXrLVx0kU7v3sFPCierDMf/CTKfgvD48kvz\nRiCYCWLsWNtp49tHiho1DF5/PZ9zz9Vp1szNu+8S8vb7hg1zadjwQRo2/AqYzsMPf8O+fSUPoRNT\nl57JK69YGTXKbEd/6KFoli0L7NBjn31m8SQEgPffD9z6f/655HDZu3f7f+ro1ElnwoTCkCcEUZIk\nBRG2unTR+eqrPDZuzKdr15J/c7nMXjFffx28CYccDgf339+CLVsuRlFGER2t06tX8QmtZk13qb2h\nvvyy5El65crAJoX69UteFaalBe6K4OQpSW02g44d5UReWUlSqEQ6dNDp1s1s09Y0ePbZooC0KQdL\nQYHZJl5Wl2pdh9tvj6FfPzs9e9rLPRewv5KTkxk48M4Sv5s9u4B//zuf8eML+OSTPFJSzn4ivuii\nkglj5UoLmzYFLql16qQzdmwhjRq5ueIKFy+95P9N8B9/VLnmGjtLl1ro2tXJiBGFfPhhXokpNo8d\ngzvuiOHSS2MZNSpKZluLcHJP4W+VpSyGAfv3a9SrZ8dmi5yy/Pyzys03x3D4sErjxm6WLcsjNdXc\nNU/dNtu2qXTrVnKG8l27ckhKCl68KSmJLF78Ptdc07ncyzid0KWLi/T06kAWkEDXrk7+85/wnZvh\nqqvs7NpVnLjefTePq6/WOXZMIStLpXVrO/36uXjvveJaz9SpBZ45FyJFZTn+wf97Cj7XX1VVQdPM\nisaJn5GsMpWlaVOFhATIzo6cssyYEcXhw2a8e/eqvPRSFFOnmlWGU7dNtWolqz82m0FsrIolgK0x\nBw4o3HNPFD//rNKtm4vnny86bf2apmKxeNHVzwLt27tJTy/+XXy84tU6gu3keQ0AjhzRWLlS5b77\noigsVLjqKnMinJMdOuTd9xIOKtPx73dZDB+53W5fFxXiNH36GIZZzzH/DRtW+vufecYwLBbDsNsN\nY+HCwMfTo0fJeKzWR0r8XVEU45NPPvF6vb/+mmUABmQZ551nGPv2BSriivHYY8XfQd26hnH4sGE0\naFDyuxk0qPj/UVGGsXlzqKMW/vD52ur4cQdWq0ZCQgzZ2fnoemRXuTRNlbKE0PDhKl98Ec3x4wp1\n67q5554CMjPN5qMzleeBB+Dee0FVzX+Bnv3rwIFooLjZZPjwaac9pJWTU+D1g1sWi9nOv3NnPikp\n5voDHXsgPfUUXHqpxpEjCl276thsBoYRw8m3I6+4ooibbnLz008KHTvqNG1qhHWZziQSj5mzKa0s\nSUmxZ1mqmM9Jwe02PB+o6+6Ib4c7QcoSGi1auNm8OZfff1dp1MhNXJzZw+hkZyqP2x34iXPcbqhX\nz+15eC4uzqBnT+dpn+3L93vimImJiZxt06VLcZwuF0yaVMj990dTUKDQsSP07OlE09x06FD8nkgV\nScdMWXwtS+Q3oIlKIzHRnPs5Li60cTz/vI0PPyx+UOuxxwpp1sxdYuBBRVEYNOg2GjRIYeTIESGM\nNvhuuMHF99/nsnVrHmvWQJTMcFqpBLajtBCVwObNJbuJ7t1rXjudOvBgVZaUZE6bqfnZozY3F+x2\nswlQhAfZFEKc4tJLS3a0v+wy6XgfaE4nDBoUTePG8Zx/fixbt8qpKFxITUGIUzz2WBE2G3z/vUr7\n9jr9+weukfybb8yT35o1Gj16BGy1Eefddy2eIVkyMszRZ/v3d/Lgg05q1IjQsVkqCUkKQpxC02DE\niMAPJPjllxr9+plzddxzTzQzZlgZPDiyHvIKlFOnXc3IUHnhhSgWL7byv/850DQ4dEjhuedsuFzw\n4INOmjatHDeAw53U2UTQHDyosGOHirNqngdZudKC2118Mvzoo6p7Tda7t5Nzzjn9JH/0qMqaNRou\nF/TpY+ett2wsXGijV68YsrJCEGgVJElBBMWSJRYuuyyWzp1j6d07hoLwHdmhwjRuXPIk2KhR1b3y\nrV4dPv3Uwbx5+UDJ5qKCAoUjRxTPDX4wk8WpM8CJiiHfsgiKSZOi0HXzKnnLFgsff1z1rpLvuad4\nEp5evVxMmBAZU9dWlLg4uOkmF927F9+zqVfPzTXXuEhONkhLK06a1au7q3QSDaaqd2SKkLBaS39d\nFWgajB9fxPz58K9/FWK3V8Ev4QzeeKOA1aud/PWXwvXXu0hIMH//7rt5zJwZhdMJDz9cRFKSOZjG\nvHlWNm/WuOwynaFDnRE1EnAkkKQggmLqVPMp2Px8hS5dSl4dCtG16+ndfhs2NJgzp2Q74/z5VsaM\nMW/WL19uJtX77quiN6kqiCQFERTdurn44YdccnIUUlMNuboTPtmypeTTcps3a5IUAkzuKYigiY+H\nOnWClxAyMhQOHpTs4y1dh7y8UEdxZm3a6KW+Fv6TpCAqpTfesHLBBbG0bBnHiBHRoQ4nYnz2mUbT\npnE0bBjPiBHRYTfH9913O3n22QJ69HAyeXKB1BIqgCQFUekUFcHYscW9nf7v/6xs2SK7enmMGBFN\nbm7x97Z2bfDmwC4PRTF7cb3+upkQpBky8ORIEZWOYZw+nPaJBCFKd+qTxg5HxX1vhgF79yr8+ads\nm3AiSUFUOlFR8NRThSiK2fZx001O2raVtufyGDmy+NmJiy/W6dKlYnqJud3mUB9t28Zx8cWxvPGG\ndM8NF9L7SFRKDz3kpEcPF3l5Cued55ZmhnJ68EEnHTvqZGQotG6tE11Bt2O++krjo4/MROB2K4wf\nH8XgwU4ZQjsMSFIQlVZa2t/TIQuvtGhR8U8OS5IOX5KXhfDCxx9/yDXXXEmjRnW48srLePvtt0Id\nkleKimDcuCi6d7czZYoNPUStaldeqdOjh9lzSNMMJk8ulFpCmJCaghDltH37NoYNu5fXXnuTa6+9\nns8/X8Pgwf0599xmtGlzeajDK5fZs23MnWsDYOtWjWrVDB56KPjdOlUVXnutgDFjComNhZQUqdGF\nC8nNQpRTZmYmjzzyD667rhuqqtK583W0aHE+mzZtCHVo5ZaeXvKQ/+mn0HU5VRRo3NiQhBBmJCkI\nUU6dOl3Lo48+7nmt6zp//vkntWunhjAq73TpUrK9qHNnGYNKlCTNR0L46OmnxxEbG0uvXn1DHUq5\nDRzoJCHBYNs2jXbtdLp1k6QgSvI5KaiqgqaZFY0TPyOZlCV8hWN5Jk4cx/Ll7/PRR6uw28vfb/Pk\nslgsoSlP375u+vY90cPI9xjCcbv4SspSTDEM30Y3MQwDRfqViUrs7bff5t577wVAURTy8vIwDIM7\n77yTb775hlWrVpGWlubVOrOzs0lMTCQrK4uEExMHCBFGfE4KGRm5WK0aCQkxZGfno+uRPSuSpqlS\nljAVTuUZNeoffPPNVt5//0MSEhK9Xt7hyKV+/drs33+Y2Ni4CogweMJpu/irqpQlKSm2zOV9bj5y\nuw3PB+q6G5crsr/IE6Qs4SvU5dm8eRNLly5hw4Zt2O3xPsUix0x4k7LIjWYhym3x4rfJycnh0kvP\nL/H7tm2vYMmSZSGKSojAkqQgRDk999wcnntuTqjDEKJCRf6tdiGEEAEjSUEIIYSHJAUhhPCBywUT\nJ0Zx/fV2Ro+OorCw7GUigdxTEEIIH7z0ko2XXjIHF9y+XSMuzmDs2KIQR+U/qSkIIYQPdu0qefr8\n+efKcTqtHKUQQogg69q15LhR111XOaZ8leYjIYTwwU03uVi4MI9NmzRatnRz442VY3BBSQpCCOGj\nLl3004Yjj3SSFIQQohLwbRS708k9BSGEiHCLFllo1CiOBg3imD/fv2t9SQpCCL/t2KEydGg0Dz4Y\nza+/ypD6wZSRofCPf0STl6dQUKDwxBM2Dh70fX3SfCSE8EtmJtx8s53MTDMZbN6ssWmTA6s1xIFV\ngIICuO++aNautdC8uZs338ynbt3QzjGdmwsuV3EidrsVsrPBbvdtfVJTEEL4Ze9e1ZMQAPbvVzl6\ntHLWFl591caqVVaKihS+/15j3Lgoz98OHlR48MFoBg6MYf16LWgxNWhg0KOH0/O6a1cXzZr5vj6p\nKQgh/NK0qZvkZDdHj5rXmOec46ZWrdBePVeUjIySye748eLX/fvHsHOnmQzWrdNYt85BgwbB+R5e\nfbWAAQOcuN3QubOBovh+apeaghDCLwkJsHx5HgMGFHHXXUW8914elkp6uXnbbU7i480TvaYZ3H23\neYVeVIQnIQDk5ytBfcJZVaFjR51OnXQ0PysplXTTCSGCqUkTg9mzK8mIcKVo3tzNF1842LpV49xz\n3VxwgTmzmc0GrVvrbN1qnpETEw0uvDAyZ3CTpCCEEF6oX9+gfv3Tn15etCiPf//bRna2wp13Oqld\nOzKb0CQpCCFEACQmwrhxMkqqEEKISkQxjEA9HC2EKEt2djaJiYlkZWWRkJAQ6nCEOI0kBSGCyDAM\ncnJyiI+PR1EqZ19+EdkkKQghhPCQewpCCCE8JCkIIYTwkKQghBDCQ5KCEEIID0kKQgghPCQpCCGE\n8JCkIIQQwuP/AQGX0Ke8jtXDAAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list_plot(Z.numpy(), figsize=4, zorder=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

異常値の検出

\n", "\t

\n", "\t\t異常値は、主成分分析を行った固有値ベクトルに投影した値と観測値の残差と井出本では定義しています。\n", "\t\tこれを再構成残差と呼んでいます。\n", "$$\n", "\t\ta_1(x') = (x' - \\hat{\\mu})^T \\left [ I_M -U_m U_m^T \\right ] (x' -\\hat{\\mu})\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tSageでの計算では、numpyのarrayにした後、カラムの和(Rのcolsumsの代わり)を求めるため、sum(1)を使っています。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 要素毎の積を計算するためnumpyデータに変換する\n", "xcm = Xcm.T.numpy()\n", "x2 = (Xcm.T*(U12.T)).numpy()\n", "# 要素の二乗を計算し、カラムの和を取り異常度a1を求める\n", "a1 = (xcm*xcm).sum(1) - (x2*x2).sum(1)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAECCAYAAAARlssoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xd4FNX6B/DvzOxushuyIZFQEkpQBAHBQOiXq1J+YAFE\n76VIsWEQBAFBhYsocGkiCggoReQi0jtRmkpUutSASlUCCcVUSMK27M6c3x8Ds2x6wmxL3s/z8IRt\nmZN3Z+adU+YcjjHGQAghhADgvV0AQgghvoOSAiGEEAUlBUIIIQpKCoQQQhSUFAghhCgoKRBCCFFQ\nUiCEEKKgpEAIIUThM0mBMYbs7GzQvXSEEOI9Gk9vMC0tp8DnTaYc1K0bicTEawgKCvZwqfwHz3MI\nCwtCZqYJkkQJtCAUo+JRjIpX3mIUHl6y86rHk0JhOI5z+VkRiSKwfr0GGRkcnn/egcjI/Dsiz3Pg\nOA48z5WLHdUdKEbFoxgVr6LGyGeSAgFGjAjEhg1aAMCiRRLi482oWrXi7IyEEO/zmT4FAmzb5szR\nqak8Dh0SvFgaQkhFREnBh9SpIyn/5ziG2rWlIt5NCCHqo6TgQ5Yts6JNGwfq1xfx8cc2NGtGSYEQ\n4lnUp+BDGjSQEBdn8XYxCCEVGNUUCCGEKCgpEEIIUVBSIIQQoqCkQAghREFJgRBCiIKSAiGEEAUl\nBUIIIQpKCoQQQhSUFAghhCgoKRBCCFFQUiCEEKKgpEAIIURBSYEQQoiCkgIhhBBFqZPC7t27Ub16\ndfTr1y/fa7/88gvatWuHkJAQPPTQQ5g2bZoqhSSEEOIZpUoKs2bNwqhRo1C/fv18ryUnJ6Nbt254\n9dVXkZmZibVr1+KTTz7B6tWrVSssIYQQ9ypVUtDr9Thy5AgeeuihfK+lpKQgNjYWsbGxEAQBLVu2\nROfOnbF3717VCksIIcS9SrXy2vDhwwt9rUWLFmjRooXLc8nJyWjatGnZSkYIIcTj3NbRPH/+fFy6\ndAlDhgxx1yYIIYSozC1rNC9YsAATJ07Ejh07EB4e7vIaz3PgeS7fZwSBV35qNDQoqjD3xokUjGJU\nPIpR8SpqjFRPChMmTMDy5cvx888/F9h0FBYWBI4rKCmIAACjUQ+jMUjtYpU7RqPe20XweUajHmvX\nAkePAk88AfTo4e0S+R7aj4pX0WKkalKYPXs21q5di8OHD6NmzZoFvicz01RgTcFksgAAsrMtEEVB\nzWKVK4LAw2jU34mT5O3i+KS7MZozJxejR+sAALNnA8uWWdGzp+jl0vkG2o+KV95iFBpasott1ZLC\npUuXMGnSpCITAgBIEoMksXzP3w26KEpwOPz/C3A3ilPxvv+ez/NYQLdudi+VxjfRflS8ihajUiUF\nvV4PjuNgt8sH1pYtW8BxHMxmM1avXg2z2ewyAokxhqioKJw9e1bdUhNSAg0bSti1y/n4kUeolkBI\ncUqVFCwWS6GvTZgwARMmTLjvAhGilvfes8NsBk6e5NGunYg33qBaAiHFccvoI0J8QUAAMGWKzdvF\nIMSvVKyxVoQQQopESYEQQoiCkgIhhBAFJQVCCCEKSgrEr0gSsGCBFoMGBWLpUq23i0NIuUOjj4hf\nWbBAh6lTAwAA336rhUYDvPIKDTUlRC1UUyB+5dgx11326FGaEoUQNVFSIH6lVSvXu5LbtKG7lAlR\nEzUfEb8ybJgdWi1w8qSANm1EDBxITUeEqImSAvErHIc701VQMiDEHaj5iBBCiIKSAiGEEAUlBULK\nodu3gR9+EPDbb3SIk9KhPYaQcubWLaBLFwP69zegU6cgLFlCN/mRkqOkQEg5s2OHBn/+6bx/Y948\nnRdLQ/wNJQVCypngYNfHRmP+5W8JKQwlBULKmW7dHOjb1w6OY6hSRcKcObTQECk5uk+BkHKG44B5\n86z49FNAS90JpJRKXVPYvXs3qlevjn79+uV7LT4+Hq1bt0ZISAiaNGmC1atXq1JIQkjpUUIgZVGq\npDBr1iyMGjUK9evXz/fa33//jeeeew5vvvkm0tLSMHfuXMTGxuLEiROqFZYQQoh7lSop6PV6HDly\nBA899FC+11atWoUGDRrg5Zdfhk6nQ6dOndCjRw8sXbpUtcISQghxr1IlheHDhyM479CGO44fP47m\nzZu7PNe8eXMcPXq07KUjhJAKJDsb+PlnAX/+yXmtDKp1NGdkZKBWrVouz4WFhSE9Pd3lOZ7nwPP5\n/2BB4JWfGg0NiirMvXEiBaMYFY9iVDxPxyg1FejaVY8rV3gIAsMXX9jQq5fnp4ZXdfQRY8WPhw4L\nCwLHFZQU5D/eaNTDaAxSs1jlktGo93YRfB7FqHgUo+J5KkbLlgFXrsj/F0UOc+cGYvBgj2zahWpJ\nITw8HBkZGS7PZWRkoGrVqi7PZWaaCqwpmEwWAEB2tgWiSKtpFUYQeBiN+jtxkrxdHJ9EMSoexah4\nno4Rz2sABCiP9XoRN29aVfv9oaElu9hWLSm0aNECy5cvd3nu6NGjaN26tctzksQgSflrFHeDLooS\nHA7aSYtDcSoexah4FKPieSpGvXvnYvduAbt3a1C1qoSPPrJ65btRrbGsf//+uHz5MpYtWwabzYYd\nO3Zg586deOONN9TaBCGElFs6HfDNNxYkJubgt99MiI72TrIuVU1Br9eD4zjY7fKqV1u2bAHHcTCb\nzQgPD8d3332Ht956C8OGDUNUVBRWrVqFxo0bu6XghBBSHgV5uUu1VEnBYrEU+Xr79u1x8uTJ+yoQ\nIYQQ76HxaIQQQhSUFAghhCgoKRBCCFFQUiCEEKKgpEAIIURBSYEQQoiCkgIhhBAFJQVCCCEKSgqE\nEEIUlBQIIYQoKCmQfBgDHA5vl4IQ4g2UFIiLPXsE1K9fCTVrVsIHHwQU/wFCSLlCSYG4GDYsEFlZ\nHCSJw+LFOhw8SAseEVKRUFIgCkkCbt92XRUvO9tLhSGEeAUlBaLgeWD48FzlcXS0iCee8PzC4YQQ\n71FtOU5SPowbl4v/+z8HsrI4tG0rQk/ruhNSoVBSIPnExNCavYRUVNR8RAghRKFqUkhISECnTp0Q\nGhqKiIgIDBw4EOnp6WpughBCiBuplhREUcSzzz6Ldu3aIS0tDX/88QdSU1MxbNgwtTZBCCHEzVRL\nCjdu3MCNGzcwYMAAaDQahIaG4oUXXsDJkyfV2gQhhBA3Uy0pREZGolmzZliyZAlMJhNSU1OxadMm\ndO/eXa1NEEIIcTPVkgLHcdi4cSO2bt0Ko9GIGjVqQBRFTJ8+Xa1NEEIIcTPVhqTm5uaie/fu6NOn\nD8aPH4/bt29j6NCh6NevHzZt2qS8j+c58DyX7/OCwCs/NRoaFFWYe+NECkYxKh7FqHgVNUYcY4yp\n8Yt27tyJXr164fbt28pzp0+fRnR0NDIzM1G5cmUAAGMMHJc/KWRnZyMkJARZWVkwGo1qFIkQQkgp\nqVZTEEURkiRBkiTwvJxZrVZrvgSQmWkqsKZgMlkAANnZFogiTcJWGEHgYTTq78SJbjIrCMWoeBSj\n4pW3GIWGBpXofaolhXbt2qFSpUqYOHEixo8fD7PZjOnTp+OJJ55QagkAIEkMkpS/cnI36KIoweHw\n/y/A3ShOxaMYFY9iVLyKFiPVGsvCwsKwe/duHDhwADVr1kSTJk1gMBiwevVqtTZBCCHEzVSd+6hZ\ns2aIj49X81cSQgjxoIrVrU4IIaRIlBQIIYQoKCkQQghRUFIghBCioKRACCFEQUmB+IQ//uCxZo0G\n58/TLkmIN9FynMTr4uMFDBigh8PBISCAYd06C9q1E71dLEIqJLosI163cqUWDoc89YnNxmH1aq2X\nS0RIxUVJgXhdeLjrtCdVq1acKQUI8TXUfES87j//seHSJR4nTgho00bE6NG53i4SIRUWJQXidZUr\nAxs2WLxdDEIIqPmIEELIPSgpEEIIUVBSIIQQoqCkQAghREEdzR525QqHtWu1MBoZXn3VjsBAb5eI\nEEKcKCl4UFoah6efNiA9Xa6g7d+vwapVNOqGEOI7qPnIg44f55WEAAA//ijA4fBigQghJA+3JIVp\n06YhIiICwcHB6NKlC65cueKOzfidunUZBIG5PNZQXY2UUUoKh8OHBdy65e2SkPJE9aTw+eefY/Xq\n1di7dy9u3LiBRo0aYc6cOWpvxi81aCBh8WIrYmJEdOjgwMqVZm8XifipX38V0KZNEHr0MKB9+yAk\nJnLeLhIpJ1S/Tp09ezZmz56NevXqAQDmzp2r9ib8Wo8eDvToQW1G5P7Mm6eDySQngtRUHkuX6jBt\nms3LpSLlgao1hevXryMxMREZGRlo3LgxqlSpgl69eiE9PV3NzRBS4QUEsCIfE1JWqiaFq1evAgA2\nbtyI+Ph4nD59GlevXsXgwYPV3AwhFd7779tQq5Y8m+yjj4oYNszu5RKR8kLV5iPG5KuVsWPHolq1\nagCAyZMn45lnnkFubi50Oh14ngPP52//FARe+anR0KCowtwbJ1KwihCjBg2AEycsuHkTeOABgOM4\nACXvV6gIMbpfFTVGqiaF6tWrAwBCQkKU56KiosAYQ2pqKmrWrImwsKA7O7ArQZBX2jIa9TAag9Qs\nVrlkNOq9XQSfVxFiVKXK/X2+IsToflW0GKmaFGrWrAmj0YiEhARER0cDABITE6HVahEREQEAyMw0\nFVhTMJnkm7iysy0QRUHNYpUrgsDDaNTfiRMtRlMQilHxKEbFK28xCg0t2cW2qklBEAQMGjQI06ZN\nwz//+U8EBwdjypQpGDhwIHheroJJEoMk5e8Uuxt0UZTgcPj/F+BuFKfiUYyKRzEqXkWLkepDUmfM\nmIHc3Fy0atUKDocD//73v/HZZ5+pvRlCCCFuoHpS0Ol0mD9/PubPn6/2ryaEEOJmFatbnRBCSJEo\nKRBCCFFQUiCEEKKgpEAIIURBSYEQQoiCkgIhhBAFJQVCyonERA5nzvBgNGFquXTqFI9evfT417/0\nOHLEfaduWveLkHJgzhwdZswIAAD06GHHl19aUcAUY8RPmUxAnz56ZGbKyeDUKQHHjt1G5crqb4tq\nCoT4ObMZ+OgjnfI4Lk6LY8fo0C5PUlM5JSEAQHY2h2vX3PMd055DiJ/jefnfvWjt7/KlZk2Ghg1F\n5XHduhIefNA98zFRUrjjxg0OgwcH4oUX9IiLoyOqpI4e5fHjjwIsFm+XpOIKDASmTbNBEOTOhJde\nykWzZhVnAreKQKsFNm+2YPRoG95+24a4ODP0bprRm85+d7z2mh7Hj8tTdh86JKBuXTOaNKEDqyjT\np+swd67cjh0dLWLbNvftqKRor71mx/PP22GzcahenXqay6MHHmAYNy7X7duhmsIdZ886QyGKHM6d\no9AURZKABQuc7dgJCQL27aN1MLwpNBSUEMh9ozPfHU8+6VD+HxTE0KqVWMS7Cc8DwcGuzxmN3ikL\nIUQ91Hx0x6JFVixZIiEtjUOfPnbUqeMfV1yiCAheukD/4gsLhgzR4/Zt4I037GjThhIpIf6OksId\ngYHAiBHub69Ty08/8YiN1SMri8Prr9sxZYrN7du8fRsYMSIQx44JaN1axNy5Vly8eBsOB412IaS8\noOYjPzVkSCAyM3mIIofFi3Ueac//+OMAfPedFn//zWPbNi3mzJH7FCghEFJ+UFLwUzk5ro+zstx/\n++r1667buHqVdh9Cyhs6qv3U8OF25f+PPiqiQwdHEe9WR69edvC83NciCAy9etmL+QQhxN+4reL/\n9ttv47PPPoMk0Vh/dxg/3o4OHey4dYtDu3YigoLcv82uXUVs327GyZMCYmLEfDdIpaVx2LFDg7Aw\nhm7dHDT3DiF+yC1JISEhAd988w04Oiu4VcuWnk+4MTESYmLybzczE+ja1aA0Kb30Ui4++cT9nd+E\nEHWp3nzEGMPQoUMxZswYtX818WH792tc+hjWrdN6sTS+68IFHi++qMdzz+mxZw/d7Ed8j+pJYdGi\nRdDr9ejXr5/av5qoxGoFfvpJwMmT6n391atLeR77x30ensQY8OKLeuzZo8GhQxq88ooeSUlUmya+\nRdWkkJKSgkmTJmHhwoVq/lqiIosF6NHDgD59DOjaNQizZumK/1AJtGolYeJEK2rUkNC4sYivvqIZ\n8vIymYDkZOchZ7NxSEyksR7Et6japzBmzBgMGjQIDRo0wJUrVwp8D89z4Pn8V0eCwCs/NRo6UApz\nb5zKYt8+AQkJzmaLefN0GDdOnU7hkSNFjBx5bzLwzvd4vzFyl8qVgbZtRRw6JMe/enUJzZszj+zv\n165x+OsvDo8+KiEszHdj5EvcESOzGZg6VYcLFzg8/bSIQYPcP2qwtFRLCnv27MHBgwfx5ZdfApD7\nFgoSFhZUYAe0IMhTJBiNehiNHhhK4+eMxrJNRxoR4fo4OJhDWJh74p2VBfz1F1CvnnfmRSprjNxp\n925g/nz57vDBg3lERbl/X//5Z+DZZ+UTUtWqwP79wMMPy6/5Yox8jZoxGjcOWLJE/n98vAZRUQH4\n979V+/WqUC0prFq1CqmpqahduzYAQJIkMMZQtWpVLFiwAL179wYAZGaaCqwpmEzyFWZ2tgWiSB1w\nhREEHkaj/k6cSj/6KDoaGDxYh6VLNahUCfj8cxtu3lR/zqIzZzj07KlHejqHatUkxMVZ8fDDnuln\nuN8YudvQoc7/37zp/u1NnRoAs1k+1FNTgU8/tePjjx0+HSNf4I796MiRQADO89vBg7no1Mkz9/uE\nhpbsAkS1pDBnzhxMnTpVeZycnIy2bdvi1KlTCA0NVZ6XJAZJyn9yuBt0UZTgcNBOWpz7idPUqVZ8\n+KG8cAfHAQ431GBnzw5Eerqc/FNSeHz2mQZz53p2iCrtSzK9nuV5LNHxVgpqxqh9e4fSfMtxDG3b\nOnwu/qolhZCQEISEhCiP7XY7OI5DjRo11NoEUZFOnf7lQuWdD8lbM7kSYMIEG377TcCVKzyaNRPx\n5pu5oMkMvOP993NRtSrDxYs8unZ14MknfW9mYbfd0VynTh2Iou/9wcQzRo+24eBBAVev8qhTR8Ko\nUf4zA21hsrIAg0GuYZXGhQs8zpzh0by5iNq1PT9U98EHGY4cMSEnh9a8KI0LFzikpACPPCIvYKQG\nQQCGDPHt6WFofkviFg8+yHDokAkpKfLykO6umbiTJAFvvhmIzZu1CApi+OorCzp2LNkFz549Al56\nSQ+7nUNQEMPWrWY89pjnmws4jhJCaWzfrkFsbCAcDqBKFQO2bzehbt2Kce8N1SELsXatBv/4hwFP\nP23Ab79RmMoiIACoXdu/EwIA7NqlwebNcvXAZOIwZkxgiT+7bJkOdjunfHblSu/f6W02AwsXajB9\nOpCSQjfPFWThQi0cDjk26ekc1qzx/vfmKVRTKMC5czxGjQqEJMk7xcCBeiQkmLxcKuItljz34ZnN\nJT+Rhoa6Xl2GhXn/anPAAD3275cP/SVLAvHTTyZUquTlQvmYvLUqo9H735un0CVwAa5e5ZSEAAA3\nbnDI9f8mcVJGTz/tQHS03FzEcQxjx5Z8FNUHH9gQEyNCo2F48kkH3nrLuzvSrVtQEgIAXLnC4/Rp\nGgWQ15QpVtSrJzfzderkwKBBvt0PoCaqKRSgRQsRNWtKygRvzzzj8PsmEFJ2BgMQF2dGQoKAKlUk\n1KtX8qvGatUYdu40u7F0pRMcDISHS0hLk/dtnY6hVi3fGhJZWvv3Czh3jkf79iIeeUSdv+WhhxiO\nHLEgODgIOTk2twzb9lWUFApQuTKwc6cZmzZpYDQCffpUnKsEUrDAQKBNG/8fTScIwJo1FnzwQSBs\nNgGjR9tQq5b/No2sWqXF22/LfTyBgQxxcWZER6uX5ApbatZsBn7/nUeNGsyv41cQSgqFqFaN4c03\nKRmQ8qdpUwnbt1sRGhqEmzdFv74KXrvWeQqzWjls3apFdLSzeS8uToONGzWIjGQYP96G4OD73+at\nW0D37gacPy9Aq2VYuNCKHj38OIh5UFIghPityEjXq/SICGct4fBhAbGxgWDs7p31HJYts973Njdu\n1OL8ebkfxm7nMHOmjpICIYT4gilTbLh5k8O5czw6dXLgtdectfuEBF5JCPJjdTrU8/Yvlrf+RkoK\nbnbrFjBrVgAyMji89JId7dr5f7s0Ib4iPJxh3bqC1+5o3VqEIDCIopwY2rZV59jr08eOuDgN9u7V\noHJlhhkzyteys5QU3Oy115xjwnfs0CA+3lSq0SuEkLJp1kzCmjUWbNmiRc2aUr7hwDYbsGmTBnY7\nh3/9y17iezUCAoCNGy1IS+NgNDIEBLih8F7k90khJ0e+Jb1SJXnoKO9jd14cOeKsslqtHE6dElCv\nXvlpfyTlz+rVGvz+u4Ann3SgSxf/rtk++aRY4KRzjAH9++uxd698ClyxQosdO8ylOsGHh5fPizsf\nO4WWjskEdOtmwIgRerz2mh4jRpR8+gE1/fyzgBdf1OONNwJx9arr3a6tWjl3yMBAhsce8++DjJRv\n8+bpMGqUHkuX6jBggAG7d5fPG9v+/ptTEgIA/PabgLNn/fp0qBq/jsKRIwLOnnXutOvXa2G9/8EF\npZKYyGHgQHkx9i1btOjXz3WVpq++suD113Pxwgt2rFljcVvT0bffatCtmx4DBuhx6RLNZ0PK5uef\nXZPATz/5fWOCQhSBBQu0GD48EPv3C6hUyXksarUMVauWzyv/0vK7b/zSJQ7vvBOItDQOTz3lAMcx\nZYRBaKjn2/cuXOBhszlPwufOCbDbndMrh4YC06e7tyPq7FkegwcHKh1qV67osW+f79xFS/xHo0YS\n9u93fextN28CixfrYLNxGDQoFzVrlu3kPWOGDvPmySeI9eu1GDvWii1btMjN5fD++zZs2aLBhg1a\n1K4tYdYsGyIjS/67r1/nkJTEo3FjUZV7IbzJ75LC4MF6Za6W8+cFxMbmYssWuU9hzhyrKgvQl8Zj\nj0moXJnh1q27IxwcpZ5v/3799RevJAQAuHiRhyTB5/pXfF16OodLlzjUry+hcmVvl8Y7JkyQL2B+\n/53Hk0+KGDiw7Ddw3r4NHD8uoEYNhvr1y5ZcJAn4978N+O03+ZjfulWDvXtNZTrxHjjgerrLyeGx\nf7988fT99wImT5Zr+WfOCLBaOWzeXLKLufh4Aa+8oofVyqFmTQnbt5tRo4b/1jr8LikkJrqe6Zo0\nETFtmveGhFWvzrBtmxkrVmhhNDIMH166Cc927RLw668axMSI6NatbB3QLVuKCAuTkJkpx6ZTJ5ES\nQimdPMmjVy8DsrM5VKkiIS7OXCFHiQUGAlOn3v/xdOsW8OyzBly8KIDnGWbNspUpwaSlcUpCAIBr\n13icP8+jRYvSJ5nmzUUcPy64PL7rr79cD5hLl0p+AM2erYPVKl+UXb3K4+uvtRg3zn9n0PS7pNCz\npx3ffCPfLRIWJuHxx73fcduwoVSmscpbtmjwxhvOPojPPrPgxRdLnxiqVWPYvt2MtWu1CAlheP11\nmp6jtObN0yE7++78+TwWL9Zh1qzyNf7ck+LitLh4UT4BSxKHOXN0ZUoKoaEM1atL+Ptv+SQdFMTK\nvHrdhx/aEBTEcO4cj86dRZe7kDt2FDFzJlOmRX/22ZIfh3mbrPX6gt/nL1RPCklJSRg1ahT27t0L\nrVaLp556Cp999hmMKi37NGuWDS1bikhN5dGjhx0REf57Nbd7t2v4v/9eU6akAMizOr7/vv9enXhb\n3gO7vI099zSDwfW4DAoq23Gq0wHr1lkwbVoArFZgzJjcMncIBwQA48cXfIw0aCA3++zYoUGtWhL6\n9HGgpONwJk+2oW9fHqmpPFq0EDFokH8fh6onhe7du6Nly5ZITk7GzZs30bNnT7zzzjtYsmRJkZ/L\nyZF/MgasW6fB3LkBCApi+PhjK5o3d1YVeR7o29e3x/mvX6/B4cMCYmIk9O9f+NVR3ml+GzTIXyW+\ncYODwcAQEqJ6Mck9xo614dgxAUlJPB55RPT6ugf+7vnnHdi1y464OC3CwqT7qnU1bChh5cqC71pW\nU+PGEho3Lv33/uijEhISTMjK4hAWxjzer6k2VZNCVlYWWrZsiRkzZkCv10Ov1+Pll1/G/Pnzi/zc\n9u0aDBliAAD07RuAo0edq54NGKDHH3+YVAv0339zCAxkbutIXLtWgxEj5PrjypXyql2FNecMH56L\nW7c4/PqrgJgYEaNHO3dIxoChQ+V1gbVahjlzrOjd27eToT+rW5fh8GETMjI4hIcz6pMpgevXObz8\nsh5nzvDo0EHEkiUWGOTDGIIALF1qhclkhV5f/gc9aDTAAw+UrQZz65bcoR4WpnKhykjVryokJARL\nly5FeHi48lxSUhIiixnbNX58gDKs89dfNS6rnqWn8zCrNLpyxIhANG1aCY0aVcKqVe4ZInTvDTEA\nsG9f4Tf/aDTApEk27NxpxtSpNpeJtfbuFZR1ge12Du+9Fwjmvy1lfkGjkftnyvsJTC0TJwbg1CkB\ndjuH77/XYPHi/DPDBQWV/4RwPxYs0KJBg0p45JFgfPRR/vjFxWnw6ac6nDrluSC6dUvHjh3DggUL\nMGHChCLfl3c+9/BwZzPKU0/ZERR0/2U5fFjA2rXaO9vjMG5cAEQ39FHnvWP5scfKNhQvb9kkCZQU\n7pPDIdfkli7VIj3dN+r4JpM8pPH33/3vzJmZyRX5mBQtLY3DlCkByn1Ws2cHIDHRGcN583R4/XU9\nZs4MwLPPGnDihGf2EbeNPjpw4AB69OiBjz/+GB06dFCe53kOPO+680yebMdbbzGIojwtxFdfWRAX\nJ9970LevAxrN/Qfj3toHIJ90eZ4vdGWlsho6VERubi4OHuQREyNhzBgHBKH05e/YkaFTJwf27NGA\n4xgmTsyFTscrv6ssv7OiKCxGr78egLi4uwvWS/j5Z0u+Bdo9KTsbePppPc6dk8s5fboNQ4Z4pomw\nNPuR1Son1LwTxg0e7MChQwIcDg7BwQwDBqhzrPoK9x9rnMvU3gAgSTw0Gvnq79tvnSen3FwOP/yg\nRatWHhhZyNwgLi6OhYSEsJUrV+Z7TZKkAj9z9mwWA8AyMrLcUSQmioz16MEYwBjHMfbJJ27ZjKpE\nkbGEBMZw2IbwAAAV30lEQVQSE71dEv9ntcrf/b3/du70bpm+/tq1POHh3i1PQZYtY0yrlcs3blz+\n13/7jbG1axm7csXzZSsPRo50fv+vvur6Wp8+rvvH8uWurycmMrZ4MWM//KBumTjG1G2UOHjwILp3\n747169ejU6dO+V7PyLidr6YAACbTbdSqVR3JyX8jKKiEc9iWEmPAmTMcKlUC6tTxz7YYQeBhNOqR\nnW2BKHp/CgJfVFiMGjbUIyVFvurjeYZ9+yxo2NB7+8F33wl46SXnJI5RURJOnHD/KBugZPuRxQJE\nRRlgtzuP1337zGjc2D+PndLy1LF25gwHSeLw6KOu28jMBN5+OwAXL/J49lkHxo+3KwNu/vqLQ+fO\nemRlyU/89782DB9edC0zNLRk7fCqNp6IoojY2FjMnDmzwIQAAJLEIEn5d6q7QRdFCQ6H+76ABg3k\nn/68Li1QfJz+/FOeprtpUwkPP1x0PFes0OLIEQGtWol46SX3VE8lSe40u3WLQ/fujjKP1MgrNZXD\nrFk65OTI8+K0bOn8W/PGaMUKC959NxA5ORxGjrTh4Ye9uz5x164S+vQRsG6dFqGh8gize8t78ybw\n+ec6mM0cXnst1y13WBe1H1mtcEkIAJCTw8p0fN64wWHQIHmkUseODnzxhRWBRUxqHB8vz1r6+OMi\nmjTx7sWPu89J9evLP/Pui0ajPKGmsxzO17Zt0ykJAQBWrtRgyBB1hlGrWlPYv38/nnjiCQQEBIAx\nBo7jlJ/nz59HrVq1kJaWU+BnzebbiIqKwOXL12EwuKemUB5oNPydBddNhe6oBw8K6NNHD5uNQ0AA\nw5o1FrRvX3Cv+rJlWowb5zw6Z8604tVX1U8Mb70ViHXr5I7+OnUk/PijqcT3XthswIkTAsLCWL57\nOTp3NihzYQUFMezbZ0JUFFdsjHyJ1SrfWJV32HXXrgacPCn/bVWqSPjpJzPS0zk88ABD9er3d9iW\nZD8CgEmTAvDFF/KomC5dHFixwlKm0USxsYHYts054m/CBBtGjCj4JLZ8uRbvvSfvkzodw5YtZpdk\n7ykljZE3rF+vwfDhzlun//lPBzZtKrqWGR5esgmjVO1Bad++PURRhNlshsVicflZq1YtNTdFirBi\nhVYZ4muzcfjmm8KH3x486Dpk9tAh9efPdziADRucldIrV3gcPlyy7VgsQM+eBjz3nAGPP27A4sXO\nvyU3F0pCAACTicOZM/7X0RkYmD8hZGdDSQiAPDS7Vy89OnQIQkxMEDZt8swMNZMm2bBnjwnffWcq\nc0IAkG+0V1pa4SOVNm927WD99lsPzzDpB3r1cmDQoFxUqSIhJkbE7NnqrRngf0eQl1gs8uyjJpO3\nS1K8KlVcryKLaqpp1kws8rEa7o7/v1dJZ5H84QeNMokZYxxmznTOP6HTAdHRzvIGBTE0buxbV3Rl\nFRws9zHcpdMxnDsnx8Fu5zB5sufm4WjSREKrVtJ93W/wyit28Lz8nQcFMfTpU3httFYt132jdu3y\n8Z2qieOAGTNsOHPGhJ07zahTh2H/fgHLl2vvez0Vv5sQzxuSkjg8/7wByck8qlWTsGmTpcxTAXvC\nO+/YcPYsj6NHBbRsKeLddwufYmDoUDtEkcPRo3KfwhtvuKdP4X//s2D06EBkZXEYNiwXTZuWLH56\nvesJIu+cOqtWWfDpp/Jkdq+/novISAbA/8fLcxywfr0ZU6cGwGzm0LSpiNmznYlA7aHU7vbccw7U\nrWvGuXM8WrUSERVV+EXBlClW5OQAZ88K6NTJ4ZbmzPLmq6+0+M9/5Ca3SpUYduww55tGp6RUH31U\nHH/sU3j33QB8/bXzbsMXXrBj0SIPL/F2hy+3c7oDY3J/xPr1WhgMDF9+acH//V/htZkLF3hs3qxF\n7do69O1rAs/7T4wcDvmGpVOn5A7WQYOcJ0O7XZ7y5aefNAgMZFi0yIpnnil7L3l53Y/sdnlAg90O\ndO/uuK8bX/0pRvf2rQHAmDE2jB3r2mdT0j4FP7ve8A10Z7HncBywYIEVM2bIo1WKWsDo6lUOzz5r\nUEZl7N0bgEWL1B/iabEAH34oT/HQvr0D77+fC+GeLpL//U+Lr7/WompVhlmzrCUe/vzppzp8+qlc\nG9i5U06Cd2fN1WqBtWstSE7mULky89pNd5cvc9i1S4MaNRiee873hvC98ooeP/wgn9aWLRPx3Xdm\nl+ljyqsaNRhOn3Z9XFblLils367BgQMCoqNF1SaQGzYsF3v2aHD1Ko+qVSWXieuIZ5Rkpa1DhwSX\nYXq7dpW90/yPP3iMHx+AnBwOI0bkomdP57700UfOmmNCgoDwcIahQ+Wr+sOHBYwdK1fjz5yRVwrc\nvbtkk3fduwAMABw7JrhMpc5xKPNaAkU5cYLHpk1aVK/O8MYbuYWeRJOSOHTtGoSbN7k7n8vF5Mm+\ns+ZESgqnJARA/m7++INHs2a+fZWvhpkzrcjODsTFizyeecaBAQPK3uRWrpLCtm0axMY6h2nl5Fhd\nquBlFRXFsH+/CcnJPCIjpXy3+5dn8fECvvpKh9BQhg8+sOXrMPYlDz8suazZXdBU5CXVv78e16/L\nPatvvhmIRx81KfcJXLjg2uN67+PLl137M0qzglfbtiJ+/tl5SLZp4/4FpM6d49Czp15ZOezcOR6f\nf15w0+iePRolIQDApk0an0oKwcEMQUEMJpNcRo2GITzcd/dXNUVEMGzbpk6tuFwlhfh4Tb7HaiQF\nADAY7u8k422pqRxu35aniC7pNOQXLvAYOFCv3MD05588du1SacpaN4iOljB/vhXLl+sQESFgypSy\nnbAsFigJAZAnUExK4lGvnnySfuopeU4qAOA4hi5dnCfv9u1FlzW7u3Ur+f43cmQu9HqG06cFPP64\nA716ub955tAhQUkIAPDTT4XXrurUcd3/844Sul8mE+5rmm2DQb7Za+zYQNjtwPjxNtSsWTGSgpp8\nPimYzfK6AgcOaPDYYyKWLLEWOsSycWMRgLPRuU4dEUOGBCItjcPLL9tdlt8DgIQEHkeOyHf9euKq\nzFtWr9ZgzJhAiCKHp56y43//s7q0gRfm7Fne5Y7W06c9M4J53z4Bb70lj1QaOTIXo0aVvLmud28H\n+vWT7nQQsjLdsazXA507O/Djj/LhERkpuazn+/LLdoSFMZw6xaNdOxEdO4o4dozH7dsc2rUTsWuX\nCVu3yn0KL75Y8qTA88CQIXYAnhtt07Cha+2qqCG9HTuKGD/ehtWrtahRQ8LcueoMtnA4gMGDA/Hd\nd1o88ICE5cutaN26bMdjx44ijh71g3HjPsznRx/NnOnsfAOAAQNyMXt2wVeAkiR31h04IOCxxyQc\nO8bj6FH5wBYEht27zcpQyF9+EfDii3o4HBw4juHLL635koYvKsuIiLp1KylVagBYvdqMzp2LP+iS\nkjg88USQ8tknnnBgwwb3z83zyCNByMx0JqDdu02lahdWY9SIzQasXKlFTg6Hvn3tRd5BPHFiABYu\nlBvi27WTY1RUh7gvuDdGa9bwWLtW7lOYPNmW7z4Xd9uwQYNhw5zNvg0bivjlF+/XSNUafXTpktzX\nUbs2w9NPe+8cU25GH6WmcnkeF361yvPAu+/m4t135cdRUc7kIooczp7llaSwebMWDof8uxnjsHGj\nxueTgt0urxznqT6N2rUZtm0z45tvtAgLYx5ZotJuh9L0cpc31j4ICECJmh6tVmDRImcGOHhQXor1\nn//0n5pnr16eaaoqzL0XLAU9zmvixABs2KBBZCTDwoUWt8wJpZa//uLQpUsQcnLkv+m992x45x3f\nHqji83c09+1rR2Cg/KULAityzeO8OnRw7ugGA0OrVs4DNTLSNfPLNz35ruRkDv/4RxCaNjWgQQPg\nypWSnyinTbNCEOS/76mn7OjQoeQnrKZN5fV1//OfXI8kI61Wbp65q2FDEe3a+e4JVqORm5vuFRzs\n2/uSr+nZ06701/A8w9tvF37S/PZbDRYu1CE9ncepUwLeektf6Ht9wa5dGiUhAMDGjT5ehYQf1BRa\ntpSwZ48Zx47xePRRqVQzJi5caMWiRRLS0zn07m1H3brOg3XEiFxcvszj0CF5+Or48aXrlBRFee6d\nvCcEd5k7V4fLl+UcnpgIfPyxFp99VrKTZb9+DnTubCp1R7O3zJxpQ9euDmRnc+jc+f5uQCrK7t0C\nVqzQoUoVhgkTbGUaqaLRAPPnW/HWW4GwWuX9Kjrafwck3A+7HdiyRQObjUPPnvYihxFbLMCaNVrY\nbECfPnZ8/70ZJ04IqF6dFTlbQEoKV+RjX5O3ozvvxagv8vk+BV8UHy8gNlaPnBwO/fvLfRzuPtGO\nHBmINWucVxm9ejnw+eeemXvfHZKTOWzZosUDDzD07WsvUcd3SZWkLfj333l06WJQmhDbtnXc15A+\nSZI7TP3lRil33K07YIAe338vX2c2aiRi505zgRdNjAH/+pce+/fL7334YRE//GCGwVD8Nq5e5dCl\niwHp6fIF0tixNowZ457mmLLG6JdfBHz4YQDsduD993ORkHD3LnsJ8+ZZVR+1VVIl7VOgpFAGTZsG\n4e+/nS1va9aY0amTe5s4Llzg8fzzeqSl8QgPB7ZsMaN+fd9tVilKSgqHDh2cB3avXvZCx8aXRUkO\n5nXrNC5NDwYDw+XLt1Urg69TOylkZHBo2ND1uP32W3OBo4jS0jg0buz63u++M6FVq5KV4/p1Dnv2\naBAZKaFjR/cdA2WJ0e3bQNOmlXD7tnyxERDAcOyYySfu7/HK1NkVxd0vvLDH7lC/voRDh0zYs8eC\nCxeARo28v5OV1cGDgpIQACjrJntSixaiy2R7ha03QUqmUiXm0pciCAzVqhV8Ig0JYQgNdb5Xp2Ol\n6tOLiGAYONDu1oRQVjdvci7nA5uNyzdYxtdRUiiDMWOc/Q9Nm4ro3NkzIzeMRqBZMwmVK3tkc24T\nFSWPjb/3sac99BDD5s1mDByYi9GjbVi82H+b4nxBQACwfLkF9eqJqFlTvoehsJlQdTrgm2/MiI4W\n0bChiMWLrT4/0OP6dQ7jxgVgzJiAIqemjoxkaNfOeT5o0kT06RmVC0LNR2X0++88MjI4tGoleqyz\nGfCvmRuLsmKFFkuXykNdZ82yFbtkaGmUlxi5E8WoeHdjlJpqQtu2emXKkmrVJBw8aCq0I91iATZs\n0MLhkJtGSzJvlyeUm/sUfFXeRbZJ6bz0kt1t60EToqaUFM5lDquUFB6XLvF47LGCzwF6Pfx636bm\nI0IIKULVqsxlKGlYmJRvHqjyhGoKhBBSBJ0O2LTJjFmzAuBwyBMX+nu/XlEoKRBCSDEefJBh4ULv\nrLboaR7vaC5MdnY2QkJCkJWVBaO3lpUihJAKzmeSAmMMOTk5CA4OBufr8zAQQkg55TNJgRBCiPfR\n6CNCCCEKSgqEEEIUlBQIIYQoKCkQQghRUFIghBCioKRACCFEQUmBEEKIgpICIYQQBSUFQgghCkoK\nhBBCFB6dJfXu/EaEEEI8ryRzy3k0KeTk5CAkJMSTmySEEHJHSWah9uiEeFRTIIQQ7ylJTYFmSSWE\nEKKgjmZCCCEKSgqEEEIUlBQIIYQoKCkQQghRUFIghBCioKRACCFEQUmBEEKIgpICIYQQhU8khaSk\nJHTr1g1VqlRB3bp1MW7cOG8XySckJSXhhRdeQJUqVVCjRg28+uqryM7OBgDEx8ejdevWCAkJQZMm\nTbB69Wovl9a73n77bfC8c3em+LiaNm0aIiIiEBwcjC5duuDKlSsAKE53JSQkoFOnTggNDUVERAQG\nDhyIjIwMABUwRswHxMTEsCFDhrCcnBz2559/svr167M5c+Z4u1he17RpUzZo0CBmNpvZtWvXWMuW\nLVlsbCy7ceMGq1SpElu+fDmz2Wzsxx9/ZAaDgR0/ftzbRfaKkydPsgceeIDxPM8YY+z69esUn3ss\nWLCANWrUiF28eJHl5OSwkSNHspEjR9J+dIfD4WARERFswoQJzG63s8zMTNalSxfWu3fvChkjryeF\no0ePMq1Wy7KyspTnFi1axBo2bOjFUnnfrVu32KBBg1hqaqry3IIFC1iDBg3YJ598wmJiYlze37dv\nXzZ06FBPF9PrJElibdq0YdOnT1eSwqxZsyg+93jwwQfZ1q1b8z1P+5EsOTmZcRzHzp07pzy3aNEi\n9vDDD1fIGHm9+ejEiROIiopymbmvefPmOH/+PEwmkxdL5l0hISFYunQpwsPDleeSk5MRGRmJ48eP\no3nz5i7vb968OY4ePerpYnrdokWLoNfr0a9fP+W5EydOUHzuuH79OhITE5GRkYHGjRujSpUq6N27\nN9LT02k/uiMyMhLNmjXDkiVLYDKZkJqaik2bNqFbt24VMkZeTwoZGRkIDQ11eS4sLAwAkJ6e7o0i\n+aRjx45hwYIFeP/99wuNWUWLV0pKCiZNmoSFCxe6PE/xcbp69SoAYOPGjYiPj8fp06eRnJyM2NhY\nitMdHMdh48aN2Lp1K4xGI2rUqAFRFDF9+vQKGSOvJwVAnlKbFO7AgQPo2rUrZs6ciY4dOwKgmAHA\nmDFjMGjQIDRo0CDfaxQf2d04jB07FtWqVUNERAQmT56MuLg4cBxHcQKQm5uL7t27o0+fPsjKysK1\na9cQEhKC/v37A6h4+5JHF9kpSHh4uNLLf1dGRgY4jnNpOqmovv32WwwcOBCff/65spMWFrOqVat6\no4hesWfPHhw8eBBffvklANcDl+LjVL16dQBwWdwqKioKjDHY7XaKE+R96fLly5g+fToAoFKlSpg0\naRKio6Px9NNPV7gYeb2m0KJFCyQlJSEzM1N57siRI2jUqBEMBoMXS+Z9Bw8exCuvvIJNmzYpCQGQ\nY3b8+HGX9x49ehStW7f2dBG9ZtWqVUhNTUXt2rURHh6OmJgYMMZQtWpVNGnSBMeOHXN5f0WLz101\na9aE0WhEQkKC8lxiYiJ0Oh2eeeYZihMAURQhSRIkSVKes1qt4DgOnTt3rngx8loX9z3atm3LYmNj\nWXZ2Njt79ix78MEH2cKFC71dLK9yOBysUaNG7Msvv8z3WmpqKgsJCWFfffUVs1qtbPv27SwoKIj9\n/vvvXiipd9y6dYtdu3ZN+Xf48GHGcRy7fv06S0pKqvDxudfo0aNZvXr12J9//slSUlLYP/7xD/b6\n66/TfnRHRkYGCw8PZxMmTGBms5mlp6ez5557jnXo0IGlpaVVuBj5RFK4du0ae+aZZ5jBYGA1atRg\n//3vf71dJK/bt28f43me6fV6FhgY6PIzKSmJ7du3j0VHR7PAwED2yCOPFDjksCK5fPmyMiSVMUbx\nuYfNZmPDhw9nYWFhzGg0stdee42ZTCbGGMXprhMnTrAOHTqwsLAwVqNGDfbiiy+yGzduMMYqXoxo\nOU5CCCEKr/cpEEII8R2UFAghhCgoKRBCCFFQUiCEEKKgpEAIIURBSYEQQoiCkgIhhBAFJQVCCCEK\nSgqEEEIUlBQIIYQoKCkQQghRUFIghBCi+H8kppMjxHOpqgAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list_plot(a1, figsize=4, zorder=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

異常度の高い車種

\n", "\t

\n", "\t\t異常度a1をXに追加して、異常度の高い順にソートしてトップ5を取ってみると、\n", "\t\t井出本の通りの結果となりました。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 異常度をXに追加\n", "X['a1'] = a1" ] }, { "cell_type": "code", "execution_count": 14, "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", "
a1PriceMPG.cityHorsepowerLength
Make
Chevrolet Corvette13.59583038.017300179
Honda Civic11.82974212.142102173
Geo Metro11.1563678.44655151
Mercedes-Benz 300E10.58602561.919217187
Volkswagen Eurovan9.97114819.717109187
\n", "
" ], "text/plain": [ " a1 Price MPG.city Horsepower Length\n", "Make \n", "Chevrolet Corvette 13.595830 38.0 17 300 179\n", "Honda Civic 11.829742 12.1 42 102 173\n", "Geo Metro 11.156367 8.4 46 55 151\n", "Mercedes-Benz 300E 10.586025 61.9 19 217 187\n", "Volkswagen Eurovan 9.971148 19.7 17 109 187" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 異常度の高い順に6個出力\n", "cols = ['a1', 'Price', 'MPG.city', 'Horsepower', 'Length']\n", "a5 = X.sort_values(by=['a1'], ascending=False)[cols]\n", "a5.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

カーネル主成分分析

\n", "\t

\n", "\t\t井出本では、カーネル主成分分析の式の導出をしていますが、潜在変数の導入理由や考え方は、\n", "\t\tPRMLの説明の方が分かりやすいのでここでは省略します。\n", "\t

\n", "\t

\n", "\t\t井出本ではRのkernlabパッケージを使っていますが、ここではsklearnのKernelPCAを使って\n", "\t\t計算してみます。\n", "\t

\n", "\t

\n", "\t\t最初にsklearnのPCAを使って第1主成分と第2主成分にXcをプロットしてみます。\n", "\t\t結果は、先ほど求めた結果と同じになりました。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3WlgE9XawPH/zGRtS0srUIpsCgiIonBlrQp6kYK7RRZZ\nioqKIAKC4hVBQNxAQd6r6AWvG4si+yICegEFWRVxQRQXFLAoa+mSpk2TzPthJKVsTdO0Sdrn96Wd\nNpmck5k5z5xlzlF0XdcRQgghADXUCRBCCBE+JCgIIYTwkaAghBDCR4KCEEIIHwkKQgghfCQoCCGE\n8JGgIIQQwkeCghDlSNd1srKykMeDRLgyBfrGI0eyUVWFhIRojh934PVG9kkueQlfFSk/Dkc2F110\nIb/9lk50dJVQJ6dUKtJxqSx5qV69+HOuVDUFVVVQFAVVVUqzm7AgeQlfFSk/iqIU+RnJKtJxkbyc\n8v4gp0cIISqtr79W2bhRw+UKdUoCJ0FBCCFKqKDgzL89/7yFzp2j6dYtijvvtEdsYJCgIIQoNY8H\nPv1UY9MmLdRJKVMOB9x5p50LL6xCcnIUv/1mNNG43fDvf1t8r9u61cTnn0fmdyFBQQhRKl4v9O1r\np0ePKO64I4qhQ22hTlKZmTnTwoYNxvicn3/WGD/eCoCqgu20bFepEpkd1hIUhBClsmuXytq1hQMZ\n580zc+hQ5HfYnk1mZtF8ZWUZ26oKr7ySR1SUjqLoPPigi1atvOWWrvnzTTRsGMPFF8cwa1bAg0oB\nCQpCiFKKiTEKwpPMZp2oqMi8Sy5Ov34uEhKMwt5sNgr/k26+2c2vv+awb18OTz+dX25pOn4cHnnE\nRlaWQk6OwsiRFg4eDHx/EhSEEKVy8cU648blYzbr2O06U6fmUSWyH8E4pwYNdD77LJfZs3PZsMFB\nSoqnyP817cxmpLKWna1QUFBYg/F4FE6cCHx/EhSEEKU2eHAB+/fn8NtvOfTs6Q51cspUYqJOSoqH\nBg3CozZUt65O166Fw6E6dXLTpEng+ytd45MQQvxNi8zBNhFPUeCtt/JYt64ArxdSUnRUNfCiXYKC\nEEJEOE2DG24wmrJMptI1AEnzkRBCBCg7G775Ri1VG364kaAghBAB+PVXhfbto7nhhmjato1m166K\nUZxWjFwIIUQ5e+01C4cOGUXo8eNqkSeaI5kEBSGECIDJdP7tSCVBQQghAjBsmIuLLzYeZKtd28uj\nj5bfA2tlqYLENiGEKF+1auls3Ojg0CGFGjV0LBWj9UiCghBCBMpshtq1w+MhtmCR5iMhSuiPPw5w\n9919aNKkPpdd1oihQweRnZ0V6mQJERQSFIQoob59exIfH8/OnT/wv/9tYM+eHxg/fkyokyVEUEhQ\nEKIEsrIyadGiJU8+OR673U7Nmkn06NGbLVs2hTppQgSF9CkIUQKxsXG8/PKrRf6Wnv4HSUm1QpQi\nIYJLgoIQpfD111/x1lszmTNnfqiTIkRQBBwUVFVB04zWp5M/I5nkJXyFa362bt1Cnz49GT9+Ih07\ndvTrPafmpbQTl4VauB6XQEheCim6rgc0nkrXdRSlYi65J0RxVqxYQb9+/Zg+fTp9+vTx+31ZWVnE\nxcWRmZlJbGxsGabw3I4ehV9+gcaNIT4+JEkQYSzgmsLx4w7MZo3YWDtZWU48nvJbj7QsaJoqeQlT\n4Zafbdu20r//3bzzzhw6dLiOjAyH3+91OJwAf+el/Bcg+PZbldtvt3HihEL16jorVji55JLAxtmH\n23EpjcqSl/j46GLfH3BQ8Hp13wd6PF7c7sj+Ik+SvISvcMiPx+Nh+PAhjB07geTkDiVOT6ivmZdf\ntnDihFHDP3JE4dVXTUydWrrpGcLhuASL5EWGpApRIl98sZ2ff/6JJ58cRd26NahXL9H3Mz39j1An\nr1inT8VgNocmHSJ8yegjIUqgbdt2/PVX5K6o8thj+WzbprF/v0qjRh6GD3eFLC2LFpmYPt1CbKzO\n88/n07RpxbhDj3QSFISoROrX19m61cGxYwrVqukhW1d5zx6VIUNseDxGU1a/fipfful/34woO9J8\nJEQlYzJBYmLoAgLAvn2KLyAAHDigUFAQuvSIQhIURIW0c6dKSkoUV18dxcKFUiH219dfq7RvH0Xj\nxjFMnlx2c0FfdZWHpKTC5qKUFLf0b4QJuVpEhaPr0LevnSNHjHuehx+2ceWVDho2rFhTHJeFBx6w\n8/vvxvf20ktW2rf3cPXVnqB/TkICrFqVy/z5ZqpU0enXT6oJ4UKCgqhwcnPxBQQAj0chPV2lYcPg\nF24VzeHDynm3g6lWLT2kHd3i7KT5SFQ40dHQubPbt12vnpcWLSQg+CMtrfCOvW5dLx07us/z6tD4\n7TeFFStM7N0rMyqUBakpiArprbeczJtnxuGAO+90E6IZJSLO00/nc+21bo4eVbjhBg8JCaFOUVHb\ntmn06GHH6VSw2XQ++MBJu3YS8INJgoKokCyWone9wn+dOoVvIfvuu2acTqOGkJen8PbbZgkKQSbN\nR6JcHDyo0L+/ja5do3j/fbkXEYGJj9fPuy1KT65OUS7uu8/Ol18aA+O/+spGw4a5tGolT7AKg9sN\n+flGf9BJXi/8738aLpfCDTe4sVrh0Ufz2bVLZft2jX/8w8OoUdJRHWxSUxDl4qefCk81XVf4+efK\neep98YWR748/DuGTY2Fm6VKN+vVjuOiiKowYYfX9/YEHbPTtG8W999rp2dOO221M9b1smZODB3P4\n8EMnF1wgNYVgq5xXpih3KSmFo1hiYnTat6987cCffqrRo4cNMAq8t9+Wp7VcLhg0yI7LZfQTzJlj\nYcMGjUOHFJYvL/x+Nm82sXt3YXElS7mUHQkKolxMm5bHs8/mMWxYPitX5lK/fuW7w1u1yoSuF5Zm\nH35YeVtvdR2efNJKkyYxRaa7AKP/KTpax2YrPEdUVScurvKdM6FQec9KUa7MZrj//so9GqhBA+95\ntyuTZctMvPHGmdNoqKpOSoqbmBiYPj2PUaOsuFwKY8fmU6+eBIXyIEFBiLP46SeV3btVWrTwBK0w\nGjCggN9/L+C//4XU1AKeeqrytoEcOlQ076qqU726zksv5fmWCL3lFje33BJ+D89VdNJ8JMRp1q7V\nuO66KB54wE7HjtF8801wLhNNgzFjjNEyU6e6iIkJym4j0k03ualWrbCmNHq0i+++c5CSUvn6msKN\nBAUhTvPWWxYKCow7WYdDYc6coh3C69b9j2bNGvLgg/eGInkVQu3aOmvX5jJtmpOFC3MZOlSGloYL\naT4S4jSnPxCVkFC4/eqr/8f778+mQYOG5Z2sCicpSad3b2keCjdSUxDiNGPH5nPppR7ACAZbtmg4\n/l4UzG63sWbNeurXvyh0CQyxP/5Q6N3bzrXX2njllVCnRgSbBAURFnJyYOBAG23aRDNypBVXCFsT\nEhN1qlXTAaMJaetWk++ZggEDBhITUyV0iQsBlwsyMwu3Bw2y8b//mdi1S2PoUNiwQSUvD9LTFTzS\nJRDxJCiIsPDss1aWLDHz228qs2dbeP3186/69fXXKvfdZ+Ohh2zs2xf8UTy5uUX36XAE5zP0v1ui\n0tKsPP+8JeyXoFy3TqNJkxgaNarCwIE2vF7Yu7dosfHZZxotW0bTokUMnTtHceJEiBIrgiLgPgVV\nVdA04+Q4+TOSSV5Ca//+omk9cEDFZCqaj5M/jx6F7t2jyMw0CuovvtD44gtnUNccHjGigLvvVnG5\nFJKSvPTv7/GlB0BRFBRFKfI3f8yda9Q4NmwwsWGDFU1TGD06fCPDyJE2cnKM73nJEjOpqR5uucXD\nW28Z+a5SxWheO3rU2P7uO4233rIyalT45ulsIvGaOZfS5iXgoJCQEI3y97PmsbH2QHcTdiQvoZGW\nBp98YvyuadCnj5n4+KKjfk7m5/vvizZn/P67itcbTbVqwUvPXXdB27awdy+0bKkSHx9V5P9Wqwnw\nEB8fffYdnMPnnxdtX9m920J8fNmthVxaeXlFt1XVxhtvQPv2kJ4O3bvD4MFFo7GmhXeezieSrpni\nBJqXgIPC8eMOzGaN2Fg7WVlOPJ7IfjpT09QKkZd580y8/rqZCy5QmTw5L2KWoLzpJliyROWbbzTa\ntfPQqpWXjAzjf6cfm1q1oFq1KI4eNW5KGjf2omlO3+uDpWpVaNnS+P30fefnu3G53GRkODhxArKy\nFOrU0c87J8+HH2qsXl10ZFPbtvlkZARvBM6uXSqLF2skJencc48bUynHF/7rXyZGjbKg6wrNmnlI\nSHBx+LCX1NTC4/Loo/l8+aWFnByFevW89O6dR0ZGZD19XFGufzh/Xvy5iQn4lPF6dd8Hejxe3O7I\n/iJPiuS87NqlMmSIBa/XKJl697awebMjxKnyX3Kyl+Rko4B0n6WcPHlsqlSBpUtz+c9/zNhsMGyY\nC69Xx1uOh03XdXRdZ+FClSFDbLhcCl26FPD223nnbMaaPdsKFPagt2nj5sEHXWfNayB+/VWhSxeb\nrz/k228VXn45v1T77N/fRXJyAatXm5g0yUrXrnaaNvWwfHkuF1xgvKZ9ezfbthXwxx8KjRp5iYk5\n+/GLBJF8/Z8u0LzIcwoVyL59qi8ggLGWrdcLaoQ2kx49qvDWW2Y0TeGxx4rm45JLvEydWroCLxB1\n6jRA16/F4/kDRdnB4sU5gFFNX73azJo1bm688ewlYs2aRS/Qa64Jbi1u40ZTkQ7yNWtMQOm/o4YN\ndRYtMpOfb+z7hx803n/fzJAhhemvXt2YpkJEvggtLsTZtG7toUaNwoLnpps8ERsQ8vPhttvsvPSS\nlUmTLHTsSEiHqQKcOAF16vyJy7UQj2crY8ZkU7Vq0Y6M8w3JfPJJY/1jgK5d3Tz8cHAzdMklRYNO\n48bBu+M9vfZT2mYpMIYh79ypcvx46fclgidCiwxxNtWr66xalcuYMS5eew3eeKP876SDZe9elZ9/\nLiyJdu2C/ftDO4Hc6tUmfvmlME2vvmphwoQ8NM24Q77mGjddupy73SQhAWbNMo7J66/nYw9yn2b7\n9h6mTMnjH//wcPPNBbz+el7xb/LTuHH5xMQY+WzVysNdd5VudNG+fQpXXx1NSko0bdrEsHOnFEXh\nQpqPKpg6dXRGjCggPt5CRkbktu3WquUlNlYnK8sIBPHxxkNloVTltGfW4uLgrrvcdOzo4MQJhUsu\n8QZ1WGwg+vUroF+/4A4Hzc2F556zkpOjYLHoDBrkKrJsZiBmzrRw8KARCDIzFaZOtTJ7tjMIqRWl\nJeFZhKW4OJg710lysptrr/Xw0UdnFsrl7cYb3fTt60JVdWrU8DJtmnEnnpSk07Rp6ANCWVm82MyO\nHUbmXC6FiROtxbyjeKc3P5lM0h8RLqSmIMJWmzYelixxYjKpxMdHk5FhtOtPmWLl2DGFtLQC2rYt\nvyG3igJTp+YzeXJ+UNrUI8Xpw2yDsRTmQw+5WLdOY88ejaQkL088IbOkhgupKYiwkp6ucOONUTRq\nFMOQIbYzmr/uvdfOjBkWFi4006OHnb17y7+foTIFBDAWBGrXzjgQdrvOxIml76uoUUNn/fpcdu7M\nYft2xxmd5CJ0KtnpLcLdE09Y+fJLo6li/nwzLVp4GDiwsDawbVthG01ensI332hcfHGEdpxECLsd\nFi92sn+/QkKCTlxccPZrMsGFF0qzUbiRmoIIK4cPFz0ljxwpWhO46qrCAGG16lx+eWQ8sR3pNA0u\nuih4AUGELwkKIqz061eAohh3j1Wq6Nx+e9FawNtvO7nnHhe3317A++87adiw+DvNF16wkJwcRb9+\ndg4frrzrIgvhD2k+EmGlT58CGjXy8MsvKu3be6hfX+fUe5eEBJg0yf/nLxYvNjF1qjFa5uefQVGs\nzJoVvPH7QlQ0EhRE2Gnd2kvr1sHpePztt6KV4X37pHIcCVauNJGertC5s/vvGwNRXuQKERVaSoob\nm62wULn1VumUDnfPPWfhnnvsjBljo3Pn6DJZREmcmwSFCDFpkoXLLoumU6cofvxRDpu/LrvMy6pV\nuTz5ZD5vvulk5EgZDx8M+/YpLF+u8euvwd/3ggWF62icOKHwySfSoFGepHSJAOvWaUyZYuXwYZVv\nv9UYNMgW0vS8+aaZ1q2j6do1it27w/8UatbMy7BhLm65JXi1hKVLTYwda2XlyspXYH3xhUqHDtHc\nfbeNyy+HzZuDew5ceGHRpsPatSP/GYY//1S45RY7jRrF8NBDtrBehjX8r+hKTtfhjTeKrkD211+h\nq07v2KHyxBM2fv9dZccOjQEDKs5KVf6aPdvMAw8YD9Hdc4+dRYsqV2B45x2Lb4pup9O4STjpyy9V\nunSJokOHKFasCOx7efXVPNq1c1O3rpeRI/Pp0iXyhx2PHm1l2zYTmZkKCxaYefttc/FvCpHKdTZH\noLlzzaxdW/QEKu0MlaVx4EDR+4g//qh87b1r12qnbZvo1q3y9FVUraqfdbugAPr2tXP8uHGOPPig\njSuvdFCnTsk6iuvX11m2rGJNjnfokHradvheN1JTCHOnT+PQoIGHp54KXbt4crKnyGIxqamVozDM\nyTFWtsvOhqZNizZnnL5d0Y0cmU+bNm5UVaddO3zzFmVn4wsIAAUFCunpUsQApKW5fM/fREfrYX3d\nSE0hzKWkePjPf3TcbiM49O1bvrUEXYdNmzS8Xrj6ag/Vq+usWZPLihUm4uP1SnGH/MsvCrffHsXh\nwyo1aniZP99Jbq7CV1+ptG3rYfDgytV5nZAAK1Y40TSVhIRo3xTtCQnGmhIbNxrFSoMGXnni/G+9\nerlp0CCXn37SaN/ezUUXBW+Y7bZtGrpuTCAZjMkKJSiEuTZtjPVw16830aSJN6idpf4YNMjG4sVG\n89WNNxprECcl6TzwQBj3lAXZv/9t9U2/cfiwyvTpFqZPL/4BuM8/1/jgAzOJiV4eeaT0axCEm7MV\nQHPnOnnvPTN5edCzp9uvPO/erfL441ayshQGD3bRs2fFvNFo1cpLq1bBrVU+/LCNDz4wrs/bbitg\n5szSP5hZqYPCsWMKEyZY+fNPhb593dx3X9l91qpVJubMMVO9upcxY1xUq+b/ncJVV3m56qoz70Y3\nbND4+muN1q09ZTKF9P79ii8gAHz0kZmff3ZVuhktVbXosfJn3YTdu1V69bLjchkl5y+/qLzzTvk/\nSe12w7vvmjl4UOH2291cfnnZHjubDe69t2Q3DGlpdvbvN4LusGE2mjfPrXRNcoE4cEDxBQSAZcvM\nPPaYi0svLd1+K3VQGDTIxqefGl/Bhg0azZpBkybB/5zvvlO5914bHo9RQOzfr7J4cek60pYsMTFw\noDHyR1V1Zs1y0rlzyQPD6tUa336rkZzsITm56PujokDTdF+6VVX3LclYmQwf7mLjRhMHDqjUru3l\nkUeKn2Zjxw7NFxAAtm4NzQo8jz9uZfZsCwBvvmnh449zwyqou91G4XaS16uwf79C06YhTFSEsNuN\na9LrNb4/RdGJitKB0rUhVepeoF27CrOv6wrffFN2n3OyYAX45pvSFxBLlxbGc69XYfnykg9xmz3b\nTFpaFC+9ZKVbN/sZo2qqVdOZNCkfq1XHYtF55pl8atWqGEGhoMDoLzmd1wuvvGIhLc3GK69Y8HqN\n0TCbNzvYti2HzZsd7NmjMmiQjUmTLOSd4+a/eXOPb+1mgBYtQlMQf/xx4XmSm6uwcWN4LQ9nMhkr\n2p2UlOSlVSvph/BHtWo6zz2Xj8mko2k6Tz2VT+3apb8+A6op6LqOw5GNyaShaR4cDiceT/jcffgr\nObmAZcuMr8Bi0WnRQimTvFx2mYLV6iI/3wgM7dq5yc31f1K3s6ld2wUUBoLYWBfZ2QVoGmiaiqZ5\nyMpysmaNgssFnTt7sJ62iuLy5VbAaJbyemH58gLatSvaTHXnnZCaavyuqsZ6veXtZH6CdWzGj7fw\n7rsmYmPh3//Op0OHwkLojTdMPPuscWGtXg0eTz4PPGAUWomJsGWLSlqaDTD+duBAAZMnn9m016gR\nzJjhYMECE4mJOo8+6iI3F3JzHUDhz7LWsKGbQ4cKA0G9ek5yc4NzfgfruPzf/+WQnGwiK8to4rLZ\n9HI/z4J9jpWXXr2gWzfjBsdiMa7P8+UlK0unSpUqKOfpkVZ0/Wz3S+eXlZVFnEysLoQQESczM5PY\n2Nhz/j+goKDrOvv2/YnJpBEbaycrK7Ki69lomho2eRk92sJ77xm1AEXRmT8/75yjFoYNs/pqOwD9\n+xfwzDNucnPtNG5c9LXLlztp3rxwP7m58NRTFnbtUmnf3suTT7rCcvH5YB6bNWs0Bg4snCbEYtHZ\nsyfXN5Jm9mxj+oqTnn46n7S0wuaNzZtVeve2cbLdtlu3AqZM8X9Iam6ug0svbcTu3T8TFRXZw5HC\n6ZoJxCefaEybZsZqhfHjC+jQwRaxeTnV+Y5LfHx0sTWFgJqPFEUhOroKJpNKbGw0Ho+G2x3ZX2Q4\n5WX79ijAKJ11HXbutNGhw9kLHl23cWozktfrIjq6gKpVo7HbdZzOwk7iWrXUvzuiDFFR8Nprp+7N\nEuScBEcwj01KCjRvHsW33xrf79Ch+URHF/YtDRwIdruZ7ds1WrXy0L9/0YunUyeYOdPEsmUm6tfX\neewxhaiokn9vUVHRREXFlCovoRZO10xJpacrDBkS7WvSvecenYMHlYjMy+nOd1xiY6sU//6ySpgI\nXIsWXn75pfCW/corz93xNnSoMTImI0OhRg3v3w9SKcTEwMyZ+YwcaaGgQGHs2Hzq1i37TuIFC0x8\n+aVGmzaesHxqMyoKli/PZdUqjbVrTWRnKxw4oBSZiiEtrYC0tHMPq7z9dvcZK8KJyJKervgCAhjD\n00+cAHP4TklUbiQohKHJk/NISNDZt0/hllvcdOx47qDQvLmXrVtz+O03lYYNvRhNhcbJftNNHlJS\n/OvQzM6GXbs0atf2lniumpNmzTLz6KNG08zbb4PL5aRXr/ArPC0WmDbNyp49RuD98EMTGzc6iIns\nG3dRAs2aebn4Yi979xq1xPbtPVSvrnHiRIgTFgYkKISh6GiYONH/0Unx8RAfH3iV99AhhZtuimL/\nfhWrVee//3WSklLyYYGffqqdtm0Ky6Dw55+KLyAApKer/PqryhVXRHazgfBfdDR8+GEu779vxmrV\nuftuD4oS2X08wVKpn1MQhrlzzb4nSvPzFd+axiV12WVFC9VmzcKzkK1eXScxsTBtcXE6deuGZ1pF\n2alWTefhh1088EABUVGhTk34kJqCKLJc5dm2/TV0qIv8fPjyS2PqjXCdKM5mgwULnLzwggW3W2HE\niHzi40OdKiHCgwQFQf/+BaxZY2LLFhPVq3tL1HR1KpOpcBrlcNekiTckcxEJEe4kKAiio2HZMicZ\nGRAb69+Eb0KIikmCQgT57DMNpxOuu+7MKSuCoTRNKJmZMHOmhfx8uPvugqDMwSKEKH+l7miePn06\nrVu34OKLLyQ5+SrefHNGMNIVEn/++Sd9+/aifv0kLrusEc8+OyHUSfJ55BEr3btHkZYWRbdudlx+\nttI4HA5atLiUYcMGl1nadB169IjixRet/PvfVm6+OYqsrODtf8WKZVx55ZXUrVuT5OSrmDPn3eDt\nXATkjz8OkJZ2F9WqVaNp0wYMHTqI7OwgHvQQeOSRR7jgguIf7gp3U6ZMplatWtStW5Pu3W/jwIH9\nJXp/qYLCJ5+s4fHHH2fGjDfZuzedV175DxMnjmft2o9Ls9uQueOOO6hXrz67d//KsmUf8fnnn7Fp\n08ZQJ4vsbJg7t/Cp2e3bTXz1lX9tPJMmPYvDkXPW/331lcr775v4/ffAptpdvNjE8OFWXnnFzM6d\nhek5eFDlhx+C0wa1c+cOHnzwPp555hl+//0gTz/9HP/610i2b98WlP2XlUOHFJ55xsKzz1o4fDh8\n1+MNVN++PalaNYEDBw6wbt3n7NnzA+PHjwl1sgL23XffMnv27PNO/xAJ3nxzJgsXzmfDhg3s3v0L\njRs3YcaM6SXaR6maj7755msuv/xyWrRoidvtxWptTWzsK7zzjpWOHSOrbXrTpo389ttvrFixBlBp\n0KARq1atC3WyALBaISpKJze38IQ9ffH0s/n++10sXbqInj37kJWVWeR/CxaYGDLEhq4rxMTorFiR\nW6IhpEuWmHjwQbtvOybGS06OcY8RFaVTr15whnhmZGQwYsRj3HzzzWRkOPjnPztz6aXN2Lp1E61b\ntwnKZ5xq0yaN5ctN1K3r5YEHCgJ6wjUvD267Lcr3YNTKlSbWr88tkya/UMjKyqRFi5Y89dR47HY7\nSUlJ9OjRO2JbCXRdZ+TI4YwcOZIxYyI3sAH85z+v8uyzL9CwYUMyMhw888ykEu+jVDWF66/vxO7d\nu9m0aSM//uiha1crf/11D2vW3MyoUaG7AjIy4JtvVHLOfoN8Vtu2beXyyy9n4sTxNG16Ea1bX8Hr\nr79adoksAYsFpk/Po2pVHZtNZ+zYfJo0Kb7QHTlyOKNHP3XWGRFnzTKj60aQyclRWLiwZKXf558X\njfiXXurln/90k5zsZvZsJzVrBqdP4frrOzFy5Cjftsfj4dChQ9SsmRSU/Z9q506V7t3tvP22hQkT\nbIweHdg5vG+f6gsIAL/8ohVZSCbSxcbG8fLLr1KtWnXf39LT/yApqVYIUxW4d955E7vdRu/evUOd\nlFL5668/2b9/H8ePH6dZs2Y0bFiXAQPSOHbsWIn2U6qg0LLlP5gyZQqpqbfSocOT5OUVVjzWrAlN\nH/aOHSqtW8dwww3RtG4dzZYt/mXx4MF0Nm/eTI0aiezc+QMvvPASzz//NKtXf1TGKfbPTTe5+emn\nHPbty+Hhh4vvUJgxYwaaptKrV5+z/r9GDf207ZLd2V95ZdHXJyd7eP99J0uWOLnmmrJbJGXChLFE\nR0dz++3dgr7vzz834XYXFt6ffRbYOVyrlpf4+MLvNyHBG7QgGY527vyKt96aySOPPBbqpJTY4cOH\nefHF53nppWmhTkqpHTyYDsDy5UtZt24dGzdu4+DBdEaOHFqi/ZSq5N648TOeeOIJFi1ahtPZnp49\nC//XqFGv3YZdAAAWAklEQVRonhCdMsVKZqZxYR89qnLHHVF+LVWp6zo1atTgoYcexu32cv31N3Dj\njTezbNliunS5sTyS7hd/mjyPHDnMuHHjWLLkw3O+5pln8vnrL5U9e1Q6dXIzYEDJ1tXt16+AnBzY\nuNFE8+YeRo4s++cTxo8fy7Jli1myZCUWS/BndL38cs95t/1VpQrMm5fLCy9YURR44on8Cjuv0qZN\nm7jzztsZO/Zprr762lAnp8TGjRtNnz5pNGp0CVlZR0OdnFI5uQrCsGGPkJiYiMUSw6hRo+nd+05c\nLpff10yJgsKcOXO4//77fZ0xKSld6datG1dffS1ut5fJk/N46aW/UJSDvP568xJmKThMpqJ3ZF6v\nwquvWujcueiayAsWzGPkyKG+vAwbNoKqVasWeU2dOvX46qsdZZvgIDg9L7feejv9+/enSZOm55wG\nOClJZ+XK0i1vNWhQAYMGlSyYFOdkXsCYon3fvkPouk7//v3Ztm07K1d+Qu3adYL6mSd17Ojh//7P\nyeLFZurW9TJuXOCr47Vo4eWDD0q3Dne4W736IwYNup/Jk6dyxx3dQ52cEtuw4VO++GI7U6e+AhQW\nqpGqRo1EwGjeO6lu3brous7Ro0eoVetCv/YT0CI7AEeOZDNgQD8SE6szefI0X+Hz+OMj+PPPP5k1\n6/1AdltqP/ygcuONUTgchbfUnTq5ee+981+gH364lKFDB/Pjj3uxWIyZPh988F5sNjvTppWs9z7U\nEhPjiI+P9wUJp9OJ1+slJiaG3bv3hjh1JTdmzCh27tzBggVLiYqK7CGDubk51K9fi99/PxjR6yls\n376Nfv16smDBfFq2bBuRaxAMGzaYZcuWYLcb17uu62RkZHDBBRfw/PMvcdttqSFOYcl4PB4aN67P\nc89N4qGHBpKR4eDjjz+mf/+72LfvEKqqUr168ddPqfoUuna9iYULF7J16xY8Hg87d+5g+fIl3HTT\nLaXZbak0beplx44crrrKmJ2zbl0vEyYUf8fXpcuNxMfHM27ck+Tm5rJx42esXv0Rd93Vr6yTHHS7\ndv3Ed999x4YNW1m/fjNpaffSpcuNrF+/OdRJK7Ft27ayYMF8Vq5cWeQOSISOx+Nh5MiHGTfuaf75\nz3+GOjkBmzjxebZu/Yr16zezYcNWPvrI6D9cv34zKSnh02TsL03T6N27H1Onvsivv/7KkSOHmTp1\nMt2790JV/S/qS9Wn0KtXb9zuPIYNG8zBg3+SlJTEsGEj6dkztL34CQnw0UdOcnPxe/ZDm83G6tWr\nGTDgPpo2vYhq1arz4ovTaNOmbdkmtgwkJSURHx+N3e7A7fZSpUoVMjOjSEysGeqkldi8eXPIzs6i\nXr16Rf7etm17PvhgSYhSVbl98cV2fv75J5544jH+9a9HURQFXddRFIXNm3dw4YW1Q51Ev8TGxvlu\nNEwmlehoM4qiROR1ctKYMePxeApo3bo1BQVubrnlthIPSy1V85HJpBIfH01GhiMiq4+nkryEr4qU\nn4rSfAQV67hUlryUefOREBXF1q0a771nYt++8w/vcrvdjBv3JDVrVmX9+rXllDohyo9MiCcqvVOX\nEY2N1fnoo1wuueTMu8Xc3Fy6dbuZxo2blncShSg3UlMQYeHYMYUhQ2ykptqZN69871VmzSp8mjsr\nS2Hp0rN/vsPhoHfvNKZNmx7xwxeFOBepKYiw8NBDNtatM07HTZs06tVz0q5d2T0Zfaozn+4+e4Ff\nvXp1+vW7uxxSFN7WrdOYMMGKoihMmQJtgj8FlQihgIOCqipomlHROPkzkkleQmv37sK06rrCjz9q\nXHONUTiXdX6mTHExYIDCL7+o3Hijm7vv9mAyFf9Zmqb69brT3xPoe8PBiRNw77123+SM3brB99+r\nxEX4aOFIvGbOpbR5CTgoJCRE+x6Oio21F/PqyCF5CY2uXeHtt43frVbo2tVKfHzRCenKKj/x8bBp\nE0ycCN9+a2bxYjMDBxb/vipVbMTHR5foszTNqP3ExtqJjS3Ze8PB4cOQe8qD8Lm5kJdnp379M197\n4ADMmGEMC3/4YWP6j3AXSddMcQLNS8BB4fhxB2azRmysnawsJx5PZA/j0jRV8hJCL7wA9eubOHhQ\nJTXVTe3aXjIyjP+VR37GjrUwfbrRt7B8OWhaHh7PXIYPHwIYU26kpx8p8p7s7DwyMhwl+hyHw3iy\n3shLBM0t/7eEBGjZ0uZbz6NVK6he3UlGRtHjkpkJycl2Dh407laXL/ewalX4rokdidfMuZwvL/7c\nxAQcFLxe3feBHo834sf2niR5CQ1FgUGDCifVc7vPfE1Z5ufLL5XTtlUmTuxBamqPU9JU9LMDSU9F\nuGYWLszlgw/MqKrC4MFWCgrOzMvOnZovIABs26Zx9KiX06YXCzuRfFxOF2hepKNZCKB9ew/bthVe\nDm3blk8ndySKiYEBAwowmVRiYqy+Gt2p6tf3YrPp5OUZwTYpyctZlvUQYSjye1WECIJRo1yMG5dH\njx4FzJjh5KabzqyqLFgwj7p1a1CvXiKKopCW1ot69RIZOXJYCFIc3mrX1nnnHSdt2ri57jo38+Y5\nKcH0OyKEpKYgBMbSsQ89VACceyrw7t170b17r/JLVIS7/noP119fsacPr4gkdgshhPCRoCCEEMJH\ngoIQQggfCQpCCCF8JCgIIYTwkaAghBDCR4KCEEIIHwkKQgghfCQoCCGE8JGgIIQQwkeCghBCCB8J\nCkIIIXwkKAghSmTfPoXt21WcZTDXna7D1q0aW7Zo6GdfKluUMQkKQkSwPXsUli3TyMwsn89buNBE\nu3bR3HxzNDfcYA/65w4aZOPWW6O47bYoHnzQFtydC79IUBAiQr32mplrronm/vujuPTSGL77Tin+\nTaU0aZIVt9v4nN27Vd5/P3j73rtXYfFis297yRIze/eWfZ5EURIUhIhQL75oBYxCs6BA4bHHyv7O\n2mYr2qZjD+I699HRoKqF+1dVnejilxQWQSZBQYgS+PDD5Vx3XTIXXVSL5OSrmDVrFosWmZg1y1xu\nTTgnaVrRbaUcbqpfeCGf2Fij4L7hBje9ewdv34mJOs88k4/JpGMyGb8nJkZGx0J2Nvz4Y9n0s5Q3\nWXlNCD/t3LmDhx66nzfeeIdOnVJYv/5/9Onjxes1bpdnzjSzenUuMTHlk57nn89jyBAbuq5gNuu8\n8EJ+mX9mcrKH77/PITtboWZNBbP5zCLk8GGFJUtMxMXp3HmnG1MJSpn77iugf39j9TuzuZgXh4nv\nvlPp0cPOsWMq9ep5WbYsl1q1IiOYnY3UFITwU0ZGBsOHP0rnzl1RVZV27Trj9Xb3/f+nnzR27Ch6\n+75pk0bHjlEkJ0exYkVw78G6d3eze3cOixY5+P77HK64whvU/Z+L1QrVqp290DtxArp2jWLsWBtD\nh9oZPLjkTVpmc+QEBIAXX7Rw7JhRlO7bp/Kf/1hCnKLSkZqCEH66/vpOXH99J9+22exBUY6j6wmA\n0QZes2ZhYel0Qlqanexso11n0CAbGzYEt43nggvgmmvKJxj4Y/t2jQMHCu81ly0zMX16ZBXypVUe\nzXhlKeCgoKoKmmYc/JM/I5nkJXyFa37Gj3+SCy88it0+m9xchUcfLaBZMzhZAc/JUXwBAcDlUjhy\npDAvJlN45aekznZc6tUDRdHRdSPfSUk6dnv457M059jo0W6+/NLE0aMKF1/sZcgQd0iPbWmvF0XX\nA3tERNd1lEgPiUIE6PHHH2fu3Ll8+umnNGzY8Kyv0XXo1AnWrTO2mzeHTz7JIjExjszMTGJjY8sx\nxeXnv/+FyZMhLg5mzICWLUOdorKXkwPp6VC/vtG8FskCDgrHjuVgNmvExtrJynLi8YRPFTYQmqZK\nXsJUqPIzf/48hg8fAoCiKKSnH0HXdR56aCA7d37FggVLqF27znn3kZcH8+aZKCiAnj3daFoOderU\n5MCBv4iOLqce6TJSkc6zypKX+Pjix/gG3Hzk9eq+D/R4vLjdkf1FniR5CV/lnZ/U1B6kpvbwbbvd\nXkaPfow9e35k5cpPiI2NKzY9JhP07evybefmyjUTziQv0tEshN+2bdvKokXz2bRpB7GxcaFOjhBl\nQoKCEH6aN28O2dnZ/OMfzYr8vW3b9nzwwZIQpUqI4JKgIISfXn75VV5++dVQJ0OIMhX+Y8WEEEKU\nGwkKQgSJrsOOHSo7d8plJSKXnL1CBMngwTa6do0mJSWaESMifLC6qLQkKAgRBD//rLJoUeFcDnPm\nWDhwQB7uFJFHgoIQQXD6OgOqqkf8k62icpKgIEQQ1KmjM2ZMPoqio6o6EybkU6NG5E6fLCovGZIq\nRJAMHepiwAAXigJRUaFOjYhUHo+xNsdPP6l07uyha1d3uX6+BAUhgkiWjzTk58P8+Wby8+HOOwuo\nWjXUKYoczz1n4ZVXjLbHuXNh3rxcrr/eU26fL0FBCBF0/frZ+fRTo3h5910za9bkSu3JT59/XrRY\n3rRJK9egIH0KQoigOnJE8QUEgD17NHbtkqLGX82bFw0A5bWi3klSUxBCBFVcnE7VqjonThhDcs1m\nnaQk6XT319NP5xMVZQxzvuEGN7feKn0KQogIZrHArFlOnnjCSn4+PP64izp1JCj4y26HCRPyQ/b5\nEhSEEEHXtq2H9etz/X59VhYMHGhn+3aN1q09zJjhpIIuTBf2pKFPCBFyL75oZe1aE9nZCmvXmnjx\nRXnyL1QkKAghQu7IEeW826L8SFAQQoTcXXcVYLEY/Q5ms85ddxWEOEWVlwQFIUrgrbfeoF27llx0\nUS3atm3Ba6+9EuokVQgdOnj4+ONcpk1z8sknuXToUH7j8kVR0tEshJ9WrVrJiy8+x7x5i7niihZs\n3bqFnj1vp0GDhqSkdA118iLepZd6ufTS8h2TL84kNQUh/FSrVi1mznyHK65oAUDbtu1o1KgxP/64\nO8QpEyJ4pKYghJ9OBgMAt9vNRx+tYP/+30lJuTGEqRIiuAIOCqqqoGlGRePkz0gmeQlf4ZafKVMm\n88ILz3LBBRfw2mszueyyZn6/99S8mEzhkZ9AhdtxKQ3JSyFF1/WAHjXUdR1FkWFjonJyu92sXr2a\nu+++mzlz5tClSxe/3peVlUVcXByZmZnEytNZIgwFHBSOHcvBbNaIjbWTleXE44nsDiJNUyUvYSpU\n+Zk/fx7Dhw8BQFEU0tOPnPGaESOG8tdff/Hee/P92qfDkUOdOjU5cOAvoqNjgpre8laRzrPKkpf4\n+OLndg+4+cjr1X0f6PF4cbsj+4s8SfISvso7P6mpPUhN7eHbHjlyOFWqxDJmzPhTXqWgaSa/0yXX\nTHiTvMjoIyH81r791bzzzpts3vw5Xq+X7du3sWTJIhmOKioUGX0khJ9uuy2VzMxMhg4dxNGjR0hK\nqsWIEaPo1atPqJMmRNBIUBCiBNLS7iEt7Z5QJ0OIMiPNR0IIIXwkKAghhPCRoCCEEMJHgoIQQggf\nCQpCCCF8JCgIIYTwkaAghBDCR4KCEEIIHwkKQgghfCQoCCGE8JFpLoQQogJYs0ZjyxYTrVt76d8/\n8P1IUBBCiAi3dKmJBx6wA/Daa+B2Q2pqYPuS5iMhhIhwH39c9P7+ww8D35cEBSGEiHBNmhRdTKeZ\n/8uGn0Gaj4SIQA4HzJljJj9foXfvAqpVC2hVXVFBDB7s4sQJ2LLFxFVXeXnqKTMOR2D7kqAgRITR\ndbjrLjtbtxqX73vvmVm71kF08cvvigrKZIKnnnIBLkwmFYsl8KAgzUdCBMDhcNCixaUMGza43D/7\n6FHFFxAA9u5V2b1bLmURHHImCRGASZOexeHICclnx8XpXHBBYRuy1apz4YXSfCSCQ4KCECX0/fe7\nWLp0ET17hmZtZosF5s510qqVh+bNPbz5ppNatSQoiOCQPgUhSmjUqEcYPfopDhzYT1ZWZkjS0LKl\nl5Urc0Py2aJiCzgoqKqCphkVjZM/I5nkJXyFU37eeedNTCaNvn37MWnScyiKgsnkf7pOzUtJ3heO\nwum4lJbkpVDAQSEhIRpFUQCIjbUHupuwI3kJX6HOz+HDh5k06TnWrVtHfHw0drsFq9VEfLz/w340\nzQMYeYmNrRjDhUJ9XIJJ8lKKoHD8uAOzWSM21k5WlhOPx1v8m8KYpqmSlzAVqvzMnz+P4cOH+G5+\nbr31dnr27E1SUj0yMhw4nS7y891kZPg/9s/hcAL8nRetTNJdXirSeVZZ8uLPDUzAQcHr1X0f6PF4\ncbsj+4s8SfISvso7P6mpPUhN7eHbTkyMo2rVqsydOwsAp9OJ1+vl449XsXv3Xr/2KddMeJO8SEez\nEH775psfi2y/9tor/PXXQSZOfCFEKRIi+CQoCOGnmjWTimxXqVKFzMwoEhNrhihFQgSfBAUhAvTY\nY0+EOglCBF3kj78SQggRNBIUhBBC+EhQEEII4SNBQQghhI8EBSGEED4SFIQQQvhIUBBCCOEjQUEI\nIYSPBAUhhBA+EhSEEEL4SFAQQgjhI0FBCCGEjwQFIYQQPhIUhBBC+EhQEEII4SNBQQghhI8EBSH8\nNG/eXGrWrEq9eonUq5dI3bo1qFcvka+//irUSRMiaGTlNSFKoH37q1m8+MNQJ0OIMiM1BSGEED4S\nFIQogT/+OED37rdxySV1ad36ChYu/CDUSRIiqAJuPlJVBU0zYsrJn5FM8hK+wiU/iYk1aNToEp56\nagKXXNKYFSuWMWjQ/dSufSFXX32tX/s4NS8mU2Qfn3A5LsEgeSmk6LquB/JGXddRFCWgDxWioujV\nqxdWq5V3333Xr9dnZWURFxdHZmYmsbGxZZw6IUou4JrC8eMOzGaN2Fg7WVlOPB5vMNNV7jRNlbyE\nqVDlZ/78eQwfPgQARVFITz9yxmtq1ryQr7/eSUaGw699OhxOgL/zogUvsSFQkc6zypKX+PjoYt8f\ncFDwenXfB3o8XtzuyP4iT5K8hK/yzk9qag9SU3v4tt9887/Ex8dz6613+P62Z8+P1K1b3+90yTUT\n3iQvMiRVCL+5XPk88cRj1KtXn2bNLmf58iWsXfsJq1evD3XShAgaCQpC+On++wfhcDgYMKA/R44c\nom7desya9T6XX9481EkTImgkKAhRAsOHP8rw4Y+GOhlClJnIH38lhBAiaCQoCCGE8JGgIIQQwkeC\nghBCCB8JCkIIIXwkKAghhPCRoCCEEMJHgoIQQggfCQpCCCF8JCgIIYTwkaAghBDFWLnSRJ8+dkaM\nsHLsWMVeR0bmPhJCiPP4+muV++6z4fEYwWD/fpWFC50hTlXZkZqCEEKcx65dmi8gAHz7bWQvjlQc\nCQpCCHEerVp5sFoLVy2++mp3CFNT9qT5SAghzqNxYy8LFzqZP99EYqLOkCGuUCepTElQEEKIYrRp\n46FNG0+ok1EupPlICCGEj6Lrul78y4QQwZCVlUVcXByZmZnExsaGOjlCnEGCghDlSNd1srOzqVKl\nCopSsce7i8gkQUEIIYSP9CkIIYTwkaAghBDCR4KCEEIIHwkKQgghfCQoCCGE8JGgIIQQwkeCghBC\nCJ//B40Jke3pVIwWAAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# sklearnのPCAを使ってみる\n", "from sklearn.decomposition import PCA, KernelPCA\n", "pca = PCA()\n", "X_pca = pca.fit_transform(Xc)\n", "list_plot(X_pca[:, 0:2], figsize=4, zorder=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

KernelPCAを使ってみる

\n", "\t

\n", "\t\tsklearnのPCAの使い方が分かったので、今度はKernelPCAで同様の計算をしてみます。\n", "\t\tカーネルにはrbfを使いgamma=0.001と0.1の結果を出力してみます。\n", "\t

\n", "\t

\n", "\t\t今度は、井出本と同様の結果が求まりました。KernelPCAは、計算に使用する固有値が少ない場合、\n", "\t\tPCAよりも速く計算することができるので、データ数が多い場合には有効です。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# sklearnのKernelPCAを使ってみる\n", "# σ=0.001(通常のPCAに近い)、σ=0.1の違いをみる\n", "kpca = KernelPCA(kernel=\"rbf\", gamma=0.001)\n", "X_kpca = kpca.fit_transform(Xc)\n", "g0_001 = list_plot(X_kpca[:, 0:2], figsize=4, zorder=2)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
σ=0.001σ=0.1
" ], "text/plain": [ "<__main__.Table2Html instance at 0x7fcd74c10908>" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kpca = KernelPCA(kernel=\"rbf\", gamma=0.1)\n", "X_kpca = kpca.fit_transform(Xc)\n", "g0_1 = list_plot(X_kpca[:, 0:2], figsize=4, zorder=2)\n", "Table2Html([[_to_png(g0_001), _to_png(g0_1)]], header=['σ=0.001', 'σ=0.1'])" ] }, { "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 }