Skip to content

Instantly share code, notes, and snippets.

@shotahorii

shotahorii/Deviance_and_AIC.ipynb

Created March 12, 2019 21:36
Show Gist options
  • Save shotahorii/99342b1b4d63b2e743301b7ab8e56f19 to your computer and use it in GitHub Desktop.
Save shotahorii/99342b1b4d63b2e743301b7ab8e56f19 to your computer and use it in GitHub Desktop.
Download ZIP
Deviance and AIC
Display the source blob
Display the rendered blob
Raw
Deviance_and_AIC.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Deviance and AIC\n",
"- 「データ解析のための統計モデリング入門」(久保) 4章\n",
"- [モデル選択と統計学的検定](http://hosho.ees.hokudai.ac.jp/~kubo/stat/2012/kobe/k3/kubostat2012k3.pdf)\n",
"- [AIC(赤池情報量基準)を使った線形回帰のモデル評価](https://www.iwanttobeacat.com/entry/2018/02/24/123626)\n",
"\n",
"[MLE and least squares](https://gist.github.com/shotahorii/be3bda5e42eb11fe7b721254768408ab), [GLM](https://gist.github.com/shotahorii/cc2508accc0b2e6815663b0531b1edb9), [Poisson Regression](https://gist.github.com/shotahorii/c4f793a200f6b56683e383eae4352ea8), [Logistic Regression](https://gist.github.com/shotahorii/63b10aa9337a4e3893ed3dac28af84e5)などの記事で、最尤法によるパラメータ推定について見てきた。 \n",
"これらは、尤度(実際の計算上は対数尤度)を最大化するパラメータの値を探す、つまり与えられた観測データへのあてはまりが最も良いパラメータの値を探す手法であった。では、複数のモデルがあった時に、この最大対数尤度の最も大きいモデルが最も良いモデルなのだろうか？ \n",
" \n",
"例えば[ここ](https://gist.github.com/shotahorii/c6033a4cc8de5b799b23d0bd386b1370)では、最小二乗法(=等分散正規分布を誤差分布に仮定した最尤推定)でモデルを複雑化していった際の過学習問題について見ていった。この例でいえば、より複雑なモデル(過学習しているモデル)ほど対数尤度は大きい。実際に見てみよう。"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import math\n",
"import numpy as np\n",
"from functools import reduce\n",
"\n",
"from matplotlib import pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 229,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# generate random dataset\n",
"x = np.linspace(0, 1, 100)\n",
"y = 0.5 + 0.4*np.sin(2*np.pi*x) + np.random.normal(0.0, 0.1, len(x))"
]
},
{
"cell_type": "code",
"execution_count": 230,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# create matrix of X \n",
"# for example, if fitting with 3 degrees, return [[1, x, x**2, x**3],[1, x, x**2, x**3],...]\n",
"def create_matrix(x,d):\n",
" \"\"\"\n",
" x: predictor variable data(array)\n",
" d: dimension (int>0) including intercept\n",
" \"\"\"\n",
" X = np.zeros((len(x), d), float) # make a matrix size of (data size - dimension)\n",
" for i in range(d):\n",
" X[:,i] = x**i\n",
" return X\n",
"\n",
"def fitting(x,y,d):\n",
" d = d+1 # +1 means adding intercept\n",
" X = create_matrix(x,d) \n",
" (beta, residuals, rank, s) = np.linalg.lstsq(X, y) # numpy function for least square \n",
" return lambda x: reduce(lambda a,b:a+b, [beta[i]*(x**i) for i in range(d)])"
]
},
{
"cell_type": "code",
"execution_count": 231,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# 次元1,3,6,10としてそれぞれモデルを作成\n",
"model1 = fitting(x,y,1)\n",
"model3 = fitting(x,y,3)\n",
"model6 = fitting(x,y,6)\n",
"model10 = fitting(x,y,10)"
]
},
{
"cell_type": "code",
"execution_count": 232,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"ここで、対数尤度は以下のように表される。<br><br>\n",
"$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}\n",
"=-\\frac{n}{2}\\log{(2\\pi\\beta)}-\\frac{1}{2\\beta}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$\n",
"<br><br>\n",
"$$\\beta = \\frac{1}{n}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$"
],
"text/plain": [
"<IPython.core.display.Latex object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%latex\n",
"ここで、対数尤度は以下のように表される。<br><br>\n",
"$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}\n",
"=-\\frac{n}{2}\\log{(2\\pi\\beta)}-\\frac{1}{2\\beta}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$\n",
"<br><br>\n",
"$$\\beta = \\frac{1}{n}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ということで、モデル1,3,6,10それぞれの対数尤度を計算してみる。"
]
},
{
"cell_type": "code",
"execution_count": 233,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def log_likelihood(model,x,y):\n",
" n = len(x)\n",
" b = reduce(lambda a,b: a+b, [(y[i]-model(x[i]))**2 for i in range(n)])/n\n",
" l = -n*math.log(2*math.pi*b)/2 - \\\n",
" reduce(lambda a,b: a+b, [(y[i]-model(x[i]))**2 for i in range(n)])/(2*b)\n",
" return l"
]
},
{
"cell_type": "code",
"execution_count": 234,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Log likelihood for model with 1 param: 9.12883332651\n",
"Log likelihood for model with 3 param: 83.9990766856\n",
"Log likelihood for model with 6 param: 87.6848248313\n",
"Log likelihood for model with 10 param: 89.9458475187\n"
]
}
],
"source": [
"print('Log likelihood for model with 1 param: ', log_likelihood(model1,x,y))\n",
"print('Log likelihood for model with 3 param: ', log_likelihood(model3,x,y))\n",
"print('Log likelihood for model with 6 param: ', log_likelihood(model6,x,y))\n",
"print('Log likelihood for model with 10 param: ', log_likelihood(model10,x,y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ということで、やはり対数尤度は複雑なモデルほど大きくなっていることがわかる。d=10のモデルは明らかに過学習していることから、やはり対数尤度を最大化するモデルが最も良い、とは言えないだろうということがわかる。 \n",
"では、どのように良いモデルを選んだらよいのだろうか？\n",
" \n",
"ここからは、久保本のデータを使って話を進めていく。"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>y</th>\n",
" <th>x</th>\n",
" <th>f</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>6</td>\n",
" <td>8.31</td>\n",
" <td>C</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>6</td>\n",
" <td>9.44</td>\n",
" <td>C</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" y x f\n",
"0 6 8.31 C\n",
"1 6 9.44 C"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pd.read_csv('./data3a.csv')\n",
"df.head(2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Deviance (逸脱度)\n",
"どのように良いモデルを選んだらよいか、という問いに答える前にまず、当てはまりの良さを表す最大対数尤度を変形した統計量である逸脱度についてふれる。"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"対数尤度$$logL$$を最大にするパラメータを探すのが最尤推定法である。ここで、この最大対数尤度を$$logL^*$$\n",
"と表記すると、「当てはまりの悪さ」を表す指標である逸脱度は以下のように表される。<br><br>\n",
"$$D = -2 logL^*$$<br><br>\n",
"つまり、当てはまりの良さを表す最大対数尤度$$logL^*$$に-2をかけたもの。"
],
"text/plain": [
"<IPython.core.display.Latex object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%latex\n",
"対数尤度$$logL$$を最大にするパラメータを探すのが最尤推定法である。ここで、この最大対数尤度を$$logL^*$$\n",
"と表記すると、「当てはまりの悪さ」を表す指標である逸脱度は以下のように表される。<br><br>\n",
"$$D = -2 logL^*$$<br><br>\n",
"つまり、当てはまりの良さを表す最大対数尤度$$logL^*$$に-2をかけたもの。"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"試しに上のデータにおいて、$$y_i$$(種子数)が$$x_i$$(体サイズ)だけに依存するモデル\n",
"$$\\lambda_i = exp(\\beta_1 + \\beta_2 x_i)$$を作成し、逸脱度を求めてみる。"
],
"text/plain": [
"<IPython.core.display.Latex object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%latex\n",
"試しに上のデータにおいて、$$y_i$$(種子数)が$$x_i$$(体サイズ)だけに依存するモデル\n",
"$$\\lambda_i = exp(\\beta_1 + \\beta_2 x_i)$$を作成し、逸脱度を求めてみる。"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import statsmodels.api as sm\n",
"# add constant to the data\n",
"x_c = sm.add_constant(df.x)\n",
"# modeling with GLM (log link function is automatically selected for Poisson)\n",
"model = sm.GLM(df.y, x_c, family=sm.families.Poisson())\n",
"result = model.fit()\n",
"# get params\n",
"b1,b2 = result.params"
]
},
{
"cell_type": "code",
"execution_count": 137,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"ポアソン分布におけるパラメータλの対数尤度は以下である。<br><br>\n",
"$$logL(\\lambda) = log(exp(-n\\lambda) \\frac{\\lambda^{(y_1+y_2+...+y_n)}}{y_1!y_2!...y_n!})$$<br>\n",
"$$= -n\\lambda + (y_1+y_2+...+y_n)log\\lambda - log(y_1!y_2!...y_n!)$$\n",
"<br><br>\n",
"ここで、モデルのλは$$\\lambda_i = exp(\\beta_1 + \\beta_2 x_i)$$であったので、\n",
"このモデルの対数尤度は以下のように計算できる。\n",
"<br><br>\n",
"$$\\sum_{i=1}^n (-\\lambda_i + y_ilog(\\lambda_i) - log(y_i!))$$<br>\n",
"$$=\\sum_{i=1}^n (-exp(\\beta_1+\\beta_2x_i) + y_ilog(exp(\\beta_1+\\beta_2x_i)) - log(y_i!))$$"
],
"text/plain": [
"<IPython.core.display.Latex object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%latex\n",
"ポアソン分布におけるパラメータλの対数尤度は以下である。<br><br>\n",
"$$logL(\\lambda) = log(exp(-n\\lambda) \\frac{\\lambda^{(y_1+y_2+...+y_n)}}{y_1!y_2!...y_n!})$$<br>\n",
"$$= -n\\lambda + (y_1+y_2+...+y_n)log\\lambda - log(y_1!y_2!...y_n!)$$\n",
"<br><br>\n",
"ここで、モデルのλは$$\\lambda_i = exp(\\beta_1 + \\beta_2 x_i)$$であったので、\n",
"このモデルの対数尤度は以下のように計算できる。\n",
"<br><br>\n",
"$$\\sum_{i=1}^n (-\\lambda_i + y_ilog(\\lambda_i) - log(y_i!))$$<br>\n",
"$$=\\sum_{i=1}^n (-exp(\\beta_1+\\beta_2x_i) + y_ilog(exp(\\beta_1+\\beta_2x_i)) - log(y_i!))$$"
]
},
{
"cell_type": "code",
"execution_count": 128,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log likelihood: -235.38625077\n",
"Deviance: 470.77250154\n"
]
}
],
"source": [
"# function to define factorial\n",
"def factorial(n):return reduce(lambda x,y:x*y,[1]+list(range(1,n+1)))\n",
"\n",
"# function to calculate log likelihood in poisson\n",
"def poisson_log_likelihood(l,y):\n",
" return -l+y*math.log(l)-math.log(factorial(y))\n",
"\n",
"# calculate log likelihood\n",
"n = len(df)\n",
"ys = df.y.tolist()\n",
"lambdas = list(map(lambda x: math.exp(b1+b2*x), df.x.tolist()))\n",
"llh = reduce(lambda a,b:a+b,[\n",
" poisson_log_likelihood(lambdas[i],ys[i]) for i in range(n)\n",
" ])\n",
"print('log likelihood:',llh)\n",
"\n",
"# calculate deviance\n",
"dev = -2*llh\n",
"print('Deviance:',dev)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Statsmodelsの機能で見てみる。"
]
},
{
"cell_type": "code",
"execution_count": 127,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"84.992996490728402"
]
},
"execution_count": 127,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"result.deviance"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"結果が異なっている。。が、これは問題ない。 \n",
"ここで計算されたdevianceは上で定義した逸脱度ではなく、Residual Deviance (残差逸脱度)と呼ばれるものである。 \n",
"Residual Devianceは以下のように定義される。 \n",
" \n",
"**残差逸脱度 = 逸脱度 - (ポアソン分布モデルで可能な最小逸脱度)** \n",
" \n",
"「ポアソン分布モデルで可能な最小逸脱度」とは、データ数の数だけパラメータを使って当てはめたモデルの逸脱度のことである。これはつまり、λ<sub>1</sub>=y<sub>1</sub>, λ<sub>2</sub>=y<sub>2</sub>,...,λ<sub>n</sub>=y<sub>n</sub>と、データ数分のパラメータを設定して作成したモデルである。これは統計モデルとしての意味はなさないが、他のどのモデルよりも「当てはまりのよさ」である対数尤度は大きくなる。つまり逸脱度は最も小さくなる。\n",
"\n",
"詳しくは久保本のp72参照。"
]
},
{
"cell_type": "code",
"execution_count": 135,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"fullmodel log likelihood: -192.889752524\n",
"Possible Minimum Deviance: 385.779505049\n",
"Residual Deviance: 84.9929964907\n"
]
}
],
"source": [
"# calculate log likelihood as λ=y\n",
"fullmodel_llh = reduce(lambda a,b:a+b,[\n",
" poisson_log_likelihood(ys[i],ys[i]) for i in range(n)\n",
" ])\n",
"print('fullmodel log likelihood:',fullmodel_llh)\n",
"\n",
"# calculate possible minimum deviance\n",
"min_dev = -2*fullmodel_llh\n",
"print('Possible Minimum Deviance:',min_dev)\n",
"res_dev = dev - min_dev\n",
"print('Residual Deviance:',res_dev)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上のStatsmodelsで出したdevianceの値と、Residual Devianceの値が一致していることがわかる。 \n",
"統計モデルのパラメータを増やしていけば、逸脱度は小さくなっていき、この残差逸脱度の値も小さくなっていく。 \n",
" \n",
"ここで先ほどとは反対に、(パラメータを最尤推定するGLMにおいて)最も逸脱度が大きくなるモデルについて考えてみる。すなわち最もあてはまりの悪いモデルである。この観測データに対するポアソン回帰の場合、最もパラメータ数の少ないモデル、すなわち切片β<sub>1</sub>だけのモデル(パラメータ数k=1)である。これは null model と呼ばれる。"
]
},
{
"cell_type": "code",
"execution_count": 141,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# modeling with GLM (log link function is automatically selected for Poisson)\n",
"nullmodel = sm.GLM(df.y, np.ones(len(df)), family=sm.families.Poisson())\n",
"result_null = nullmodel.fit()\n",
"# get params\n",
"b1_null = result_null.params"
]
},
{
"cell_type": "code",
"execution_count": 149,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"nullmodel log likelihood: -237.643221309\n",
"Deviance of nullmodel: 475.286442619\n",
"Null Deviance: 89.5069375696\n"
]
}
],
"source": [
"# calculate deviance of null model\n",
"lambda_null = math.exp(b1_null)\n",
"nullmodel_llh = reduce(lambda a,b:a+b,[\n",
" poisson_log_likelihood(lambda_null,ys[i]) for i in range(n)\n",
" ])\n",
"print('nullmodel log likelihood:',nullmodel_llh)\n",
"\n",
"# calculate deviance of null model\n",
"nullmodel_dev = -2*nullmodel_llh\n",
"print('Deviance of nullmodel:',nullmodel_dev)\n",
"null_dev = nullmodel_dev - min_dev\n",
"print('Null Deviance:',null_dev)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上のように、一定モデルで計算した残差逸脱度(Nullmodelの逸脱度 - Fullmodelの逸脱度)、つまり最大逸脱度-最小逸脱度を、Null Devianceと呼ぶ。 \n",
"Statsmodelsの機能でも以下のように出すことができる。"
]
},
{
"cell_type": "code",
"execution_count": 150,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"89.506937569581211"
]
},
"execution_count": 150,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"result.null_deviance"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ここまでの議論から、最尤推定するパラメータ数(k)さえ増やしていけば、残差逸脱度はどんどん小さくなり、当てはまりがよくなることがわかる。 \n",
"しかし、これは「たまたま得られたデータへのあてはめ向上を目的とする特殊化」であり、その統計モデルの「予測の良さ」を損なっているかもしれない。ここで登場するのが**AIC**である。 \n",
" \n",
"#### AIC (Akaike's information criterion)\n",
"複数の統計モデルの中から何らかの基準で良いモデルを選択することをモデル選択と呼ぶ。よく使われているモデル選択基準の一つが**AIC**である。これは、統計モデルのあてはまりの良さ(goodness of fit)ではなく、予測の良さ(goodness of prediction)を重視する選択基準である。AICは以下のように定義され、この値が最も小さいモデルが良いモデルとなる。"
]
},
{
"cell_type": "code",
"execution_count": 156,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"AIC = -2(最大対数尤度 - 最尤推定したパラメータ数)<br><br>\n",
"$$= -2(logL^* - k)$$<br>\n",
"$$= D + 2k$$"
],
"text/plain": [
"<IPython.core.display.Latex object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%latex\n",
"AIC = -2(最大対数尤度 - 最尤推定したパラメータ数)<br><br>\n",
"$$= -2(logL^* - k)$$<br>\n",
"$$= D + 2k$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 試してみる\n",
"まず、上のポアソン回帰モデルのAICを見てみる。"
]
},
{
"cell_type": "code",
"execution_count": 163,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# calc AIC\n",
"def aic(d,k): return d+2*k"
]
},
{
"cell_type": "code",
"execution_count": 170,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AIC Null model: 477.286442619\n",
"AIC Full model: 585.779505049\n",
"AIC x model: 474.77250154\n"
]
}
],
"source": [
"# null model \n",
"print('AIC Null model:', aic(nullmodel_dev,1))\n",
"# full model\n",
"print('AIC Full model:', aic(min_dev,100))\n",
"# model with x\n",
"print('AIC x model:', aic(dev,2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ちなみに、Statsmodelsのresultでは以下のようにもAICの値を出せる。"
]
},
{
"cell_type": "code",
"execution_count": 173,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"474.77250154\n",
"477.286442619\n"
]
}
],
"source": [
"print(result.aic)\n",
"print(result_null.aic)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"次に、この記事の上の方の、異なる次数の最小二乗法モデルでAICの違いを見てみる。"
]
},
{
"cell_type": "code",
"execution_count": 240,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# models with 1~10 degrees\n",
"models = [fitting(x,y,i) for i in range(1,11)]\n",
"llhs = [log_likelihood(m,x,y) for m in models]\n",
"devs = [-2*llh for llh in llhs]\n",
"aics = [devs[i]+2*(i+2) for i in range(10)] # k=i+2 because of intercept"
]
},
{
"cell_type": "code",
"execution_count": 241,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[-14.25766665301343,\n",
" -12.267023221608298,\n",
" -159.99815337126734,\n",
" -158.00136831437828,\n",
" -163.17197064106125,\n",
" -161.3696496626319,\n",
" -160.55170023114999,\n",
" -158.6132787166801,\n",
" -159.86845008863213,\n",
" -157.89169503737139]"
]
},
"execution_count": 241,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"aics"
]
},
{
"cell_type": "code",
"execution_count": 242,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAGnCAYAAAAUmzt8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUnXV97/H3N5lcIAmTBDGBBDQRggECZLToqbZOi6L0\noniWB6ldi1Zdtktty7FaK9JVcrq6qta6qu1Sz+k61kuXVsXWViyiWBlbu45ARUhCuIQimHANJJlw\nC5Dke/549iQ7yUwyM/vy25f3a6295tnP3vPs784wm8/8ft/n90RmIkmSpPabUboASZKkfmUQkyRJ\nKsQgJkmSVIhBTJIkqRCDmCRJUiEGMUmSpEKKBbGIeF1E3BERd0XEH5aqQ5IkqZQosY5YRMwA7gLO\nBx4AbgIuycw72l6MJElSIaVGxM4DNmfmfZn5HPBl4A2FapEkSSqiVBBbBmypu7+1tk+SJKlvDJQu\n4EgiwusvSZKkrpGZMZXnlxoRux84pe7+8tq+w2Rm392uvPLK4jX4vn3fvm/ft+/b9+37ntptOkoF\nsZuAUyPiBRExG7gE+EahWiRJkoooMjWZmXsj4neA71CFwc9k5u0lapEkSSqlWI9YZl4LnF7q9TvZ\n8PBw6RKK8H33F993f/F995d+fd/TUWQdscmKiOzk+iRJksZEBNklzfqSJEl9zyAmSZJUiEFMkiSp\nEIOYJElSIQYxSZKkQgxikiRJhRjEJEmSCjGISZIkFWIQkyRJKsQgJkmSVIhBTJIkqRCDmCRJUiEG\nMUmSpEIMYpIkSYUYxCRJkgoxiEmSJBViEJMkSSrEICZJklSIQUySJKkQg5gkSVIhA6UL6BeZ8NRT\nsH37gduOHQff374dHn8cPv5xWLKkdMWSJKnVDGJTtG8fjI4eHqAmClb1+wcGYNEiWLz48NuiRbBi\nBXz+8zAyAm9+c+l3KkmSWq1vg9izzx45PE302K5dsGDBxIFq+XJYs2b8oDV37tHreuQRuPlmg5gk\nSf2gq4NYJjz55NRGpcZuzzxzcEg6NDitXj3+YwsXwsyZrXtPQ0PwiU+07viSJKlzRGaWrmFCEZF/\n8id5xGA1a9bhI08TTf3V358/HyJKv8PDPfAAnH02bNvWmfVJkqTxRQSZOaX/e3f8iNgzz8DJJ8M5\n5xweqiY73ddNTjyx6iXburV635IkqXd1fBD70z8tXUF7RcDatVWfmEFMkqTe5jpiHWhoCH7849JV\nSJKkVjOIdaCxETFJktTbDGIdyBExSZL6g0GsA61YUa2wv21b6UokSVIrGcQ60FjDvqNikiT1NoNY\nh7JPTJKk3mcQ61D2iUmS1PtaFsQi4s8j4vaIuCUi/iEijqt77PKI2Fx7/IJW1dDNHBGTJKn3tXJE\n7DvAmZl5LrAZuBwgIs4ALgZWAxcCn4rwYj6HOv306nJHu3aVrkSSJLVKy4JYZn43M/fV7v4QWF7b\nfj3w5czck5n3UoW081pVR7caGIA1a+CWW0pXIkmSWqVdPWJvA66pbS8DttQ9dn9tnw5hn5gkSb2t\noWtNRsR1wJL6XUACV2Tm1bXnXAE8l5l/P53XWLdu3f7t4eFhhoeHp1tu1xkagn//99JVSJKk8YyM\njDAyMtLQMSIzm1PNeAeP+E3gHcAvZuYztX0fADIzP1K7fy1wZWbeMM73Zyvr63Q/+hG89a2wfn3p\nSiRJ0tFEBJk5pb73lgWxiHgd8DHg5zPzsbr9ZwBfBF5GNSV5HXDaeImr34PYM8/AokXw2GNwzDGl\nq5EkSUcynSDWyh6xvwbmA9dFxM0R8SmAzNwEfBXYRNU39q6+TltHMGcOrFoFGzeWrkSSJLVCS6cm\nG9XvI2IAb3sbvOxl8Nu/XboSSZJ0JJ02IqYm8JqTkiT1LoNYhxsacoV9SZJ6lVOTHe6JJ2DJEti5\nE2bNKl2NJEmaiFOTPWj+fDj5ZLjjjtKVSJKkZjOIdQH7xCRJ6k0GsS5gn5gkSb3JINYFHBGTJKk3\n2azfBR57DFauhB07YIbRWZKkjmSzfo86/nhYuBDuuad0JZIkqZkMYl3CPjFJknqPQaxLrF1rEJMk\nqdcYxLrE0JAN+5Ik9RqDWJcYGxHz3AVJknqHQaxLnHRSdcbk/feXrkSSJDWLQaxLRNgnJklSrzGI\ndRH7xCRJ6i0GsS7iiJgkSb3FINZFHBGTJKm3GMS6yIoVMDoKjz5auhJJktQMBrEuMmOGFwCXJKmX\nGMS6jH1ikiT1DoNYl7FPTJKk3mEQ6zKOiEmS1DsiO/iaORGRnVxfCXv2wOAgPPggHHdc6WokSdKY\niCAzYyrf44hYlxkYgLPOgltvLV2JJElqlEGsC9knJklSbzCIdSH7xCRJ6g0GsS7kiJgkSb3BZv0u\ntHs3LFoEO3bA3Lmlq5EkSWCzft+YOxdWrYKNG0tXIkmSGmEQ61L2iUmS1P0MYl3KPjFJkrqfQaxL\nOSImSVL3s1m/Sz3+OCxdCqOj1SKvkiSprI5s1o+I90bEvohYXLfv8ojYHBG3R8QFra6hFy1YAMuW\nwR13lK5EkiRNV0uDWEQsB14D3Fe3bzVwMbAauBD4VERMKT2qYp+YJEndrdUjYn8J/MEh+94AfDkz\n92TmvcBm4LwW19GT7BOTJKm7tSyIRcTrgS2ZueGQh5YBW+ru31/bpylyREySpO7WUJt3RFwHLKnf\nBSTwR8AHqaYlG7Ju3br928PDwwwPDzd6yJ6xdm0VxPbtgxme/ypJUluNjIwwMjLS0DFactZkRJwF\nfBd4iiqcLaca+ToPeBtAZn649txrgSsz84ZxjuNZk0dxyilw/fXwoheVrkSSpP7WMWdNZubGzFya\nmSszcwWwFVibmY8A3wDeHBGzI2IFcCpwYyvq6AdDQ/aJSZLUrdo1oZVUI2Nk5ibgq8Am4BrgXQ57\nTd/Y9KQkSeo+Luja5a6+Gj75Sbj22tKVSJLU3zpmalLtM7aEhXlVkqTuYxDrcstqC3888EDZOiRJ\n0tQZxLpchH1ikiR1K4NYD/DMSUmSupNBrAc4IiZJUncyiPUAR8QkSepOBrEesHIl7NwJjz1WuhJJ\nkjQVBrEeMGMGnHuu05OSJHUbg1iPsE9MkqTuYxDrEfaJSZLUfQxiPcIRMUmSuo/XmuwRe/bA4CA8\n/DDMn1+6GkmS+o/XmuxjAwNw5plw662lK5EkSZNlEOsh9olJktRdDGI9xD4xSZK6i0GshzgiJklS\nd7FZv4fs3g2LF8OOHTBnTulqJEnqLzbr97m5c+HUU2HjxtKVSJKkyTCI9Rj7xCRJ6h4GsR5jn5gk\nSd3DINZjHBGTJKl72KzfY3btghNPhNHRapFXSZLUHjbri+OOg5NOgjvvLF2JJEk6GoNYDxoacnpS\nkqRuYBDrQWvX2rAvSVI3MIj1IEfEJEnqDjbr96Bt2+C006oV9mNKLYOSJGm6bNYXACecAAsWwE9+\nUroSSZJ0JAaxHmWfmCRJnc8g1qPsE5MkqfMZxHqUI2KSJHU+g1iPGrvmpOc6SJLUuQxiPWr5cti7\nFx58sHQlkiRpIgaxHhVhn5gkSZ2upUEsIn43Im6PiA0R8eG6/ZdHxObaYxe0soZ+Zp+YJEmdbaBV\nB46IYeBXgTWZuScinlfbvxq4GFgNLAe+GxGnuXJr8w0NwVe+UroKSZI0kVaOiL0T+HBm7gHIzEdr\n+98AfDkz92TmvcBm4LwW1tG3xhr2JUlSZ2plEFsF/HxE/DAiro+Il9T2LwO21D3v/to+NdmLXgTb\nt1c3SZLUeRqamoyI64Al9buABP6oduxFmfnyiPgZ4Cpg5VRfY926dfu3h4eHGR4ebqDi/jJjBpx7\nbtWwf/75pauRJKm3jIyMMDIy0tAxWnbR74i4BvhIZn6/dn8z8HLgHQCZ+eHa/muBKzPzhnGOYetY\ngy67DE4+Gd73vtKVSJLU2zrtot//BPwiQESsAmZn5mPAN4A3R8TsiFgBnArc2MI6+pp9YpIkda5W\nBrHPAisjYgPwJeBSgMzcBHwV2ARcA7zLYa/WWbvWtcQkSepULZuabAanJhv33HOwcCE8/DDMn1+6\nGkmSelenTU2qA8yaBWecAevXl65EkiQdyiDWB+wTkySpMxnE+oB9YpIkdSaDWB9wREySpM5ks34f\nePppOP542LED5swpXY0kSb3JZn2N65hjqssd3XZb6UokSVI9g1ifsE9MkqTOYxDrE/aJSZLUeQxi\nfcIRMUmSOo/N+n1idBSWLau+zpxZuhpJknqPzfqa0OAgLF0Kd91VuhJJkjTGINZH7BOTJKmzGMT6\niH1ikiR1FoNYH3FETJKkzmKzfh955BE4/XTYvh1iSq2EkiTpaGzW1xE9//kwbx7ce2/pSiRJEhjE\n+o59YpIkdQ6DWJ+xT0ySpM5hEOszjohJktQ5DGJ9xhExSZI6h0Gsz5x8Mjz7LDz4YOlKJEmSQazP\nRFSjYk5PSpJUnkGsD9knJklSZzCI9SH7xCRJ6gwGsT7kiJgkSZ3BSxz1oX37YHAQfvpTWLSodDWS\nJPUGL3GkSZkxA845B265pXQlkiT1N4NYn7JPTJKk8gxifco+MUmSyjOI9SlHxCRJKs9m/T713HNV\nw/62bTBvXulqJEnqfjbra9JmzYIzzoD160tXIklS/zKI9TH7xCRJKssg1sfsE5MkqSyDWB/z4t+S\nJJXVsiAWEedExP+LiB9HxI0R8dK6xy6PiM0RcXtEXNCqGnRka9bA7bfDs8+WrkSSpP7UyhGxPweu\nzMy1wJXARwEi4gzgYmA1cCHwqYiY0hkGao5jj4WVK+G220pXIklSf2plENsHDNa2FwL317ZfD3w5\nM/dk5r3AZuC8FtahI7BhX5KkcloZxN4D/EVE/JRqdOzy2v5lwJa6591f26cCbNiXJKmcgUa+OSKu\nA5bU7wISuAJ4NXBZZv5TRLwJ+FvgNVN9jXXr1u3fHh4eZnh4uIGKdai1a+FrXytdhSRJ3WdkZISR\nkZGGjtGylfUjYmdmLjz0fkR8AMjM/Eht/7VUvWQ3jHMMV9ZvsZ07YflyGB2FmTNLVyNJUvfqtJX1\n74+IVwFExPlUvWAA3wAuiYjZEbECOBW4sYV16AgWLoQlS2Dz5qM/V5IkNVdDU5NH8Q7gryJiJrAb\n+C2AzNwUEV8FNgHPAe9y2KussT6xF7+4dCWSJPUXL/ot/uzPYMcO+OhHS1ciSVL36rSpSXUJz5yU\nJKkMR8TEww/D6tXw2GPg0rqSJE2PI2KaliVL4Jhj4L77SlciSVJ/MYgJcIV9SZJKMIgJsE9MkqQS\nDGICHBGTJKkEg5gAR8QkSSrBICYATjkFdu+Ghx4qXYkkSf3DICagWrZiaMjpSUmS2skgpv3sE5Mk\nqb0MYtrPPjFJktrLIKb9HBGTJKm9vMSR9tu7FxYuhC1bqq+SJGnyvMSRGjJzJpx9NtxyS+lKJEnq\nDwYxHcQ+MUmS2scgpoPYJyZJUvsYxHQQR8QkSWofm/V1kGefrRr1H30Ujj22dDWSJHUPm/XVsNmz\nYfVqWL++dCWSJPU+g5gOY5+YJEntYRDTYewTkySpPQxiOowjYpIktYfN+jrMk0/CCSfA6CjMmlW6\nGkmSuoPN+mqKefPghS+ETZtKVyJJUm8ziGlc9olJktR6BjGNyz4xSZJazyCmcTkiJklS69msr3Ht\n2AGnnFI17M8wrkuSdFQ266tpFi2qzpzcvLl0JZIk9S6DmCZkn5gkSa1lENOE7BOTJKm1DGKakCNi\nkiS1ls36mtBDD8GZZ8Kjj0JMqfVQkqT+Y7O+mmrpUpgzB7ZsKV2JJEm9qaEgFhFvioiNEbE3IoYO\neezyiNgcEbdHxAV1+4ciYn1E3BURH2/k9dV6a9faJyZJUqs0OiK2AXgj8P36nRGxGrgYWA1cCHwq\nYv/k1qeBt2fmKmBVRLy2wRrUQkND9olJktQqDQWxzLwzMzcDh86HvgH4cmbuycx7gc3AeRGxFFiQ\nmTfVnvcF4KJGalBrOSImSVLrtKpHbBlQ31l0f23fMmBr3f6ttX3qUI6ISZLUOgNHe0JEXAcsqd8F\nJHBFZl7dqsLGrFu3bv/28PAww8PDrX5J1XnBC+Cpp+CRR+D5zy9djSRJnWNkZISRkZGGjtGU5Ssi\n4nrgvZl5c+3+B4DMzI/U7l8LXAncB1yfmatr+y8BXpWZ75zguC5f0QHOPx/e/354rd18kiRNqPTy\nFfUv/A3gkoiYHRErgFOBGzPzIWA0Is6rNe9fCvxzE2tQC9gnJklSazS6fMVFEbEFeDnwzYj4FkBm\nbgK+CmwCrgHeVTe09W7gM8BdwObMvLaRGtR69olJktQarqyvo7r9dvjVX4W77y5diSRJnWs6U5MG\nMR3V3r2wcCFs3QqDg6WrkSSpM5XuEVOPmjkT1qyBW24pXYkkSb3FIKZJsU9MkqTmM4hpUjxzUpKk\n5jOIaVIcEZMkqfls1tekPPts1bD/6KNw7LGlq5EkqfPYrK+WmT0bXvxi2LChdCWSJPUOg5gmbe1a\npyclSWomg5gmbWjIhn1JkprJIKZJc0RMkqTmsllfk/bkk3DCCTA6CrNmla5GkqTOYrO+WmrePHjB\nC6prT0qSpMYZxDQl9olJktQ8BjFNiX1ikiQ1j0FMU+KImCRJzWOzvqZk+3Z44Qth506YYYyXJGk/\nm/XVcosXw/HHw913l65EkqTuZxDTlNknJklScxjENGX2iUmS1BwGMU2ZI2KSJDWHzfqasgcfhDVr\nYNs2iCm1JEqS1Lts1ldbnHgiDAzA1q2lK5EkqbsZxDQt9olJktQ4g5imxT4xSZIaZxDTtDgiJklS\n4wximhZHxCRJapxBTNOyYgU8/nh15qQkSZoeg5imJcJRMUmSGmUQ07TZJyZJUmMMYpo2R8QkSWqM\nQUzT5oiYJEmN8RJHmra9e2FwEB54AI47rnQ1kiSV5SWO1FYzZ1bXnLzlltKVSJLUnQxiasjQkH1i\nkiRNl0FMDVm71j4xSZKmq6EgFhFvioiNEbE3Iobq9r86Iv4zIm6NiJsi4hfqHhuKiPURcVdEfLyR\n11d5johJkjR9jY6IbQDeCHz/kP3bgF/JzHOA3wT+ru6xTwNvz8xVwKqIeG2DNaigM8+Eu++Gp58u\nXYkkSd2noSCWmXdm5mYgDtl/a2Y+VNu+DZgbEbMiYimwIDNvqj31C8BFjdSgsubMgdNPh40bS1ci\nSVL3aXmPWES8Cbg5M58DlgFb6x7eWtunLmafmCRJ0zNwtCdExHXAkvpdQAJXZObVR/neM4EPAa+Z\nboHr1q3bvz08PMzw8PB0D6UWsU9MktSPRkZGGBkZaegYTVnQNSKuB96bmTfX7VsO/CvwG5n5w9q+\npcD1mbm6dv8S4FWZ+c4JjuuCrl3gP/4D3vMeuPHG0pVIklRO6QVd979wRAwC3wT+cCyEAdT6xkYj\n4ryICOBS4J+bWIMKOOccuO02eO650pVIktRdGl2+4qKI2AK8HPhmRHyr9tDvAC8C/jgifhwRN0fE\n82qPvRv4DHAXsDkzr22kBpU3fz6cfDLccUfpSiRJ6i5ea1JN8Za3wOteB5deWroSSZLKKD01qT7m\nmZOSJE2dQUxN4ZmTkiRNnVOTaort22HFCtixA2YY7yVJfcipSRWzeDEsWgT33FO6EkmSuodBTE1j\nn5gkSVNjEFPT2CcmSdLUGMTUNI6ISZI0NQYxNc3YiJjnV0iSNDkGMTXNiSdCBNx/f+lKJEnqDgYx\nNU2EfWKSJE2FQUxNZZ+YJEmTZxBTUzkiJknS5BnE1FSOiEmSNHkGMTXVypWwaxc8+mjpSiRJ6nwG\nMTVVBJx7rtOTkiRNhkFMTWefmCRJk2MQU9PZJyZJ0uQYxNR0Q0MGMUmSJiOyg69HExHZyfVpfHv2\nwOAgPPggHHdc6WokSWqPiCAzYyrf44iYmm5gANasgVtvLV2JJEmdzSCmlli71oZ9SZKOxiCmlrBP\nTJKkozOIqSUcEZMk6ehs1ldL7N4NixfD9u0wd27paiRJaj2b9dUx5s6F006DjRtLVyJJUucyiKll\n7BOTJOnIDGJqGfvEJEk6MoOYWsYRMUmSjsxmfbXM44/D0qUwOlot8ipJUi+zWV8dZcECWL4c7rij\ndCWSJHUmg5haamjIPjFJkiZiEFNLrV1rn5gkSRMxiKmlHBGTJGliNuurpR57DFauhB07YIaxX5LU\nw9rerB8Rb4qIjRGxNyKGxnn8lIh4PCJ+v27fUESsj4i7IuLjjby+Ot/xx8PChfCTn5SuRJKkztPo\nGMUG4I3A9yd4/GPANYfs+zTw9sxcBayKiNc2WIM6nH1ikiSNr6Eglpl3ZuZm4LBhuIh4A3APcFvd\nvqXAgsy8qbbrC8BFjdSgzmefmCRJ42tJ105EzAPeD/wvDg5py4Ctdfe31vaphzkiJknS+I663nlE\nXAcsqd8FJHBFZl49wbetA/4yM5+KmFLP2uEHWrdu//bw8DDDw8MNHU/tN3apo0xo8D8HSZI6xsjI\nCCMjIw0doylnTUbE9cB7M/Pm2v1/A5bXHl4E7AX+GPhH4PrMXF173iXAqzLznRMc17Mme0AmLFlS\nTU8uc/xTktSjpnPWZDOvALj/hTPz5+uKuhJ4PDM/Vbs/GhHnATcBlwJ/1cQa1IEiDvSJGcQkSTqg\n0eUrLoqILcDLgW9GxLcm8W3vBj4D3AVszsxrG6lB3cE+MUmSDueCrmqLq66CL30Jvv710pVIktQa\nbV/QVZosR8QkSTqcQUxtsXIl7NxZXfJIkiRVDGJqixkz4NxzXdhVkqR6BjG1jSvsS5J0MIOY2sY+\nMUmSDmYQU9s4IiZJ0sFcvkJts2cPDA7Cww/D/Pmlq5EkqblcvkIdbWAAzjoLbr21dCWSJHUGg5ja\nyj4xSZIOMIiprewTkyTpAIOY2soRMUmSDrBZX221ezcsXgw7dsCcOaWrkSSpeWzWV8ebOxdOPRU2\nbixdiSRJ5RnE1Hb2iUmSVDGIqe3sE5MkqWIQU9s5IiZJUsVmfbXdrl1w0kkwOgozZ5auRpKk5rBZ\nX13huOOqIHbnnaUrkSSpLIOYirBPTJIkg5gKsU9MkiSDmApxREySJJv1Vci2bXDaadUK+zGltkZJ\nkjqTzfrqGiecAAsWwE9+UroSSZLKMYipGPvEJEn9ziCmYuwTkyT1O4OYinFETJLU7wxiKsYRMUlS\nvzOIqZjly2HvXnjwwdKVSJJUhkFMxURU05OOikmS+pVBTEWtXWufmCSpfxnEVJQjYpKkfmYQU1E2\n7EuS+plBTEW96EXVZY62by9diSRJ7TdQugD1txkz4Jxzqj6x888vXU3r7NkDO3dWgbP+NjoKJ59c\n/RuccorX3ZSkTpMJTz998Gf32ADCodvT0VAQi4g3AeuA1cDPZObNdY+dDfxv4Dhgb+3xZyNiCPgc\nMBe4JjP/ZyM1qPuNLezaDUFs9+7Df/GOdBt73hNPwOAgLF584LZoERx3HHzzm7B+PTz5JKxZU4Wy\ns8+ubmvWwLx5pd+1pKl6+unqD62Bgep3eO5c/9Aqbe/e6g/iiULUkQLWzJkHf3Yfun3KKdX2VVdN\nva7IzGm/qYg4HdgH/B/gfWNBLCJmAjcDv56ZGyNiEbAzMzMibgB+JzNviohrgE9k5rcnOH42Up+6\nw+c/D9/+NnzpS+15vcwqGB0pOE1027Pn4DB1pFv9L+vgYDX6dyTbtsGGDVUou/XW6uvtt8OyZVUo\nqw9oL3zh0Y8nqTHPPXdgJHvHjgOfD/VfJ9retw8WLqw+M558sjrWvHkH3+bPb8792bNL/0u119jo\n1GRD1Nj9xx+v/vg9UqCa6LG5cydXW0SQmVOK3A0FsboXvh54b10QuxD4tcy89JDnLQW+l5ln1O5f\nArwqM985wXENYn1gwwa4+OIqdEzF3r3VX5yTHZWqv82Zc/TwNN7t2GPb+1ftnj1w111VKKsPaKOj\n1WhZfUBbswYWLGhfbVI32Lev+n2ZapDavr36H/7ChQc+F+q/Hm37mGMO/qwYC2T1tyeeaOz+2D5o\nTrAbb99AixqYxn4ukw1R9duZRw5NjfxB3KhOCmKXAS8Bng88D/hKZn40Il4CfCgzL6g975XA+zPz\n9RMc1yDWB557rvqw+8EPqg+WyU771f91c7QAVf9LuWhRFcS62fbtB8LZWEDbtAmWLj0wajYW0Fau\ndPRM3S2z+myYTHg6dN+uXVWgmGyQqt+3YEF3/O48+2zzwt2h92fNmnq4O/bYg39e4wWqXbuq504n\nUB1zTOl/8YlNJ4gdNetGxHXAkvpdQAJXZObVRzjuK4CXAruBf42I/wR2TaU4gHXr1u3fHh4eZnh4\neKqHUIebNQt++ZfhLW8ZPzytXj3+/sHBat6+Hy1eDMPD1W3M3r1w990HRs0+97lqe/t2OOusgwPa\nmjXVv5/UTrt3Tz1Ijd1mzTpyaFq+fPxwNTjYulGdTjF7dnVbtKi5x82sfmZTCW6PPQZPPXUg/B76\n+T32c1m4sDc+v0dGRhgZGWnoGK0aEXsz8LrMfGvt/h8BTwNfBK7PzNW1/U5NSi22c+fhvWcbN8Lz\nnndw39k551TLifTCh6MOt2cPPPNMddu9+8D2ofdbtT06WtUwlWm+XhrFVn8oPTX5vsz8Ue3+QuC7\nwCuBPcC3gI9l5rUR8UPg94CbgH8B/iozr53guAYxqQX27YP/+q/De88efhjOPPPggHb22c3/S7vf\n7NtXTaXv3Fn1HZUIQ/v2VQ3Hc+Yc+Dp2q7/fiu05c6oRkHb3WErt1vYgFhEXAX9N1Qe2E7glMy+s\nPfYW4INUZ1X+S2ZeXtv/Eg5evuKyIxzfICa10a5dB0bPxgLahg1VEDv0zM3TTuv9KZ8xYz1KO3ce\nOP19bPvQ23iP7dpVhZDBwepruwJQ/fbAgCFIarViI2KtYhCTytu3D+6998Co2djtgQeq/o9DA9rx\nx5eu+HBjvS4TBajJ7B8b1Vm0qPp66G2i/QsX9kefkiSDmKQ2euKJqtfs0IC2YMHhZ26uWlU1Wzfi\nmWemHp7qbzNmTD44HfrY4GD/rdUkaeoMYpKKyoT77ju892zLFnjxiw8EtLPOqqbJpjLVt3fv+OFp\nMoFqcHDk5aWNAAAG10lEQVTyCzJK0nQZxCR1pCefhNtuOxDQNm6szs6cbKgaW9naHidJncwgJkmS\nVMh0glgXrBksSZLUmwxikiRJhRjEJEmSCjGISZIkFWIQkyRJKsQgJkmSVIhBTJIkqRCDmCRJUiEG\nMUmSpEIMYpIkSYUYxCRJkgoxiEmSJBViEJMkSSrEICZJklSIQUySJKkQg5gkSVIhBjFJkqRCDGKS\nJEmFGMQkSZIKMYhJkiQVYhCTJEkqxCAmSZJUiEFMkiSpEIOYJElSIQYxSZKkQgxikiRJhRjEJEmS\nCjGISZIkFWIQkyRJKsQgJkmSVIhBTJIkqZCGglhEvCkiNkbE3ogYqts/EBGfi4j1EXFbRHyg7rGh\n2v67IuLjjbx+rxoZGSldQhG+7/7i++4vvu/+0q/vezoaHRHbALwR+P4h+/8HMDszzwZeCvx2RJxS\ne+zTwNszcxWwKiJe22ANPadf/wP2ffcX33d/8X33l35939PRUBDLzDszczMQhz4EzIuImcCxwDPA\nrohYCizIzJtqz/sCcFEjNUiSJHWrVvWIfQ14CngQuBf4i8zcCSwDttY9b2ttnyRJUt+JzDzyEyKu\nA5bU76Ia8boiM6+uPed64L2ZeXPt/s8C7wR+Azge+HfgdbXtD2XmBbXnvRJ4f2a+foLXPnJxkiRJ\nHSQzD50lPKKBSRzwNdOo4y3AtZm5D9gWEf9B1Sv2A+DkuuctB+4/wmtP6c1IkiR1k2ZOTdaHpp8C\nvwgQEfOAlwO3Z+ZDwGhEnBcRAVwK/HMTa5AkSeoajS5fcVFEbKEKWt+MiG/VHvoksCAiNgI3AJ/J\nzNtqj70b+AxwF7A5M69tpAZJkqRuddQeMUmSJLVGR66sHxGfiYiHI2J96VraJSKWR8T3agvgboiI\n3ytdUztExJyIuCEiflx731eWrqmdImJGRNwcEd8oXUu7RMS9EXFr7Wd+Y+l62iUiBiPiqoi4vfZ7\n/rLSNbVaRKyq/Zxvrn0d7aPPtvfUFjxfHxFfjIjZpWtqh4i4rPZZ3tP/Hxsvp0TEooj4TkTcGRHf\njojByRyrI4MY8Fmg3xZ63QP8fmaeCfw34N0R8eLCNbVcZj4D/EJmrgXOBS6MiPMKl9VOlwGbShfR\nZvuA4cxcm5n99LP+BHBNZq4GzgFuL1xPy2XmXbWf8xDwEuBJ4OuFy2q5iDgJ+F1gqLaw+QBwSdmq\nWi8izgTeTnVy3rnAr0TEyrJVtcx4OeUDwHcz83Tge8DlkzlQRwaxzPwBsKN0He2UmQ9l5i217Seo\nPqT7Yo21zHyqtjmH6gOrL+bLI2I58EvA/y1dS5sFHfrZ0yoRcRzwc5n5WYDM3JOZuwqX1W6vBv4r\nM7eULqRNZlItbD5AtbD5A4XraYfVwA2Z+Uxm7gX+DfjvhWtqiQlyyhuAz9e2P88kF6zvqw/DbhER\nL6T6a+KGspW0R2167sfAQ8B1dVde6HV/CfwBfRI86yRwXUTcFBHvKF1Mm6wAHo2Iz9am6f4mIo4p\nXVSbvRn4+9JFtENmPgB8jGoFgfuBnZn53bJVtcVG4OdqU3THUv2hefJRvqeXPD8zH4ZqcAV4/mS+\nySDWYSJiPtWVCS6rjYz1vMzcV5uaXA68LCLOKF1Tq0XELwMP10ZBg8MvE9bLXlGbqvolqin4V5Yu\nqA0GgCHgk7X3/hTVNEZfiIhZwOuBq0rX0g4RsZBqdOQFwEnA/Ih4S9mqWi8z7wA+AlwHXAP8GNhb\ntKiyJvVHtkGsg9SGsL8G/F1m9t36arWpmuuprsLQ614BvD4i7qEaJfiFiPhC4ZraIjMfrH3dRtUv\n1A99YluBLZn5n7X7X6MKZv3iQuBHtZ95P3g1cE9mbq9N0f0j8LOFa2qLzPxsZr40M4eBnVRLVfWL\nhyNiCUDt2tqPTOabOjmI9dsoAcDfApsy8xOlC2mXiHje2Jkltama1wB3lK2q9TLzg5l5SmaupGri\n/V5mXlq6rlaLiGNro75jiz1fQDWd0dNq0xVbImJVbdf59NdJGr9Gn0xL1vwUeHlEzK0tXn4+fXBy\nBkBEnFD7egrwRuBLZStqqUNzyjeA36xt/waTXLD+qJc4KiEivgQMA8dHxE+BK8eaXHtVRLwC+HVg\nQ61fKoEP9sGCtycCn4+IGVR/GHwlM68pXJNaZwnw9dp1ZAeAL2bmdwrX1C6/B3yxNk13D/DWwvW0\nRa1X6NXAb5WupV0y88aI+BrV1Nxzta9/U7aqtvmHiFhM9b7f1asnpYyXU4APA1dFxNuA+4CLJ3Us\nF3SVJEkqo5OnJiVJknqaQUySJKkQg5gkSVIhBjFJkqRCDGKSJEmFGMQkSZIKMYhJkiQV8v8BiZJr\nmeJryF8AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x1124f6da0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,1,figsize=(10,7))\n",
"chart = ax.plot(range(1,11),aics)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"次数=5(切片も含め推定パラメータ数k=6)の時が最もAICが低い。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"-----\n",
"**なぜAICが小さいモデルが良いモデルと言えるのか、なぜAICは予測の良さを図る指標と言えるか等、時間があるときに追記。久保本p77~p91参照。**"
]
}
],
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment