02_mt_method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anomaly Detection Examples with R MASS road dataset by MT-Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "井出剛著「入門 機械学習による異常検知〜Rによる実践ガイド」P50-52 の例"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### データの読み込みと確認"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'https://vincentarelbundock.github.io/Rdatasets/csv/MASS/road.csv'\n",
    "df = pd.read_csv(url, index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>deaths</th>\n",
       "      <th>drivers</th>\n",
       "      <th>popden</th>\n",
       "      <th>rural</th>\n",
       "      <th>temp</th>\n",
       "      <th>fuel</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Alabama</th>\n",
       "      <td>968</td>\n",
       "      <td>158</td>\n",
       "      <td>64.0</td>\n",
       "      <td>66.0</td>\n",
       "      <td>62</td>\n",
       "      <td>119.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Alaska</th>\n",
       "      <td>43</td>\n",
       "      <td>11</td>\n",
       "      <td>0.4</td>\n",
       "      <td>5.9</td>\n",
       "      <td>30</td>\n",
       "      <td>6.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arizona</th>\n",
       "      <td>588</td>\n",
       "      <td>91</td>\n",
       "      <td>12.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>64</td>\n",
       "      <td>65.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arkanas</th>\n",
       "      <td>640</td>\n",
       "      <td>92</td>\n",
       "      <td>34.0</td>\n",
       "      <td>73.0</td>\n",
       "      <td>51</td>\n",
       "      <td>74.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Calif</th>\n",
       "      <td>4743</td>\n",
       "      <td>952</td>\n",
       "      <td>100.0</td>\n",
       "      <td>118.0</td>\n",
       "      <td>65</td>\n",
       "      <td>105.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Colo</th>\n",
       "      <td>566</td>\n",
       "      <td>109</td>\n",
       "      <td>17.0</td>\n",
       "      <td>73.0</td>\n",
       "      <td>42</td>\n",
       "      <td>78.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Conn</th>\n",
       "      <td>325</td>\n",
       "      <td>167</td>\n",
       "      <td>518.0</td>\n",
       "      <td>5.1</td>\n",
       "      <td>37</td>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Dela</th>\n",
       "      <td>118</td>\n",
       "      <td>30</td>\n",
       "      <td>226.0</td>\n",
       "      <td>3.4</td>\n",
       "      <td>41</td>\n",
       "      <td>20.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DC</th>\n",
       "      <td>115</td>\n",
       "      <td>35</td>\n",
       "      <td>12524.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>44</td>\n",
       "      <td>23.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Florida</th>\n",
       "      <td>1545</td>\n",
       "      <td>298</td>\n",
       "      <td>91.0</td>\n",
       "      <td>57.0</td>\n",
       "      <td>67</td>\n",
       "      <td>216.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Georgia</th>\n",
       "      <td>1302</td>\n",
       "      <td>203</td>\n",
       "      <td>68.0</td>\n",
       "      <td>83.0</td>\n",
       "      <td>54</td>\n",
       "      <td>162.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Idaho</th>\n",
       "      <td>262</td>\n",
       "      <td>41</td>\n",
       "      <td>8.1</td>\n",
       "      <td>40.0</td>\n",
       "      <td>36</td>\n",
       "      <td>29.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ill</th>\n",
       "      <td>2207</td>\n",
       "      <td>544</td>\n",
       "      <td>180.0</td>\n",
       "      <td>102.0</td>\n",
       "      <td>33</td>\n",
       "      <td>350.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ind</th>\n",
       "      <td>1410</td>\n",
       "      <td>254</td>\n",
       "      <td>129.0</td>\n",
       "      <td>89.0</td>\n",
       "      <td>37</td>\n",
       "      <td>196.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Iowa</th>\n",
       "      <td>833</td>\n",
       "      <td>150</td>\n",
       "      <td>49.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>30</td>\n",
       "      <td>109.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kansas</th>\n",
       "      <td>669</td>\n",
       "      <td>136</td>\n",
       "      <td>27.0</td>\n",
       "      <td>124.0</td>\n",
       "      <td>42</td>\n",
       "      <td>94.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kent</th>\n",
       "      <td>911</td>\n",
       "      <td>147</td>\n",
       "      <td>76.0</td>\n",
       "      <td>65.0</td>\n",
       "      <td>44</td>\n",
       "      <td>104.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Louis</th>\n",
       "      <td>1037</td>\n",
       "      <td>146</td>\n",
       "      <td>72.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>65</td>\n",
       "      <td>109.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maine</th>\n",
       "      <td>1196</td>\n",
       "      <td>46</td>\n",
       "      <td>31.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>30</td>\n",
       "      <td>37.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maryl</th>\n",
       "      <td>616</td>\n",
       "      <td>157</td>\n",
       "      <td>314.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>44</td>\n",
       "      <td>113.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mass</th>\n",
       "      <td>766</td>\n",
       "      <td>255</td>\n",
       "      <td>655.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>37</td>\n",
       "      <td>166.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mich</th>\n",
       "      <td>2120</td>\n",
       "      <td>403</td>\n",
       "      <td>137.0</td>\n",
       "      <td>95.0</td>\n",
       "      <td>33</td>\n",
       "      <td>306.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Minn</th>\n",
       "      <td>841</td>\n",
       "      <td>189</td>\n",
       "      <td>43.0</td>\n",
       "      <td>110.0</td>\n",
       "      <td>22</td>\n",
       "      <td>132.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Miss</th>\n",
       "      <td>648</td>\n",
       "      <td>85</td>\n",
       "      <td>46.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>57</td>\n",
       "      <td>77.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mo</th>\n",
       "      <td>1289</td>\n",
       "      <td>234</td>\n",
       "      <td>63.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>40</td>\n",
       "      <td>180.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mont</th>\n",
       "      <td>259</td>\n",
       "      <td>38</td>\n",
       "      <td>4.6</td>\n",
       "      <td>72.0</td>\n",
       "      <td>29</td>\n",
       "      <td>31.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         deaths  drivers   popden  rural  temp   fuel\n",
       "Alabama     968      158     64.0   66.0    62  119.0\n",
       "Alaska       43       11      0.4    5.9    30    6.2\n",
       "Arizona     588       91     12.0   33.0    64   65.0\n",
       "Arkanas     640       92     34.0   73.0    51   74.0\n",
       "Calif      4743      952    100.0  118.0    65  105.0\n",
       "Colo        566      109     17.0   73.0    42   78.0\n",
       "Conn        325      167    518.0    5.1    37   95.0\n",
       "Dela        118       30    226.0    3.4    41   20.0\n",
       "DC          115       35  12524.0    0.0    44   23.0\n",
       "Florida    1545      298     91.0   57.0    67  216.0\n",
       "Georgia    1302      203     68.0   83.0    54  162.0\n",
       "Idaho       262       41      8.1   40.0    36   29.0\n",
       "Ill        2207      544    180.0  102.0    33  350.0\n",
       "Ind        1410      254    129.0   89.0    37  196.0\n",
       "Iowa        833      150     49.0  100.0    30  109.0\n",
       "Kansas      669      136     27.0  124.0    42   94.0\n",
       "Kent        911      147     76.0   65.0    44  104.0\n",
       "Louis      1037      146     72.0   40.0    65  109.0\n",
       "Maine      1196       46     31.0   19.0    30   37.0\n",
       "Maryl       616      157    314.0   29.0    44  113.0\n",
       "Mass        766      255    655.0   17.0    37  166.0\n",
       "Mich       2120      403    137.0   95.0    33  306.0\n",
       "Minn        841      189     43.0  110.0    22  132.0\n",
       "Miss        648       85     46.0   59.0    57   77.0\n",
       "Mo         1289      234     63.0  100.0    40  180.0\n",
       "Mont        259       38      4.6   72.0    29   31.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = df.to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "各列をdrivers列で割って運転者1人あたりの数値に変換"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "drivers = df['drivers'].copy()\n",
    "for label, content in df.iteritems():\n",
    "    df[label] = content / drivers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>deaths</th>\n",
       "      <th>drivers</th>\n",
       "      <th>popden</th>\n",
       "      <th>rural</th>\n",
       "      <th>temp</th>\n",
       "      <th>fuel</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Alabama</th>\n",
       "      <td>6.126582</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.405063</td>\n",
       "      <td>0.417722</td>\n",
       "      <td>0.392405</td>\n",
       "      <td>0.753165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Alaska</th>\n",
       "      <td>3.909091</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.036364</td>\n",
       "      <td>0.536364</td>\n",
       "      <td>2.727273</td>\n",
       "      <td>0.563636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arizona</th>\n",
       "      <td>6.461538</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.131868</td>\n",
       "      <td>0.362637</td>\n",
       "      <td>0.703297</td>\n",
       "      <td>0.714286</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arkanas</th>\n",
       "      <td>6.956522</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.369565</td>\n",
       "      <td>0.793478</td>\n",
       "      <td>0.554348</td>\n",
       "      <td>0.804348</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Calif</th>\n",
       "      <td>4.982143</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.105042</td>\n",
       "      <td>0.123950</td>\n",
       "      <td>0.068277</td>\n",
       "      <td>0.110294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Colo</th>\n",
       "      <td>5.192661</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.155963</td>\n",
       "      <td>0.669725</td>\n",
       "      <td>0.385321</td>\n",
       "      <td>0.715596</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Conn</th>\n",
       "      <td>1.946108</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.101796</td>\n",
       "      <td>0.030539</td>\n",
       "      <td>0.221557</td>\n",
       "      <td>0.568862</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Dela</th>\n",
       "      <td>3.933333</td>\n",
       "      <td>1.0</td>\n",
       "      <td>7.533333</td>\n",
       "      <td>0.113333</td>\n",
       "      <td>1.366667</td>\n",
       "      <td>0.666667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DC</th>\n",
       "      <td>3.285714</td>\n",
       "      <td>1.0</td>\n",
       "      <td>357.828571</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.257143</td>\n",
       "      <td>0.657143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Florida</th>\n",
       "      <td>5.184564</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.305369</td>\n",
       "      <td>0.191275</td>\n",
       "      <td>0.224832</td>\n",
       "      <td>0.724832</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Georgia</th>\n",
       "      <td>6.413793</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.334975</td>\n",
       "      <td>0.408867</td>\n",
       "      <td>0.266010</td>\n",
       "      <td>0.798030</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Idaho</th>\n",
       "      <td>6.390244</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.197561</td>\n",
       "      <td>0.975610</td>\n",
       "      <td>0.878049</td>\n",
       "      <td>0.707317</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ill</th>\n",
       "      <td>4.056985</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.330882</td>\n",
       "      <td>0.187500</td>\n",
       "      <td>0.060662</td>\n",
       "      <td>0.643382</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ind</th>\n",
       "      <td>5.551181</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.507874</td>\n",
       "      <td>0.350394</td>\n",
       "      <td>0.145669</td>\n",
       "      <td>0.771654</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Iowa</th>\n",
       "      <td>5.553333</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.326667</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>0.726667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kansas</th>\n",
       "      <td>4.919118</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.198529</td>\n",
       "      <td>0.911765</td>\n",
       "      <td>0.308824</td>\n",
       "      <td>0.691176</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kent</th>\n",
       "      <td>6.197279</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.517007</td>\n",
       "      <td>0.442177</td>\n",
       "      <td>0.299320</td>\n",
       "      <td>0.707483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Louis</th>\n",
       "      <td>7.102740</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.493151</td>\n",
       "      <td>0.273973</td>\n",
       "      <td>0.445205</td>\n",
       "      <td>0.746575</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maine</th>\n",
       "      <td>26.000000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.673913</td>\n",
       "      <td>0.413043</td>\n",
       "      <td>0.652174</td>\n",
       "      <td>0.804348</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maryl</th>\n",
       "      <td>3.923567</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>0.184713</td>\n",
       "      <td>0.280255</td>\n",
       "      <td>0.719745</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mass</th>\n",
       "      <td>3.003922</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.568627</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>0.145098</td>\n",
       "      <td>0.650980</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mich</th>\n",
       "      <td>5.260546</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.339950</td>\n",
       "      <td>0.235732</td>\n",
       "      <td>0.081886</td>\n",
       "      <td>0.759305</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Minn</th>\n",
       "      <td>4.449735</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.227513</td>\n",
       "      <td>0.582011</td>\n",
       "      <td>0.116402</td>\n",
       "      <td>0.698413</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Miss</th>\n",
       "      <td>7.623529</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.541176</td>\n",
       "      <td>0.694118</td>\n",
       "      <td>0.670588</td>\n",
       "      <td>0.905882</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mo</th>\n",
       "      <td>5.508547</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.269231</td>\n",
       "      <td>0.427350</td>\n",
       "      <td>0.170940</td>\n",
       "      <td>0.769231</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mont</th>\n",
       "      <td>6.815789</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.121053</td>\n",
       "      <td>1.894737</td>\n",
       "      <td>0.763158</td>\n",
       "      <td>0.815789</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            deaths  drivers      popden     rural      temp      fuel\n",
       "Alabama   6.126582      1.0    0.405063  0.417722  0.392405  0.753165\n",
       "Alaska    3.909091      1.0    0.036364  0.536364  2.727273  0.563636\n",
       "Arizona   6.461538      1.0    0.131868  0.362637  0.703297  0.714286\n",
       "Arkanas   6.956522      1.0    0.369565  0.793478  0.554348  0.804348\n",
       "Calif     4.982143      1.0    0.105042  0.123950  0.068277  0.110294\n",
       "Colo      5.192661      1.0    0.155963  0.669725  0.385321  0.715596\n",
       "Conn      1.946108      1.0    3.101796  0.030539  0.221557  0.568862\n",
       "Dela      3.933333      1.0    7.533333  0.113333  1.366667  0.666667\n",
       "DC        3.285714      1.0  357.828571  0.000000  1.257143  0.657143\n",
       "Florida   5.184564      1.0    0.305369  0.191275  0.224832  0.724832\n",
       "Georgia   6.413793      1.0    0.334975  0.408867  0.266010  0.798030\n",
       "Idaho     6.390244      1.0    0.197561  0.975610  0.878049  0.707317\n",
       "Ill       4.056985      1.0    0.330882  0.187500  0.060662  0.643382\n",
       "Ind       5.551181      1.0    0.507874  0.350394  0.145669  0.771654\n",
       "Iowa      5.553333      1.0    0.326667  0.666667  0.200000  0.726667\n",
       "Kansas    4.919118      1.0    0.198529  0.911765  0.308824  0.691176\n",
       "Kent      6.197279      1.0    0.517007  0.442177  0.299320  0.707483\n",
       "Louis     7.102740      1.0    0.493151  0.273973  0.445205  0.746575\n",
       "Maine    26.000000      1.0    0.673913  0.413043  0.652174  0.804348\n",
       "Maryl     3.923567      1.0    2.000000  0.184713  0.280255  0.719745\n",
       "Mass      3.003922      1.0    2.568627  0.066667  0.145098  0.650980\n",
       "Mich      5.260546      1.0    0.339950  0.235732  0.081886  0.759305\n",
       "Minn      4.449735      1.0    0.227513  0.582011  0.116402  0.698413\n",
       "Miss      7.623529      1.0    0.541176  0.694118  0.670588  0.905882\n",
       "Mo        5.508547      1.0    0.269231  0.427350  0.170940  0.769231\n",
       "Mont      6.815789      1.0    0.121053  1.894737  0.763158  0.815789"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 対数変換"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "各列に１を加えた上で対数変換"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>deaths</th>\n",
       "      <th>popden</th>\n",
       "      <th>rural</th>\n",
       "      <th>temp</th>\n",
       "      <th>fuel</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Alabama</th>\n",
       "      <td>1.963832</td>\n",
       "      <td>0.340082</td>\n",
       "      <td>0.349051</td>\n",
       "      <td>0.331033</td>\n",
       "      <td>0.561422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Alaska</th>\n",
       "      <td>1.591089</td>\n",
       "      <td>0.035718</td>\n",
       "      <td>0.429418</td>\n",
       "      <td>1.315677</td>\n",
       "      <td>0.447014</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arizona</th>\n",
       "      <td>2.009762</td>\n",
       "      <td>0.123869</td>\n",
       "      <td>0.309422</td>\n",
       "      <td>0.532566</td>\n",
       "      <td>0.538997</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Arkanas</th>\n",
       "      <td>2.073992</td>\n",
       "      <td>0.314493</td>\n",
       "      <td>0.584157</td>\n",
       "      <td>0.441056</td>\n",
       "      <td>0.590199</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Calif</th>\n",
       "      <td>1.788779</td>\n",
       "      <td>0.099883</td>\n",
       "      <td>0.116849</td>\n",
       "      <td>0.066047</td>\n",
       "      <td>0.104625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Colo</th>\n",
       "      <td>1.823365</td>\n",
       "      <td>0.144934</td>\n",
       "      <td>0.512659</td>\n",
       "      <td>0.325932</td>\n",
       "      <td>0.539761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Conn</th>\n",
       "      <td>1.080485</td>\n",
       "      <td>1.411425</td>\n",
       "      <td>0.030082</td>\n",
       "      <td>0.200126</td>\n",
       "      <td>0.450351</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Dela</th>\n",
       "      <td>1.596015</td>\n",
       "      <td>2.143980</td>\n",
       "      <td>0.107359</td>\n",
       "      <td>0.861482</td>\n",
       "      <td>0.510826</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DC</th>\n",
       "      <td>1.455287</td>\n",
       "      <td>5.882845</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.814100</td>\n",
       "      <td>0.505095</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Florida</th>\n",
       "      <td>1.822056</td>\n",
       "      <td>0.266486</td>\n",
       "      <td>0.175024</td>\n",
       "      <td>0.202804</td>\n",
       "      <td>0.545130</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Georgia</th>\n",
       "      <td>2.003342</td>\n",
       "      <td>0.288913</td>\n",
       "      <td>0.342786</td>\n",
       "      <td>0.235870</td>\n",
       "      <td>0.586691</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Idaho</th>\n",
       "      <td>2.000161</td>\n",
       "      <td>0.180287</td>\n",
       "      <td>0.680877</td>\n",
       "      <td>0.630233</td>\n",
       "      <td>0.534923</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ill</th>\n",
       "      <td>1.620771</td>\n",
       "      <td>0.285842</td>\n",
       "      <td>0.171850</td>\n",
       "      <td>0.058893</td>\n",
       "      <td>0.496757</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ind</th>\n",
       "      <td>1.879645</td>\n",
       "      <td>0.410701</td>\n",
       "      <td>0.300396</td>\n",
       "      <td>0.135989</td>\n",
       "      <td>0.571913</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Iowa</th>\n",
       "      <td>1.879974</td>\n",
       "      <td>0.282670</td>\n",
       "      <td>0.510826</td>\n",
       "      <td>0.182322</td>\n",
       "      <td>0.546193</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kansas</th>\n",
       "      <td>1.778187</td>\n",
       "      <td>0.181095</td>\n",
       "      <td>0.648027</td>\n",
       "      <td>0.269129</td>\n",
       "      <td>0.525424</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kent</th>\n",
       "      <td>1.973703</td>\n",
       "      <td>0.416739</td>\n",
       "      <td>0.366154</td>\n",
       "      <td>0.261841</td>\n",
       "      <td>0.535020</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Louis</th>\n",
       "      <td>2.092202</td>\n",
       "      <td>0.400888</td>\n",
       "      <td>0.242140</td>\n",
       "      <td>0.368252</td>\n",
       "      <td>0.557657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maine</th>\n",
       "      <td>3.295837</td>\n",
       "      <td>0.515164</td>\n",
       "      <td>0.345746</td>\n",
       "      <td>0.502092</td>\n",
       "      <td>0.590199</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Maryl</th>\n",
       "      <td>1.594033</td>\n",
       "      <td>1.098612</td>\n",
       "      <td>0.169501</td>\n",
       "      <td>0.247059</td>\n",
       "      <td>0.542176</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mass</th>\n",
       "      <td>1.387274</td>\n",
       "      <td>1.272181</td>\n",
       "      <td>0.064539</td>\n",
       "      <td>0.135490</td>\n",
       "      <td>0.501369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mich</th>\n",
       "      <td>1.834267</td>\n",
       "      <td>0.292633</td>\n",
       "      <td>0.211664</td>\n",
       "      <td>0.078706</td>\n",
       "      <td>0.564919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Minn</th>\n",
       "      <td>1.695567</td>\n",
       "      <td>0.204990</td>\n",
       "      <td>0.458697</td>\n",
       "      <td>0.110111</td>\n",
       "      <td>0.529694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Miss</th>\n",
       "      <td>2.154494</td>\n",
       "      <td>0.432546</td>\n",
       "      <td>0.527162</td>\n",
       "      <td>0.513176</td>\n",
       "      <td>0.644945</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mo</th>\n",
       "      <td>1.873116</td>\n",
       "      <td>0.238411</td>\n",
       "      <td>0.355820</td>\n",
       "      <td>0.157807</td>\n",
       "      <td>0.570545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Mont</th>\n",
       "      <td>2.056146</td>\n",
       "      <td>0.114268</td>\n",
       "      <td>1.062894</td>\n",
       "      <td>0.567106</td>\n",
       "      <td>0.596520</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           deaths    popden     rural      temp      fuel\n",
       "Alabama  1.963832  0.340082  0.349051  0.331033  0.561422\n",
       "Alaska   1.591089  0.035718  0.429418  1.315677  0.447014\n",
       "Arizona  2.009762  0.123869  0.309422  0.532566  0.538997\n",
       "Arkanas  2.073992  0.314493  0.584157  0.441056  0.590199\n",
       "Calif    1.788779  0.099883  0.116849  0.066047  0.104625\n",
       "Colo     1.823365  0.144934  0.512659  0.325932  0.539761\n",
       "Conn     1.080485  1.411425  0.030082  0.200126  0.450351\n",
       "Dela     1.596015  2.143980  0.107359  0.861482  0.510826\n",
       "DC       1.455287  5.882845  0.000000  0.814100  0.505095\n",
       "Florida  1.822056  0.266486  0.175024  0.202804  0.545130\n",
       "Georgia  2.003342  0.288913  0.342786  0.235870  0.586691\n",
       "Idaho    2.000161  0.180287  0.680877  0.630233  0.534923\n",
       "Ill      1.620771  0.285842  0.171850  0.058893  0.496757\n",
       "Ind      1.879645  0.410701  0.300396  0.135989  0.571913\n",
       "Iowa     1.879974  0.282670  0.510826  0.182322  0.546193\n",
       "Kansas   1.778187  0.181095  0.648027  0.269129  0.525424\n",
       "Kent     1.973703  0.416739  0.366154  0.261841  0.535020\n",
       "Louis    2.092202  0.400888  0.242140  0.368252  0.557657\n",
       "Maine    3.295837  0.515164  0.345746  0.502092  0.590199\n",
       "Maryl    1.594033  1.098612  0.169501  0.247059  0.542176\n",
       "Mass     1.387274  1.272181  0.064539  0.135490  0.501369\n",
       "Mich     1.834267  0.292633  0.211664  0.078706  0.564919\n",
       "Minn     1.695567  0.204990  0.458697  0.110111  0.529694\n",
       "Miss     2.154494  0.432546  0.527162  0.513176  0.644945\n",
       "Mo       1.873116  0.238411  0.355820  0.157807  0.570545\n",
       "Mont     2.056146  0.114268  1.062894  0.567106  0.596520"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.log(df[['deaths','popden', 'rural', 'temp', 'fuel']] + 1)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 中心化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "mx = np.mean(x, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "deaths    1.858592\n",
       "popden    0.668448\n",
       "rural     0.348965\n",
       "temp      0.367111\n",
       "fuel      0.526476\n",
       "dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "xc = x - mx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 共分散行列の計算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sx = np.cov(xc, rowvar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 異常度の計算とプロット"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.sum(xc.dot(np.linalg.inv(sx)) * xc.to_numpy(), axis=1) / len(xc.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEeCAYAAACqgP9OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcHFW99/HPLxuLuMQQFA1MjIpyBVky4ChyZREvKLigAhH0yiOC4gIu+FxUBMSLzxXxomwaQRCJiAgoREAQwyYGmIGwg+BgAEUTYhBQBML8nj9OdVLTU1Vd1dM93Wfm+369+jXT3aeqTlVX//rU2crcHRERicekTmdARESqUeAWEYmMAreISGQUuEVEIqPALSISGQVuEZHIKHCLiERGgVtEJDIK3CIikZnSjpWuv/76Pnv27HasWkRkXBoYGHjU3WeWSduWwD179mz6+/vbsWoRkXHJzJaWTauqEhGRyChwi4hERoFbRCQyCtwiIpFR4BYRiYwCt4hIZBS4W2hg6UpOXnQ/A0tXdjorIjKOtaUf90Q0sHQl+562mGdWDTFtyiQWHNDH3J7pnc6WiIxDKnG3yOLBFTyzaoghh2dXDbF4cEWnsyQi45QCd4v0zZnBtCmTmGwwdcok+ubM6HSWRGScUlVJi8ztmc6CA/pYPLiCvjkzVE0iIm2jwN1Cc3umK2CLSNupqkREJDIK3CIikVHgFhGJjAK3iEhkFLhlGI3+FOl+6lUiq2n0p0gcVOKW1TT6UyQOCtyymkZ/isRBVSWymkZ/isRBgVuG0ehPke5XuqrEzCab2S1mtrCdGRIRkWJV6rgPAe5uV0ZERKScUoHbzGYB7wBOa292RESkkbIl7hOALwBDbcyLiIiU0DBwm9nuwDJ3H2iQ7kAz6zez/uXLl7csgyIiMlyZEvd2wDvN7I/AT4CdzOzs+kTuPt/de929d+bMmS3OpoiI1DQM3O5+uLvPcvfZwD7Ab9x9v7bnTEREMmnkpIhIZCoNwHH3q4Cr2pITEREpRSVuEZHIKHCLiERGgVtEJDIK3CIikVHgFhGJjAK3iEhkFLhFRCKjwC0iEhkFbhGRyChwi4hERoFbRCQyCtwiIpFR4BYRiYwCt4hIZBS4RUQio8AtIhIZBW4RkcgocIuIRCbqwD2wdCUnL7qfgaUrO50VkajouxO3Svec7CYDS1ey72mLeWbVENOmTGLBAX3M7Zne6WyJdD19d+IXbYl78eAKnlk1xJDDs6uGWDy4otNZEomCvjvxizZw982ZwbQpk5hsMHXKJPrmzOh0lkSioO9O/MzdW77S3t5e7+/vb/l66w0sXcniwRX0zZmhSz2RCvTd6T5mNuDuvWXSRlvHDTC3Z7pOOpEm6LsTt2irSkREJioFbhGRyChwi4hERoFbRCQyCtwiIpFR4BYRiYwCt4hIZBS4RUQio8AtIhIZBW4RkcgocIuIREaBW0QkMgrcIiKRaRi4zWxtM7vRzG41szvN7OixyJiIiGQrM63r08BO7v6kmU0FrjOzS919cZvzJiIiGRoGbg93WngyeTo1ebT+7gsiIlJKqTpuM5tsZkuAZcAV7n5DRpoDzazfzPqXL1/e6nyKiEiiVOB29+fcfUtgFrCtmW2WkWa+u/e6e+/MmTNbnU8REUlU6lXi7o8Bi4Bd25MdERFppEyvkplm9qLk/3WAXYB72p0xERHJVqZXyYbAD81sMiHQ/9TdF7Y3WyIikqdMr5LbgK3GIC8iIlKCRk6KiERGgVtEJDIK3CIikVHgFhGJjAK3iEhkFLhFRCKjwC0iEhkFbhGRyChwi4hERoFbRCQyCtwiIpFR4BYRiYwCt4hIZBS4RUQio8AtIhIZBW4RkcgocIuIREaBW0QkMgrcIiKRUeAWEYmMAreISGQUuEVEIqPALSISGQVuEZHIKHCLiERGgVtEJDIK3CIikVHgFhGJjAK3iEhkFLhFRFpgYOlKTl50PwNLV7Z9W1PavgURkXFuYOlK9j1tMc+sGmLalEksOKCPuT3T27Y9lbhFREZp8eAKnlk1xJDDs6uGWDy4oq3bU+AWERmlvjkzmDZlEpMNpk6ZRN+cGW3dnqpKRERGaW7PdBYc0MfiwRX0zZnR1moSUOAWEWmJuT3T2x6wa1RVIiISmYaB28w2MrNFZnaXmd1pZoeMRcZERCRbmaqSVcDn3P1mM3s+MGBmV7j7XW3Om4iIZGhY4nb3R9z95uT/J4C7gZe3O2MiIpKtUh23mc0GtgJuaEdmRESksdKB28zWA84HDnX3xzPeP9DM+s2sf/ny5a3Mo4iIpJQK3GY2lRC0F7j7BVlp3H2+u/e6e+/MmTNbmUcREUkp06vEgNOBu939W+3PkoiIFClT4t4O+CCwk5ktSR5vb3O+REQkR8PugO5+HWBjkBcRESlBIydFRCKjwC0iEhkFbhGRyChwi4hERoFbRCQyCtwiIpFR4BYRiYwCt4hIZBS4RUQio8AtIhIZBW4RkcgocIuIREaBW0QkMgrcIiKRUeAWEYmMAreISGQUuEVEIqPALSISGQVuEZHIKHCLiERGgVtEJDIK3CIikVHgFhGJjAK3iEhkFLhFRCKjwC0iEhkFbhGRyChwi4hERoFbRCQyCtwiIpFR4BYRiYwCt4hIZBS4RUQio8AtIhIZBW4RkcgocIuIREaBW0QkMg0Dt5n9wMyWmdkdY5EhEREpVqbEfSawa5vzISIiJU1plMDdrzGz2e3PSmKHHUa+ttdecPDB8M9/wtvfPvL9D384PB59FN73vpHvf/zjsPfe8NBD8MEPjnz/c5+DPfaAe++Fgw4a+f6XvwxvfSssWQKHHjry/WOPhTe9Ca6/Hr74xZHvn3ACbLkl/PrX8LWvjXz/e9+D17wGLr4Yjj9+5Ps/+hFstBGcey6ceurI93/2M1h/fTjzzPCod8klsO66cMop8NOfjnz/qqvC329+ExYuHP7eOuvApZeG/485Bq68cvj7M2bA+eeH/w8/HH73u+Hvz5oFZ58d/j/00HAM0zbZBObPD/8feCD8/vfD399yy3D8APbbDx5+ePj7b3wjfP3r4f/3vhdWrBj+/s47wxFHhP932w2eemr4+7vvDp//fPhf597I93Xuhf/Lnnu1/WmzltVxm9mBZtZvZv3Lly9v1WpllJ741yq+d/UfGFi6stNZEZEWMXdvnCiUuBe6+2ZlVtrb2+v9/f2jy5mM2sDSlex72mKeWTXEtCmTWHBAH3N7pnc6WyKSwcwG3L23TFr1KonMwNKVnLzo/lIl6MWDK3hm1RBDDs+uGmLx4IqGy4hI92tYxy3do2oJum/ODKZNmcSzq4aYOmUSfXNmjGFuRaRdGgZuMzsH2AFY38weBo5099PbnTEZKasEXRS45/ZMZ8EBfSweXEHfnBmqJhEZJ8r0Kpk3FhmRxpopQc/tma6ALTLOqKokIipBiwgocEdHJWgRUa8SEZHIKHCLiERGgVtEJDIK3CIikVHgFhGJjAK3iEhkFLhFRCKjwC0iEhkFbhGRyChwi4hERoFbpM2qzKEuUobmKhFpI92FSNpBJW7pejGXWCfqXYhi/sxioBK3dLXYS6wT8S5EsX9mMVDglq5W9a4/3WYizqEe+2cWAwVu6WrjocQ60eZQHw+fWbczd2/5Snt7e72/v7/l65WJaWDpyglVYh0P9JlVZ2YD7t5bJq1K3NL1JlqJdTzQZ9Ze6lUiIhIZBW4Zd9QVTbKMp/NCVSUyrqgrWpzaXSc+3s6LCVXiHk+/uN2i247pRB3wErNaUD3+8nvZ97TFpc+lKufeeDsvJkyJe7z94naDbjym6ooWn2b6fVc998bbeTFhArcGBbReNx7TiTjgJXbNBNWq5954Oy8mTOAeb7+43aBbj+lE7IoWc7/pZoJqM+feeDovJtQAnJhP7m41EY9pt+3zWFRZdds+Q3fmaTQ0ACdH1V/c8XZitMN4KsWU0Y31+u2usurGfYaJd+6lTaheJVU029It8Ym9d0Kt2mCy0ZYqq27c54luQpW4q+jGhrfxopuuZMZD74R2N7x14z5PdArcOXSytke3XXaPl94J7aw26NZ9brduKmDUU+DOMVYnazefHO3QbVcyE713QlkTbZ+7rYBRT4G7QLtP1m4/Odqh265kJmppUop1WwGjngJ3B3X7ydEO3RgoJ1ppUhrrtgJGva4K3BOt2qDbT452UaCUbteNBYy0UgNwzGxX4NvAZOA0d/9/RembGYAzEasNYOL9WEk5Oi8mnpYOwDGzycDJwC7Aw8BNZnaRu981umwONxGrDUClTxlpohZipLwyA3C2Be5390F3fwb4CfCuVmek3YMIRGKhAS/SSJk67pcDD6WePwy8odUZ6fY6JZGxMlHbPqS8ljVOmtmBwIEAG2+8cVPrULWBiAox0liZwP0nYKPU81nJa8O4+3xgPoTGyZbkTmSCUiFGipSp474JeLWZvcLMpgH7ABe1N1siIpKnYYnb3VeZ2SeBXxG6A/7A3e9se85ERCRTqTpud78EuKTNeRERkRI0H7eISGQUuEVEIqPALSISmbbcLNjMlgNLm1x8feDRCZR+LLbRbenHYhuxpx+LbcSefiy2MRb7UNPj7jNLpXT3rnoA/RMpfTfmSfvc+fTdmKduS9+NeWpmH5p5qKpERCQyCtwiIpHpxsA9f4KlH4ttdFv6sdhG7OnHYhuxpx+LbYzFPlTWlsZJERFpn24scYuISAEFbhGRyChwi4hEpmsCt5ltYGYb1x6tTt9kniaZ2QsK3j8k+btdE+ueZmabJY+po8nnaJjZ2mY2otO/mc00s7UbLDvdzLY1s3+vPdqQv644Tqn8NLXPjc6lVLr/KfNa6r1Xmtlayf87mNmnzexFBenNzPYzs68kzzc2s23L7EMZZvZ+M3t+8v+XzewCM9u6VetPbedNZvYBM/tQ7VGQttIxSi23npmtVyJdpc+sFTreOGlm7wSOB14GLAN6gLvd/XUtSr828BHgdcDqQOTu/ycn/Y+BjwHPEeYifwHwbXc/LiPtEnff0sxudvfSJ6eZ7QD8EPgjYIQbVfynu1+Tk34T4DDCvq6e0dHddyrYRh9wIrApMI0wJe8/3P0FdenmA5e5+wV1r78HeJu7fzxn/QcAhxBurLEE6AN+1yBP72Dk5/DVgvQ7UOI4mdnFQO6J7O7vHE361HKV9rnKuZRaZsS5ZGa3ufvrc9IvAXqB2YQZPH8BvM7d356T/lRgCNjJ3Tc1s+nA5e6+TV26Zo/Rbe7+ejN7M/A14DjgK+6eebvDJKC+N8l/+twuOi9+BLyS8Bk8t2YR/3RO+qrHaHPgLODFhPNuOeG8uyMnfaXPrBVaduuyUTiG8AX4tbtvZWY7Avu1MP2PgHuA/wC+CuwL3F2Q/t/c/XEz2xe4FPgvYIBwAta728zuA15mZrelXjfCiZT3wR1PCIr3wurAfA4wNyf9ecB3ge+z5kRt5CTCTS/OI5y0HwI2yUg3190PrH/R3S80s68VrP8QYBtgsbvvaGavBY7NS2xm3wXWBXYETgPeB9zYYB/KHqdvNlhPvVr6PYGXAmcnz+cBfy1YrtI+U+FcMrOPAwcDc+rOpecDvy3YxpCHOfPfA5zo7iea2S0F6d/g7lvX0rj7yuQGKfWqHtOa2vn5DmC+u/+ywXn0C+DvhOPydMlt9BKObdlSZ9Vj9D3gs+6+CFYXIOYDb0onGsVnNnpjMTyzzBBR4FZgUu3/Fqa/Jfl7W/J3KuGLl5f+ziTNecBbSqz/pUleeuofBcvcVua11HsDoziut6VeuyUj3d0F6yh676bk7xJgrdqxa7TPqb/rAdc22IdKx6nZY9TotVHsc+lzCXghoUR4Tt159OIG+3AD4QfnDuAVyWt3NEg/Gbg5eT4z67xIpd+j9j0reUwXEgLfIPAiYK0G35/cvBYscx6wYYX0VY/RiPzmvNbUZ9aKRzeUuB9L6pGuARaY2TLgHy1M/2xquc2AvwAbFKT/HuHS/FbgGjPrAR7PS+zufwG2KFhfln4zO401Jb19gf6C9Beb2cHAhaRKJe7+t4Jl/pmUpJaY2TeAR8hu01hmZtu6+7DSr5ltQ7hEzPNwUk/4c+AKM1tJ8cRiT6Xy9TJgBbBhQXooeZzM7HaKL+vzrnyeZ2Zz3H0wWc8rgOcV5KfqPpc+l9z974SS5zwzmwy8hHBFvJ6ZrefuD+ZsY39Cdcx/u/sDyT78qCBP3yGcRxuY2X8Trny+XJB+b+AEMzufcPerewrSAuwF7Ap8090fM7MNCdV8ea43s83d/fYG601bH7jLzG5k+Pchs/qG6sdo0MyOSKXZj/BDNMwoPrNR64Y67ucB/yJUL+xL+BVb4O4rCtI/RQhCtfRn5wWxpF7yfOD1wBmEkt5X3P27FfI4xd1XZbz+U3ffKyNwFFaVJPV6nwDenLx0LXCKu2deKprZAxkvu7vPKchzD6ENYCrwGcJxOsXd769Lty3wU+BMwuUqrKla2cfdb8jbRmodb0nWf5m7P5OT5ghCnfvOwMmE43Waux9RsN5SxynZ11zunhlczWxXwiXwIOEz6wEOcvdfFa0vWbbhPucsl3kupd7/JHAUocpmaM0utK6+NKni2Zmwz1e6e1HVIRYaVecRAqATvkfnuPsTGWlfCTzs7k8nVQyvB85y98fq0tW+M1OAVxM+g6dpXM1YO/YjuPvVRfuRLDsd2Mjdb2uQ5miGn3dHufvKnPRt/8xGbLPTgbsmOTnSjRN5gXg3d7+07rWPVQnEJfJSqhHNzDZ090fyAkdewOg2ZvYSQl3dZslLdwInufuyjLQvLlpXg6uA2jrWAtZOSiwdleTltcnTe/J+PJO0xxCu9K5396KrvPQyVRtk7yfUQ2cWXFLp8goNtW0UBb506bCWvrB0aGYzgA8ChxLaiF4FfMfdT6xLV6ohsNkf22aY2VXAOwn7O0Ao0PzW3T/bovWX+sxaqeNVJWZ2EOHX7V+EXysjnIh5pckjzOxpd/9NsvxhwE6Exrus9Vdqta7SiObujyR/S51ko7ikx8zexMh9OKsV23D3vwJHWtIt0N2LqkgGkvUbsDGwMvn/RcCDwCvq8rOTu//GzPbMyCte15ul7v3tCCWZ+t40c+rSPUH2PtdKb0Xd8Oay5rhukeRpxHFNDBJKnt9JtnktcI27/yIn/800yD5EuPxu5JDk7+4l0qbz9CngSELp8DnWfN/yrg7fBXyYEKjPArZ192Vmti5wF+EqKq3WELgnBQ2Bte+Mhd5Pd9ZK70kBblMKqqCSdf8PocrTaPw5v9BDI/EBhNL/kXWNibX1XpS3zSTPeVUxZT+zlul4idtCr4w3unupycfNbH1CA8hhhLq01wLzCi7RL2NNq/XqHhnufnxO+lp3ptrf9YBL3X37jLSVAsYoLulLd3+qug0zM8IX+ROERiuSbZzYoGT4feBCDzeSxsx2A97t7gfVpTs6+aKckZ2d7G6ZybL3EKp56j+7lpRsqhzXuuVeSqjL/Tww3d2fn5Ou9LmUWuZ04DXALxlef/utCrtWlPdKpUMzO5NQtz2iq6qZ7ezuV9a9dgNwAvAlYI+kTvkOd9+sfvkk/S3A1p4EIjObRGggzu1em+zDHo2qeFLpbwfeRuha+iV3v8kyuutZuAHMQ4TGxhsI3+PV8qpi2v2ZZel4iRv4A/DPsond/VELfbl/TfhCv6/2oeeY5e67VshP6Ua0vC9snnTQTKonan1nb8yqlkgp3f2pbhs9wKvd/ddmtg7Zn/dngO0IJakHkuXmAKea2Wfc/X9zNtXn7h9NbfdSC42g9fk5Mvm7f6O8Z/h7fbVYi1XqVmahofTfCKXVawkl6JsLFmmmQfbB5DEteTTKU7rwMI3QpjGiv35K6dJhUqXSkxW0AeqDdqJqQ6Clj7+7D5lZo7j017JBO/FV4FfAdUnQngPcl5HupcAuhKuqDxAC8TnufmeD9Vf6zFqhGwL34YSW5RsY/ms1rNSTOkFrl3bTCNUp7zOzosukqq3WCy30HDiO8KV0wmVuQ2a2AcPrMjPrDc1sr2T9VxH250QzO8zdf5az6jsIJ9UjJfcBM/socCBhEMErCYNGvktolEr7ILBL+orH3QfNbD/gciAvcP/ZzL7M8B4ffy7IT1Z94t8JXR2X1KWtlbYWmdlxwAUMPzeKgmUVVY/rDMJVyWPA34BHvaChkSbOJXc/GsDM1nX3hgWadOEhuXp6F2GcQ55B4Coza1g6dPfnzGzIzF5Ytj3C3e8CPp16/gChWiM3P2b2aeDU5PnBZPTgqNNvZucSevek9yGz2s3dzyN0Iaw9HyRUn9anew64DLgsqWKdRzhWR7v7SXmZSX1m6yXPn2yQ/1HrhqqSG4HrgNtZ0yKLu/+wReu/i1A/9wAlW61Ty5ZqRLPqozlvJQTLZcnzmYQBRZndCs1sEbAloX60TPenWiPRtsAN7r5V8trt7r55Xbqiy9ii915MqGKpDfm+Bjja8xuVf0wo4V6cvLQ7cBuhfvk8d/9GKu2ivP0ifHa5ozOraOa4JsttShjQ9RlgsrvPKrGtsufSG4HTgfXcfWMz24LQ0+XgRttIreOW2mee8d6RWa/Xgk9G+l8AWwFXkOp2m1edZGavBr5OuDJJF2Iy26ySws53CO1UDlwJHFp0BVq22s3MvuDu3zCzE8luwM2qalyLMHhoHuHcvIhQVfSngvxsRriqqDXcPwp8qERJvWndUOKe6hVad5MGqyXu/o+kVLg1cEJe6RbYrWqGrK4hsEGDFVQfzTmp7sRcQfG8MUeVyXedp939mVAIg+TyM+tXuqgrW+57SYA+JO/9DLMIdZlPJvk5knAp+u+EKq/Vgdvdd6yw3tE4qkpiM9sd2J6Q5xcBvyFUmeSlfz+hu+AThDaZrc3sGHcvGrV3AuFH4SIAd7/VCuZDseGNvpMIP47/ykufF6ALXJA8yjqD8IP+v4RG2f0pOLeT78E+VTJUodqtVp1SNEZiNTM7i9Cz6hJCISRziHuG+Ywcafl96kZatlI3lLiPJQxSuJgSg0uS1uAtCK3gZxIuPfdy98y+nckyWxC+cBBG691akLZyg5WZ9bt7b1KS3iqpp7u1oAR9XJL/c5KX9iaMCPy/BduoUidOUt/8GKE/9qcIl6B3ufuX6tI9R/YAJiOUEKfWpT/B3Q+1nLks8kqrSUPj5u7+bPK8NqLutQ1KiJW607WTmZ1ECNTXuntutVAqfaV5O5JlbnD3N6SPSYNzKV36XEX4Ln2//vxo9nOryswG3H1u+uqu9lpdumZKw5WXqZj3IdZ8F7LGZWRWx2Z9PkWfWSt0Q4l7XvL38NRrRd0BV7m7W+imdJK7n25mH8lbuYUZ/D7KmlLD2WY23+v6n6ZUnQcBSo7mNLNXAS9x98OSklKtg//vgAUF+1C1ThzCvBgfIVRBHUQoRYyoX3X3yfWvNVBraKo6l8UC4Ibk0hvCUOofWxhQdVfWAtZcd7qGrMnug+7+ydoPaFIP3+gHtOq8HQAPJVd8bmE2xEMomFunQumzqc+tatUH8LSFniH3WRiY8ifCoLd6lUrDzSxjFbv3uXuzs6WWGmnZUt7mMfWtfgBXE4L87wkNS5OA2wvS3wY8L/X8eRTPC1J6HgRC3fl2yTonEX4I/xP4CmHypvr0CwmlzvrXNwcuLtjOrcAGqeczKZj/oS7dzDZ8BpMJo1urLtdLCESHAL0l0lee36TN5977Cf2Lf0jo0/wAoVdTXvpK83Yky6xP+JH7K6G95GwK5r4gVEFdmKRdRhglPKuF+3wdoUH7NkLbzVHAVwvSb5N8TrMI1SYXEHogdeLzWk5oFD6MUL31lvSjhduZTqinvzl5fJvQTbRt+9bxqhJYXblf/4ueWadsoQ/tBwgT/lxrYS7uHQrS3w5s4+7/Sp6vnSy7eU760g1WZrYQONzreqxYmBbyWHffo+71m7xu+sx0PgvyNOy9pERza1b6pGfBkcAnWVO32LBfdlVmdh1hatCGw70tdCu7091f2yht3XK1aoPFhJn8ViTreVVTmR6lJhqV1yWMNbjd3e+zMG/H5u5+eUbajdz9oZz17O7uC3PeuwL4McNLe/u6+y456as2Hpaq+qiqamm4mWWS867Wve/1lO/e1/U6XlWSNFLtQDiRLiE0Jl5HKNGM4GFSp2+lnj+YlzZxBuES/cLk+buBHxSkP6pk1iFUe4zoZujut5vZ7Iz0RZO3r1Pw3mVm9iuG14lfkpO21i97G6/WL7uqQeC3yZcp3dsgr1vZvWa2sVebeKfprpltUqlR2d3/mVQNvcTW3Owjb5KmK8xsV3f/Y/pFM9ufMAlUZuAmXFGl67nPNLNDC/ahUuMhJas+mgjEb6RgsEuOSst4k937ymrmx6dVOl7iTkrEWxCmltwiqUM8u6DEkHWDgCfd/YUF29ia1IQxXtyqXyXv97n7q3Peu7++ZGhm5wC/cffv171+AKEkt3fd67U68d/W1Yk/Rqiq+EPGdm+hrl928vpMwoT5mY2AVVn1bmXXELqV3cjwQF/q5LYumN+kaqOyDR9eXjj5kJm9ndCj5B3ufl/y2uGEq8vd3P3hnG1cSTLpU/LSPGB/d6/vr19LX6kEbWGWyLsJhY5jCBNrfcPdF9elqzTqsJnScJPLVO7eV1bVfW6pdtbDlHkQGnggdAl7QbLz9xSk7yfULd9CCNr7A18vSH86sGXda0cVpO8j3K3kSUJ3uOeAx3PSngN8NOP1A4BzM15/CXA9oZHx+ORxNaFx8qUZ6SvXiVM8z3DluY9LfH7rEfocN0r3lqxHTto9ix4dPl/fS7ji+xbwngZp7wdmVFj3zskymxGC+PU0qCsl1DtfRKjPXUYYlLJRQfrrCSXsCwjVae8B7m3BcZlMqBb6YfLd/Bphcqkyy65FmA9lOfDJVi1DuBK/OcnLZm04F5re51Fveyw20mDnTyH8mn+MMAz1FuCMgvSlbhCQeu9hQuPeh1Kv3Vy0fkr+MFQNxKnldiR00fsUoZ44L91NBe9lNsg22Lfc95r43DZLjtHS5DHQ6KRNjtfuyWODgnRnJI9fEiaxOj95/A1Y2OlztsIxWgRMqbjM9oQBHBcRrjCa2e6hBe+VajxMtp/7aLD9UoEuaIDIAAAGK0lEQVQ4SbcnoUPATcARwMtLrLvUMoSrnCeSx+OpxxPkFMZG8VlX/vEZzaPjVSVpSb3wC7x4rtxrgLcS6jr/Qhiu/GHPbyC6mRAozybMJ3AIISDm9Ruu9clePQmNFfQzTt7fkdSUqJ7MXDhaVatiktcr9cseRd6uJ0zYkx50cKy7Zw46yOjSuD1Q2KXRzC4n3OvvkeT5hsCZ7v4frdiHsprtPmgVJh+y4VM6rEW4Acjq2fvytpGz3QfdfVQ30G6mGqBKtYQNH+zyEy8x2KWZZdqtnVUxhdvthsCdqr91wkQwFxak7aHEDQJS6dMDGY4iBP2XeX4reqUfhnaqWic+lqzioIOqPTKSNHe7+6ap55MIP4yb5i3TTaq2A7Rwuw+5+0Z1r7W1R0bVoGpNDHZpZpl26uQPSccDt5mdQqiaSDf4/MHdP9Gi9R/tyQx1yfM9gM94/p25ewiNSdMo8cPQTklD7YWEuvb03WmmEepX/zLWeUrl7UJC/WG6G9pcd39PTvrSXRpTaU4i3B0lfW7c7+6fasEujFtZJe7RNKSlemQcRxgKPqJHRrcF1bHQyX3uhsB9D7CpJxnJK1VZ8Q0CvFUlYjOb6+4Dda/l9qMdC+2qihkNq357p6weGbe7+xcabGdP1kxXcE3R1Vi3Sa4qvsDIIfujniSrQfXNOu4+pS59V/XIkNHphsC9EPiEr7kjRg9hKHv94JWerMWBjQiDYOpvjXSdu7854wRvVC95M6Eh847k+T6EEnru/BJSTl2XxmtjCsLNSOrozyXccOFjhFG1y71gTpoxyleZEnTX1SfLGh0L3LZmspsXElq6b0yev4HQRXCHgmW3IvRvfT9h2PH5WSdfk/maA/wsWf/2hEmadvcuuD9iN7CcSYpq6utK65btYc2NHdYlTImadcPZ0dyKrGuk+kynG7pzR8+OQX6qNB5OuKqPmHRy5GTVyW42IZxw8wjdpc4l/PDkTgFqTQy19nATgX0I/WEfBN7m7k81WGwiqTq5FAA28sYOLyf7xg54xTsLdbFnk7+PWJjl8M+smbN5TFnFKUu9+QmXZAx0vKqkrKQEcC3wkVpDoZkN5vUOSS33C+BT3vgu1vV16BsQ7tDyNICXuPHCRGDVh63Xlit1Y4fxxML83dcSqvNOJAwwO8rdLy5csD15UQl6HOmGuUqyhrBn3TNvT8KE64ss3AD4J1BqfoPpwJ0W7rSz+sR193fVpat0t+wJ7OeEm1dgZue7+4hbQOUoe2OH6FkyYVSqQfvvhLEEtWA+5lSCHl86HriBkwgB+TxCV7cPAZvUJ3L3nwM/tzB/87uAQ4ENzOxUwt3GR8y4ljgi9X9t4MeIO264+9JmqlYmoPSPZeHVTp2rzeyLwDpmtgvhxg5jXvIcI81OGCVSSlf8CidVH5Pd/TkPM53l3pXd3f/h7j9Oep3MIgy7zm2lT/qnPk4oUZ9JuLfdd3PSPgfca2tmcpORPOf/Rv6LMBz4dkJd9y+97m4848hngcstTKEKUJsw6rOEOVpERqXjddztGqmY05j5eXfP6lZYn5/6WeyyqlYmpNSQeiNMRVu7E3lmXamFOxXNcveTk+c3Em7w4MAXvPguPtEys50JN1F4N2HSsW0JM/9l9nMXqaIbAnelIewV1ttsY2a6RLS6asVz7tguxczst4Tj91DyfAnhqmc9wmRimdOPjgdmtj1h5Ov1hPui5t7EV6SKjtdx1wbeAE8RRuK1SlONme5+dUY/8cyqFSllmg+/s8t1Hm4E/bekvWLcsZETRu0MLLPQMqseHDJqnRyAM1ZD2GuNmfMIJb2zyGjMbLZqRYpZziyGyXt/cPdXjnWeRGLXycBdaQh7i7Y5nVCK3rv+Er3ZqhUpZmYLgKt85AyHBxHuFTqvMzkTiVfH67ihvUPYK+Th3YSqle0I96n7CXCau79iLPMx3pjZBoS+308TZhMEmEuoQni3u/+1U3kTiVUnS9xdWTVRtmpFqjGznQiz5EGXzHAoEqtOBu6ur5ooqloREemUTgZuVU2IiDSh43XcqpoQEamm44E7TVUTIiKNdVXgFhGRxrpikikRESlPgVtEJDIK3CIikVHgFhGJjAK3iEhk/j9VfMq0bnCCpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(a, '.')\n",
    "plt.plot(np.repeat(1, len(xc.index)), 'r--')\n",
    "plt.xticks(rotation=90)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### カリフォルニア州のデータをサンプルにS/N比を計算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "xc_prime = xc.loc['Calif']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "sn = 10 * np.log10(xc_prime**2/np.diag(sx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADrRJREFUeJzt3X+sZGV9x/H3RxYFKyrIFhFcL1KwYhuJ3tCildC4WC0xiIUU66+2tlsjpNaUJhjahMYaKVWpadrC2pBF608gCwgqCLXSIip3C8IC8kO4FCiVi8ZWKyEBvv1jzsqw3t17986cO3f3eb+SyTznnGfO8zwzcz73zJkz96SqkCTt+p426Q5IkpaHgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqxKpJd2DYvvvuW1NTU5PuhiTtVDZt2vRwVa1eqN6KCvypqSlmZmYm3Q1J2qkkuXcx9TykI0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWrEivrhlSTtiKnTLp90F8Zm9sxje2/DPXxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiPGEvhJzkvyUJLNQ/P2SfKVJHd293uPoy1J0tKMaw9/A/D6readBlxdVYcAV3fTkqQJGUvgV9U1wA+2mn0ccH5XPh940zjakiQtTZ/H8Perqge78n8D+81XKcm6JDNJZubm5nrsjiS1bVm+tK2qAmoby9ZX1XRVTa9evXo5uiNJTeoz8L+XZH+A7v6hHtuSJC2gz8C/FHhnV34ncEmPbUmSFjCu0zI/A1wHvCTJ/UneBZwJHJPkTmBtNy1JmpBV41hJVb1lG4teO471S5JG5y9tJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGrOq7gSSzwI+Ax4HHqmq67zYlST+r98Dv/HpVPbxMbUmS5uEhHUlqxHIEfgFXJtmUZN3WC5OsSzKTZGZubm4ZuiNJbVqOwP+1qnoF8Abg5CRHDS+sqvVVNV1V06tXr16G7khSm3oP/Kp6oLt/CNgIHNF3m5Kkn9Vr4Cf5uSR7bSkDrwM299mmJGl+fZ+lsx+wMcmWtj5dVV/uuU1J0jx6Dfyquht4eZ9tSJIWx9MyJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNaL3wE/y+iS3J7kryWl9tydJml+vgZ9kN+DvgTcAhwFvSXJYn21Kkua3quf1HwHcVVV3AyT5LHAccGvP7UrNmDrt8kl3YSxmzzx20l3Y5fV9SOcA4L6h6fu7eZKkZdb3Hv6CkqwD1gGsWbNmpHW1vKezq4wddnz8LY99qY/ZVbQ89qXoO/AfAF44NH1gN++nqmo9sB5genq6eu6PdkFu9NLi9H1I53rgkCQHJXk6cBJwac9tSpLm0eseflU9luQU4ApgN+C8qrqlzzYlSfPr/Rh+VX0R+GLf7UiSts9f2kpSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9Jjegt8JOckeSBJDd2t9/sqy1J0sJW9bz+s6vqwz23IUlaBA/pSFIj+g78U5LclOS8JHv33JYkaTtGCvwkVyXZPM/tOOAfgYOBw4EHgY9sYx3rkswkmZmbmxulO5Kk7RjpGH5VrV1MvSQfBy7bxjrWA+sBpqena5T+SJK2rc+zdPYfmjwe2NxXW5KkhfV5ls5ZSQ4HCpgF/qjHtiRJC+gt8Kvq7X2tW5K04zwtU5IaYeBLUiMMfElqhIEvSY3o+3/paJnMnnnspLsgaYVzD1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDVipMBPcmKSW5I8kWR6q2XvT3JXktuT/MZo3ZQkjWrViI/fDLwZOHd4ZpLDgJOAlwEvAK5KcmhVPT5ie5KkJRppD7+qbquq2+dZdBzw2ap6tKruAe4CjhilLUnSaPo6hn8AcN/Q9P3dPEnShCx4SCfJVcDz51l0elVdMmoHkqwD1gGsWbNm1NVJkrZhwcCvqrVLWO8DwAuHpg/s5s23/vXAeoDp6elaQluSpEXo65DOpcBJSZ6R5CDgEOBbPbUlSVqEUU/LPD7J/cCRwOVJrgCoqluAzwO3Al8GTvYMHUmarJFOy6yqjcDGbSz7IPDBUdYvSRoff2krSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDVipMBPcmKSW5I8kWR6aP5UkkeS3Njdzhm9q5KkUawa8fGbgTcD586z7LtVdfiI65ckjclIgV9VtwEkGU9vRjR75rGT7oIkrVh9HsM/KMkNSb6W5DXbqpRkXZKZJDNzc3M9dkeS2rbgHn6Sq4Dnz7Po9Kq6ZBsPexBYU1XfT/JK4OIkL6uq/926YlWtB9YDTE9P1+K7LknaEQsGflWt3dGVVtWjwKNdeVOS7wKHAjM73ENJ0lj0ckgnyeoku3XlFwOHAHf30ZYkaXFGPS3z+CT3A0cClye5olt0FHBTkhuBC4F3V9UPRuuqJGkUo56lsxHYOM/8i4CLRlm3JGm8/KWtJDXCwJekRhj4ktSIVK2cU9+TzAH3TrofC9gXeHjSnZiQlscObY+/5bHDyh//i6pq9UKVVlTg7wySzFTV9MI1dz0tjx3aHn/LY4ddZ/we0pGkRhj4ktQIA3/HrZ90Byao5bFD2+Nveeywi4zfY/iS1Aj38CWpEU0HfpIzkpy6hMcdneRVQ9Mbkpww3t7tHJL8eNJ9WEm698Zlk+7H9iR5bpL3TLofK0WSP05yW5JPLeGxs0n27aNffWg68EdwNPCqhSpp55KBBbeJJKNeGnTSngsY+E96D3BMVb110h3pW3OBn+T0JHck+XfgJd28g5N8OcmmJP+W5Be7+W9M8s3uyl1XJdkvyRTwbuB93QXat1zN66gkX09y95a9/ST7J7mmq7d5e1f+6lt3YfnvJPlUtzdzYZJnJnltN76bk5yX5Bld/dkkZ3Xzv5XkF7r5ByW5rpv/V1u18WdJrk9yU5K/HGr3tiQf7y54f2WSPZf/GZhf17/bk3yCwTWaHx9adkKSDV15Q5JzknwTOCvJEd3zcEP3ur9kMiNYkjOBg7v35d9s53X7TjfuO7r3zdok1ya5M8kRXb0zknyyey7uTPKHEx3ZDkpyDvBi4EtJ/mf4E3+3zU515bd128GNSc5N9+/fdzpV1cwNeCVwM/BM4NnAXcCpwNXAIV2dXwH+pSvvzZNfbP8B8JGufAZw6tB6NwAXMPgDehhwVzf/TxlcGQxgN2CvCY59Cijg1d30ecCfA/cBh3bzPgH8SVeeHer7O4DLuvKlwDu68snAj7vy6xicyZDuebiMwb/JngIeAw7v6n0eeNuk3wtbPS9PAL/aTf94aNkJwIah1/gyYLdu+tnAqq68FrioKx+95blaqbduzJsX+br9cjd/U/eeCXAccPHQtvBtYE8Gv0a9D3jBpMe4g8/HbNf3rbfrzd3z8FLgC8Du3fx/GNoGZoF9Jz2Gxd529o+mO+o1wMaq+glAkkuBPRgcnrkgT16M/Rnd/YHA55LsDzwduGc76764qp4Abk2yXzfveuC8JLt3y28c62h23H1VdW1X/mfgL4B7quqObt75DEL8b7vpzwzdn92VXw38Vlf+JPDXXfl13e2GbvpZDC58859dG1vGvonBRrSS3FtV31hEvQuqassngOcA5yc5hMEf0t17612/FnrdbgZIcgtwdVVVkpt56mt4SVU9AjyS5KvAEcDFy9T/5fBaBjuL13cZsSfw0ER7tEStBf58ngb8sKoOn2fZ3wEfrapLkxzNYA9gWx4dKgegqq5JchRwLLAhyUer6hPj6faSbH0O7g+B5y2y/rbKWwT4UFWd+5SZg4/Ew8/N4ww2mJXk/4bKw2PbYzv1PgB8taqO78b4r730rH+Lfd2eGJp+gqdmx9bvh531XO/HeOph7i2vf4Dzq+r9y9+l8WrtGP41wJuS7JlkL+CNwE+Ae5KcCD/94u7lXf3nAA905XcOredHwF4LNZbkRcD3qurjwD8BrxjPMJZsTZIju/LvMLjG8NSW4/PA24GvDdX/7aH767rytcBJXXn4S64rgN9P8iyAJAck+fkx9385fC/JSzP48vb47dQbfm/8bu+9Gq/h9+84XrfjkuyR5HkMDmldP7aeLq9Zum00ySuAg7r5VwMnbHlekuzTbds7naYCv6r+A/gcg2OOX+LJN+ZbgXcl+TZwC4NjlDDYo78gySae+p/yvgAcv9WXtvM5Gvh2khsYhObHxjSUpbodODnJbQy+nzgb+D0GY7yZwZ7bOUP1905yE/Be4H3dvPd267gZOGBLxaq6Evg0cF237EIW8UdxBTqNwXHsrwMPbqfeWcCHutd2p/qkXFXfB65Nshk4htFft5uArwLfAD5QVf81zv4uo4uAfbrDV6cAdwBU1a0Mvu+6stsevgLsP7FejsBf2jai+4h+WVX90iLrzwLTVbWS/yWsJizJGQy+6P7wpPuihTW1hy9JLXMPX5Ia4R6+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJasT/A7wmRgC5RLZpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(sn.index, sn)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ワシントンDCの場合"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADflJREFUeJzt3X2MZXV9x/H3R1YRKyrKShFdh1q00jYSnWDVamhEa0sM0mKK9amtdtsoqTW1CcY2obHGLfWhpmmrqyEL1kc0KIIPCLWlRVR2C8Ii8hBZCtTK2sZGq7EBvv3jntHLdnZndmfO3N39vl/JzfzOOb85v+/vPnzm3HPP3U1VIUk6+D1g1gVIktaGgS9JTRj4ktSEgS9JTRj4ktSEgS9JTRj4ktSEgS9JTRj4ktTEulkXMO3II4+subm5WZchSQeUbdu2fbuq1i/Vb78K/Lm5ObZu3TrrMiTpgJLk9uX085SOJDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSE/vVF6+07+bOumTWJayaHZtOmXUJ0kHJI3xJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmViXwk5yb5O4k26fWPTLJ55PcMvw8YjXGkiTtm9U6wt8CvGCXdWcBl1fVccDlw7IkaUZWJfCr6grgv3ZZfSpw3tA+D3jRaowlSdo3Y57DP6qqvjm0/wM4asSxJElLWJMPbauqgFpsW5KNSbYm2bpz5861KEeSWhoz8L+V5GiA4efdi3Wqqs1VNV9V8+vXrx+xHEnqbczAvwh45dB+JfDJEceSJC1htS7L/BBwFfCkJHcmeRWwCXhekluAk4dlSdKMrFuNnVTVS3az6bmrsX9J0sr5TVtJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJamLd2AMk2QF8F7gXuKeq5sceU5L0/40e+INfqqpvr9FYkqRFeEpHkppYi8Av4NIk25JsXIPxJEmLWItTOr9YVXcleTTw+SRfr6orFjYOfwQ2AmzYsGENypGknkY/wq+qu4afdwMXAifusn1zVc1X1fz69evHLkeS2ho18JP8RJLDF9rA84HtY44pSVrc2Kd0jgIuTLIw1ger6rMjjylJWsSogV9V3wCeMuYYkqTl8bJMSWrCwJekJgx8SWrCwJekJgx8SWpirf7xNEkjmTvrklmXsCp2bDpl1iUc9DzCl6QmPMLXAe9gOcIFj3I1Lo/wJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJakJA1+SmjDwJamJdWMPkOQFwLuAQ4D3VdWmMcaZO+uSMXY7Ezs2nTLrEiQdhEY9wk9yCPA3wK8AxwMvSXL8mGNKkhY39imdE4Fbq+obVfW/wIeBU0ceU5K0iLED/xjgjqnlO4d1kqQ1lqoab+fJ6cALqurVw/LLgadX1ZlTfTYCGwE2bNjwtNtvv320eiQdXPzsbiLJtqqaX6rf2Ef4dwGPm1p+7LDuR6pqc1XNV9X8+vXrRy5HkvoaO/CvBo5LcmySBwFnABeNPKYkaRGjXpZZVfckORP4HJPLMs+tqhvGHFOStLjRr8Ovqk8Dnx57HEnSnvlNW0lqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqwsCXpCYMfElqYrTAT3J2kruSXDvcfnWssSRJS1s38v7fWVVvG3kMSdIyeEpHkpoYO/DPTHJdknOTHDHyWJKkPVhR4Ce5LMn2RW6nAn8HPAE4Afgm8Pbd7GNjkq1Jtu7cuXMl5UiS9mBF5/Cr6uTl9EvyXuDi3exjM7AZYH5+vlZSjyRp98a8SufoqcXTgO1jjSVJWtqYV+mck+QEoIAdwO+NOJYkaQmjBX5VvXysfUuS9p6XZUpSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEwa+JDVh4EtSEysK/CQvTnJDkvuSzO+y7Y1Jbk1yU5JfXlmZkqSVWrfC398O/BrwnumVSY4HzgB+FngMcFmSJ1bVvSscT5K0j1Z0hF9VN1bVTYtsOhX4cFX9sKpuA24FTlzJWJKklRnrHP4xwB1Ty3cO6yRJM7LkKZ0klwE/ucimN1XVJ1daQJKNwEaADRs2rHR3kqTdWDLwq+rkfdjvXcDjppYfO6xbbP+bgc0A8/PztQ9jSZKWYaxTOhcBZyQ5NMmxwHHAV0YaS5K0DCu9LPO0JHcCzwAuSfI5gKq6Afgo8DXgs8BrvUJHkmZrRZdlVtWFwIW72fYW4C0r2b8kafX4TVtJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmDHxJasLAl6QmVvRfHErSLO3YdMqsSzigeIQvSU0Y+JLUhIEvSU0Y+JLUhIEvSU0Y+JLUhIEvSU0Y+JLUhIEvSU2kqmZdw48k2QncPus6lnAk8O1ZFzEjnecOvefv3Pdvj6+q9Ut12q8C/0CQZGtVzc+6jlnoPHfoPX/nfnDM3VM6ktSEgS9JTRj4e2/zrAuYoc5zh97zd+4HAc/hS1ITHuFLUhOtAz/J2UnesA+/d1KSZ04tb0ly+upWt/9L8r1Z17C/GZ4bF8+6jj1J8ogkr5l1HfuLJH+Q5MYkH9iH392R5Mgx6hpD68BfgZOAZy7VSQeWTCz5mkhyoP9PcY8ADPwfew3wvKp66awLGVu7wE/ypiQ3J/kX4EnDuick+WySbUn+OcnPDOtfmOTLSa5JclmSo5LMAb8PvD7JtUmePez6OUm+mOQbC0f7SY5OcsXQb/tU3zWVZC7J15N8YDiS+ViShyR57jC365Ocm+TQof+OJOcM67+S5KeH9ccmuWpY/+e7jPHHSa5Ocl2SP5sa98Yk701yQ5JLkxy29vfA7g013pTkfGA7cO/UttOTbBnaW5K8O8mXgXOSnDjcF9cMj/uTZjODfbIJeMLwvPzLPTx2Xx/mffPw3Dk5yZVJbkly4tDv7CTvH+6LW5L87kxntpeSvBv4KeAzSf57+h3/8JqdG9ovG14L1yZ5T5JDZlPxClVVmxvwNOB64CHAw4BbgTcAlwPHDX2eDvzD0D6CH3+w/Wrg7UP7bOANU/vdAlzA5A/o8cCtw/o/At40tA8BDp/RvOeAAp41LJ8L/AlwB/DEYd35wB8O7R1Tdb8CuHhoXwS8Ymi/Fvje0H4+kysZMtwHFwPPGca9Bzhh6PdR4GWzfh4sct/cB/zCsPy9qW2nA1umHuOLgUOG5YcB64b2ycDHh/ZJC/fX/nob5rx9mY/dzw/rtw3PmwCnAp+Yei18FTiMyTdS7wAeM+s57uX9sWOofdfX9fbhfngy8CnggcP6v516HewAjpz1HJZ7O9Dfmu6tZwMXVtX3AZJcBDyYyemZC5Is9Dt0+PlY4CNJjgYeBNy2h31/oqruA76W5Khh3dXAuUkeOGy/dlVns3fuqKorh/bfA38K3FZVNw/rzmMS4n81LH9o6uc7h/azgF8f2u8H/mJoP3+4XTMsPxQ4Dvi3YYyFeW9j8gLa39xeVV9aRr8LqmrhHcDDgfOSHMfkj+kDR6tuXEs9dtcDJLkBuLyqKsn13P9x/GRV/QD4QZIvACcCn1ij+tfCc5kcLF49ZMRhwN0zrWgfdQv8xTwA+E5VnbDItr8G3lFVFyU5ickRwO78cKodgKq6IslzgFOALUneUVXnr07Ze23X62+/Azxqmf13114Q4K1V9Z77rZy8HZ6+X+5l8mLZ3/zPVHt6fg/eQ783A1+oqtOGef7jKJWNb7mP3X1Ty/dx/+zY9TlxoF7rfQ/3P8298PgHOK+q3rj2Ja2ubufwrwBelOSwJIcDLwS+D9yW5MXwow/unjL0fzhw19B+5dR+vgscvtRgSR4PfKuq3gu8D3jq6kxjn2xI8oyh/ZvAVmBu4fw88HLgn6b6/8bUz6uG9pXAGUN7+gOuzwG/k+ShAEmOSfLoVa5/rXwryZMz+fD2tD30m35u/NboVa2u6efvajx2pyZ5cJJHMTmldfWqVbq2djC8RpM8FTh2WH85cPrC/ZLkkcNr+4DTKvCr6l+BjzA55/gZfvzEfCnwqiRfBW5gco4SJkf0FyTZxv3/tbxPAaft8qHtYk4CvprkGibB+a5Vmsq+uAl4bZIbmXw28U7gt5nM73omR23vnup/RJLrgNcBrx/WvW7Yx/XAMQsdq+pS4IPAVcO2j7GMP4j7qbOYnMf+IvDNPfQ7B3jr8NgeUO+Uq+o/gSuTbAeex8ofu+uALwBfAt5cVf++mvWuoY8DjxxOX50J3AxQVV9j8pnXpcNr4vPA0TOrcgX8pm0Dw9vzi6vq55bZfwcwX1X7+z8JqxlLcjaTD7rfNutatLRWR/iS1JlH+JLUhEf4ktSEgS9JTRj4ktSEgS9JTRj4ktSEgS9JTfwfKv0QJ/H8cRsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xc_prime = xc.loc['DC']\n",
    "sn = 10 * np.log10(xc_prime**2/np.diag(sx))\n",
    "plt.bar(sn.index, sn)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}