{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 機械学習のプロトタイプづくり\n",
"\n",
"機械学習の手法を実装するにあたって、アプローチはいくつもあり、その目的に依存する部分が多く、ただ一つの正確があるわけではない。我々の目的としては、機械学習の手法を効率よく「発想」の段階から「動くソフトウェア」の段階まで持っていくことができれば良い。なぜなら、色々な手法を試したり、調整したりする必要があるからである。実装の細かいところに神経を使わないといけないのであれば、バグは間違いなく続出するであろうし、本来の目的である「機械学習の理解」がおろそかにされてしまう可能性が高い。\n",
"\n",
"そのため、簡明でわかりやすく、また多様な手法でも使い勝手の良い実装方法を心がける。ここで紹介するアプローチは、__データ__、__モデル__、そして__アルゴリズム__という3つのオブジェクトを中心としている。その基本的な性質は下記の通りである。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- __データ(オブジェクト):__ 訓練、検証などに使うすべてのデータを格納するクラスオブジェクトである。\n",
"\n",
"- __モデル(オブジェクト):__ アルゴリズムの状態とデータに依存する関数をメソッドとして持ち、訓練中および訓練後の評価に使うクラスオブジェクトである。\n",
"\n",
"- __アルゴリズム(オブジェクト):__ 制御対象となるパラメータなどの「状態」を持つイテレータである。データとモデルに依存する関数をメソッドとして持ち、繰り返して更新するごとに呼び出す。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__目次__\n",
"\n",
"- データとモデルを定義してみる\n",
"\n",
"- 自明なアルゴリズムを作ってみる\n",
"\n",
"- 三者を結びつける\n",
"\n",
"- 最適化と学習\n",
"\n",
"___"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## データとモデルを定義してみる\n",
"\n",
"簡単な具体例を拠り所として、基本的なオブジェクトの定義を説明していく。古典的な線形モデル、つまり\n",
"\n",
"\\begin{align*}\n",
"y = \\langle w^{\\ast}, x \\rangle + \\varepsilon, \\quad \\mathbf{E}\\varepsilon = 0, \\mathbf{E}\\varepsilon^2 < \\infty, x \\in \\mathbb{R}^d\n",
"\\end{align*}\n",
"\n",
"を出発点とする。データ$(x_1, y_1),\\ldots,(x_n, y_n)$を所与として、未知の$(x,y)$に対して$\\langle \\widehat{w}, x \\rangle \\approx y$という近似がなるべく正確であるように、$\\widehat{w}$を決めたい。\n",
"\n",
"段取りとして、まず必要なのは、これら$n$個のデータ点を格納すること。`DataSet`をベースクラスとして以下のように定義する。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"class DataSet:\n",
" '''\n",
" Base class for data objects.\n",
" '''\n",
" \n",
" def __init__(self, paras=None):\n",
" self.X_tr = None\n",
" self.X_te = None\n",
" self.y_tr = None\n",
" self.y_te = None\n",
" self.paras = paras\n",
" \n",
" def init_tr(self, X, y):\n",
" self.X_tr = X\n",
" self.y_tr = y\n",
" self.n_tr = X.shape[0]\n",
" \n",
" def init_te(self, X, y):\n",
" self.X_te = X\n",
" self.y_te = y\n",
" self.n_te = X.shape[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"今の学習課題は線形モデルを前提とした回帰である。我々の使うモデルオブジェクトとして、`Model`と呼ぶベースクラスを作りたいところだが、この「線形モデル」と「回帰」のどこがモデルオブジェクトの範囲内なのか。以下のように決めておこう。\n",
"\n",
"- `Model`のサブクラスは、訓練および検証等の評価に使う各種の損失(ロス)関数(必要に応じてその勾配等の情報)を実装したメソッドを持つとする。\n",
"\n",
"このように決めておくとあとは簡単である。モデルオブジェクトのベースクラスとなる`Model`を以下のように定義する。"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"class Model:\n",
" '''\n",
" Base class for model objects.\n",
" '''\n",
"\n",
" def __init__(self, name=None):\n",
" self.name = name\n",
"\n",
" def l_imp(self, w=None, X=None, y=None, lamreg=None):\n",
" raise NotImplementedError\n",
" \n",
" def l_tr(self, w, data, n_idx=None, lamreg=None):\n",
" if n_idx is None:\n",
" return self.l_imp(w=w, X=data.X_tr,\n",
" y=data.y_tr,\n",
" lamreg=lamreg)\n",
" else:\n",
" return self.l_imp(w=w, X=data.X_tr[n_idx,:],\n",
" y=data.y_tr[n_idx,:],\n",
" lamreg=lamreg)\n",
" \n",
" def l_te(self, w, data, n_idx=None, lamreg=None):\n",
" if n_idx is None:\n",
" return self.l_imp(w=w, X=data.X_te,\n",
" y=data.y_te,\n",
" lamreg=lamreg)\n",
" else:\n",
" return self.l_imp(w=w, X=data.X_te[n_idx,:],\n",
" y=data.y_te[n_idx,:],\n",
" lamreg=lamreg)\n",
"\n",
" def g_imp(self, w=None, X=None, y=None, lamreg=None):\n",
" raise NotImplementedError\n",
" \n",
" def g_tr(self, w, data, n_idx=None, lamreg=None):\n",
" if n_idx is None:\n",
" return self.g_imp(w=w, X=data.X_tr,\n",
" y=data.y_tr,\n",
" lamreg=lamreg)\n",
" else:\n",
" return self.g_imp(w=w, X=data.X_tr[n_idx,:],\n",
" y=data.y_tr[n_idx,:],\n",
" lamreg=lamreg)\n",
" \n",
" def g_te(self, w, data, n_idx=None, lamreg=None):\n",
" if n_idx is None:\n",
" return self.g_imp(w=w, X=data.X_te,\n",
" y=data.y_te,\n",
" lamreg=lamreg)\n",
" else:\n",
" return self.g_imp(w=w, X=data.X_te[n_idx,:],\n",
" y=data.y_te[n_idx,:],\n",
" lamreg=lamreg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"ベースクラスなので`Model`では何も実装していないが、ロス関数にかかわるメソッドは`DataSet`を前提とした`data`と任意のパラメータ`w`を引数とする。\n",
"\n",
"回帰モデルを実装するためには、言うまでもなくロスや勾配の詳細を定めて、それらを実装する必要がある。どのロスを使うかは我々の自由である。典型例として、「最小二乗法」として古くから知られている回帰モデルを使うことにする。この場合、`w`が線形写像を定めるベクトル$w \\in \\mathbb{R}^d$に相当する。「最小」と呼ぶのは、最終的にこの誤差が小さくなるように最適化していくことを表わす。式で示すと以下のとおりである。\n",
"\n",
"\\begin{align*}\n",
"\\min_{w \\in \\mathbb{R}^d} \\frac{1}{n} \\sum_{i=1}^{n} (\\langle w, x_i \\rangle - y_i)^2 \\to \\hat{w}.\n",
"\\end{align*}\n",
"\n",
"この最適化を行うにあたって、ロスの取る値そのもの、あるいはそのロス関数の勾配ベクトルを$w$について求める。具体的には、下記を実装する必要がある。\n",
"\n",
"- `l_imp`: 任意のペア$(x,y)$に対して、$(y - \\langle w, x \\rangle)^2 / 2$を計算する。\n",
"\n",
"- `g_imp`: 任意のペア$(x,y)$に対して、$(y - \\langle w, x \\rangle)(-1)x \\in \\mathbb{R}^d$を計算する。\n",
"\n",
"これほど単純なモデルはほかにないので、実装することも大して苦労しない。`LinReg`をベースクラスとして、`LinearL2`をそのサブクラスとして定義する。"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"class LinReg(Model):\n",
" '''\n",
" General-purpose linear regression model.\n",
" No losses are implemented, just a predict()\n",
" method.\n",
" '''\n",
" \n",
" def __init__(self, data=None, name=None):\n",
" super(LinReg, self).__init__(name=name)\n",
" pass\n",
"\n",
" def predict(self, w, X):\n",
" '''\n",
" Predict real-valued response.\n",
" w is a (d x 1) array of weights.\n",
" X is a (k x d) matrix of k observations.\n",
" Returns array of shape (k x 1) of predicted values.\n",
" '''\n",
" return X.dot(w)\n",
"\n",
"\n",
"class LinearL2(LinReg):\n",
" '''\n",
" An orthodox linear regression model\n",
" using the l2 error; typically this\n",
" will be used for the classical least\n",
" squares regression.\n",
" '''\n",
" \n",
" def __init__(self, data=None):\n",
" super(LinearL2, self).__init__(data=data)\n",
"\n",
" \n",
" def l_imp(self, w, X, y, lamreg=None):\n",
" '''\n",
" Implementation of l2-loss under linear model\n",
" for regression.\n",
"\n",
" Input:\n",
" w is a (d x 1) matrix of weights.\n",
" X is a (k x numfeat) matrix of k observations.\n",
" y is a (k x 1) matrix of labels in {-1,1}.\n",
" lamreg is a regularization parameter (l2 penalty).\n",
"\n",
" Output:\n",
" A vector of length k with losses evaluated at k points.\n",
" '''\n",
" if lamreg is None:\n",
" return (y-self.predict(w=w,X=X))**2/2\n",
" else:\n",
" penalty = lamreg * np.linalg.norm(w)**2\n",
" return (y-self.predict(w=w,X=X))**2/2 + penalty\n",
" \n",
" \n",
" def g_imp(self, w, X, y, lamreg=None):\n",
" '''\n",
" Implementation of the gradient of the l2-loss\n",
" under a linear regression model.\n",
"\n",
" Input:\n",
" w is a (d x 1) matrix of weights.\n",
" X is a (k x numfeat) matrix of k observations.\n",
" y is a (k x 1) matrix of labels in {-1,1}.\n",
"\n",
" Output:\n",
" A (k x numfeat) matrix of gradients evaluated\n",
" at k points.\n",
" '''\n",
" if lamreg is None:\n",
" return (y-self.predict(w=w,X=X))*(-1)*X\n",
" else:\n",
" penalty = lamreg*2*w.T\n",
" return (y-self.predict(w=w,X=X))*(-1)*X + penalty"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上記のオブジェクトがあれば、データと評価の仕組みなどができているので、学習アルゴリズムに渡す情報の準備はできている。あとはアルゴリズムそのもののクラスを用意するだけだが、その前に簡単な使用例を示しておく。"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"===============\n",
"l_tr: 63.67480525704821 l_te: 72.52161001938249\n",
"g_tr:\n",
" [-1.1240842 -6.437987 -1.79680855 -7.2479921 -4.13374471]\n",
"g_te:\n",
" [ 1.23114245 -6.7552612 -1.8668131 -8.89835267 -4.93590447]\n",
"===============\n",
"l_tr: 35.945556539267216 l_te: 40.85252548386378\n",
"g_tr:\n",
" [-0.84580072 -4.82831557 -1.36736611 -5.42541324 -3.08957585]\n",
"g_te:\n",
" [ 0.94982208 -5.06642986 -1.37297361 -6.65188944 -3.6999843 ]\n",
"===============\n",
"l_tr: 16.167935313729377 l_te: 18.30555145320523\n",
"g_tr:\n",
" [-0.56751724 -3.21864415 -0.93792366 -3.60283438 -2.045407 ]\n",
"g_te:\n",
" [ 0.66850172 -3.37759852 -0.87913411 -4.40542621 -2.46406413]\n",
"===============\n",
"l_tr: 4.3419415804346695 l_te: 4.880687927406857\n",
"g_tr:\n",
" [-0.28923376 -1.60897273 -0.50848121 -1.78025552 -1.00123814]\n",
"g_te:\n",
" [ 0.38718136 -1.68876719 -0.38529462 -2.15896298 -1.22814396]\n",
"===============\n",
"l_tr: 0.46757533938310303 l_te: 0.5779349064686516\n",
"g_tr:\n",
" [-0.01095028 0.00069869 -0.07903877 0.04232334 0.04293072]\n",
"g_te:\n",
" [1.05860995e-01 6.41498784e-05 1.08544880e-01 8.75002561e-02\n",
" 7.77620457e-03]\n",
"===============\n"
]
}
],
"source": [
"# Prepare simulated data.\n",
"d = 5\n",
"n_tr = 100\n",
"n_te = 100\n",
"w_star = np.arange(d).reshape((d,1)) + 1.0\n",
"\n",
"X_tr = np.random.normal(size=(n_tr,d))\n",
"noise_tr = np.random.normal(size=(n_tr,1))\n",
"y_tr = X_tr.dot(w_star) + noise_tr\n",
"\n",
"X_te = np.random.normal(size=(n_te,d))\n",
"noise_te = np.random.normal(size=(n_te,1))\n",
"y_te = X_te.dot(w_star) + noise_te\n",
"\n",
"# Prepare data object.\n",
"data = DataSet()\n",
"data.init_tr(X=X_tr, y=y_tr)\n",
"data.init_te(X=X_te, y=y_te)\n",
"\n",
"# Prepare model object.\n",
"mod = LinearL2()\n",
"\n",
"# Prepare random initial value, and take closer to solution.\n",
"w_init = np.random.uniform(low=-5.0, high=5.0, size=(d,1))\n",
"weight_balance = np.linspace(0.0, 1.0, 5)\n",
"print(\"===============\")\n",
"for wt in weight_balance:\n",
" w_est = wt * w_star + (1-wt) * w_init\n",
" print(\"l_tr:\", np.mean(mod.l_tr(w_est, data=data)),\n",
" \"l_te:\", np.mean(mod.l_te(w_est, data=data)))\n",
" print(\"g_tr:\\n\", np.mean(mod.g_tr(w_est, data=data), axis=0))\n",
" print(\"g_te:\\n\", np.mean(mod.g_te(w_est, data=data), axis=0))\n",
" print(\"===============\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"明らかなように、パラメータを$w^{\\ast}$に近づけていくにつれて、ロスの値が小さくなり、その点での勾配も平坦になっていく。次は学習アルゴリズムについて考える。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 自明なアルゴリズムを実装する\n",
"\n",
"まずはきわめて単純なアルゴリズムを作る。特別な働きは何もないが、自身の「状態」と「かかった反復回数」を管理する能力、そして条件に応じて終了する能力は持つ。"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"class Algo_trivial:\n",
" '''\n",
" An iterator which does nothing special at all, but makes\n",
" a trivial update to the initial parameters given, and\n",
" keeps record of the number of iterations made.\n",
" '''\n",
"\n",
" def __init__(self, w_init, t_max):\n",
" \n",
" self.w = np.copy(w_init)\n",
" self.t_max = t_max\n",
"\n",
"\n",
" def __iter__(self):\n",
"\n",
" self.t = 0\n",
" print(\"(__iter__): t =\", self.t)\n",
"\n",
" return self\n",
" \n",
"\n",
" def __next__(self):\n",
" \n",
" # Condition for stopping.\n",
" if self.t >= self.t_max:\n",
" print(\"--- Condition reached! ---\")\n",
" raise StopIteration\n",
"\n",
" print(\"(__next__): t =\", self.t)\n",
" self.w += 5\n",
" self.t += 1\n",
" # Note: __next__ does not need to return anything.\n",
" \n",
"\n",
" def __str__(self):\n",
"\n",
" out = \"State of w:\" + \"\\n\" + \" \" + str(self.w)\n",
" return out\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"その振る舞いを見てみよう。"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Printing docstring:\n",
"\n",
" An iterator which does nothing special at all, but makes\n",
" a trivial update to the initial parameters given, and\n",
" keeps record of the number of iterations made.\n",
" \n",
"(__iter__): t = 0\n",
"(__next__): t = 0\n",
"(__next__): t = 1\n",
"(__next__): t = 2\n",
"(__next__): t = 3\n",
"(__next__): t = 4\n",
"(__next__): t = 5\n",
"(__next__): t = 6\n",
"(__next__): t = 7\n",
"(__next__): t = 8\n",
"(__next__): t = 9\n",
"--- Condition reached! ---\n"
]
}
],
"source": [
"\n",
"import numpy as np\n",
"\n",
"al = Algo_trivial(w_init=np.array([1,2,3,4,5]), t_max=10)\n",
"\n",
"print(\"Printing docstring:\")\n",
"print(al.__doc__) # check the docstring\n",
"\n",
"for onestep in al:\n",
" pass # do nothing special\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"要点:\n",
"\n",
" - `StopIteration`という例外を挙げてから、ただちにループから脱出する。\n",
" - 0番目のステップでは、`iter`も`next`も呼び出される。\n",
"\n",
"上記の`for`ループと同じ働きは手動でもできる。"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(__iter__): t = 0\n",
"(__next__): t = 0\n",
"(__next__): t = 1\n",
"(__next__): t = 2\n",
"(__next__): t = 3\n"
]
}
],
"source": [
"iter(al)\n",
"next(al)\n",
"next(al)\n",
"next(al)\n",
"next(al) # and so on..."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(__iter__): t = 0\n",
"(__next__): t = 0\n",
"State of w:\n",
" [ 6 7 8 9 10]\n",
"(__next__): t = 1\n",
"State of w:\n",
" [11 12 13 14 15]\n",
"(__next__): t = 2\n",
"State of w:\n",
" [16 17 18 19 20]\n",
"(__next__): t = 3\n",
"State of w:\n",
" [21 22 23 24 25]\n",
"(__next__): t = 4\n",
"State of w:\n",
" [26 27 28 29 30]\n",
"(__next__): t = 5\n",
"State of w:\n",
" [31 32 33 34 35]\n",
"(__next__): t = 6\n",
"State of w:\n",
" [36 37 38 39 40]\n",
"(__next__): t = 7\n",
"State of w:\n",
" [41 42 43 44 45]\n",
"(__next__): t = 8\n",
"State of w:\n",
" [46 47 48 49 50]\n",
"(__next__): t = 9\n",
"State of w:\n",
" [51 52 53 54 55]\n",
"--- Condition reached! ---\n",
"One last check after exiting:\n",
"State of w:\n",
" [51 52 53 54 55]\n"
]
}
],
"source": [
"al = Algo_trivial(w_init=np.array([1,2,3,4,5]), t_max=10)\n",
"\n",
"for onestep in al:\n",
" print(al) # useful for monitoring state.\n",
"\n",
"print(\"One last check after exiting:\")\n",
"print(al) # ensure that no additional changes were made."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Exercises:__\n",
"\n",
"0. `Algo_trivial`を反復するごとに、`w`がどのように更新されているか説明すること。この算法を「自明」と呼んでいるのは、データやモデルたるものには一切依存しない更新則だからである。\n",
"\n",
"0. 反復するごとに、`w`の全要素が倍増されるように`Algo_trivial`を改造してみること。\n",
"\n",
"0. さらに、初期値が$w=(0,1,2,3,4)$で、何回か倍増を繰り返したあとの状態が$w=(0, 16, 32, 48, 64) = (0 \\times 2^4, 1 \\times 2^4, 2 \\times 2^4, 3 \\times 2^4, 4 \\times 2^4)$となるように、`al`を初期化する際に使うパラメータを設定すること。\n",
"\n",
"0. `al`を初期化しなおすこと無く、`for`ループを複数回走らせると、`w`の状態がどうなっていくか。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 三者を結びつける\n",
"\n",
"先ほどの線形回帰の具体例に戻って、その実装を完成させていこう。有用なアルゴリズムを作るには、上手に観測データに応じたパラメータ更新則が不可欠である。ここで、経験期待損失最小化(empirical risk minimization, ERM)を、最急降下法で実装した学習アルゴリズムにする。\n",
"\n",
"要点は以下のとおりである。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- データ:入力ベクトルと実数値の応答$x \\in \\mathbb{R}^{d}$, $y \\in \\mathbb{R}$.データセットの全体は$\\{(x_{1},y_{1}),\\ldots,(x_{n},y_{n})\\}$とする。\n",
"\n",
"- モデル:二乗誤差を用いた線形回帰モデル。つまり、$y = \\langle w^{\\ast}, x\\rangle + \\varepsilon$と仮定している。ロス関数は$L(w;x,y) = (y - \\langle w, x\\rangle)^{2}/2$で、勾配が$\\nabla L(w;x,y) = -(y-\\langle w, x\\rangle)x$である。\n",
"\n",
"- アルゴリズム:最急降下法を使って、固定したステップサイズ$\\alpha > 0$を用いる。数式で表わすと、$z_{i}=(x_{i},y_{i})$として以下の通りである。\n",
"\n",
"\\begin{align*}\n",
"w_{(t+1)} \\gets w_{(t)} - \\alpha \\frac{1}{n} \\sum_{i=1}^{n}\\nabla L(w_{(t)}; z_{i}).\n",
"\\end{align*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"データもモデルもできているので、あとは先ほどのイテレータの準備を踏まえて、アルゴリズムのオブジェクトを定義することである。`Algo_SimpleGD`と名づけて、以下の通りである。"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import numpy as np\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"\n",
"class Algo_SimpleGD:\n",
" '''\n",
" Iterator which implements a line-search steepest descent method,\n",
" using the sample mean estimate of the gradient.\n",
" '''\n",
" \n",
" def __init__(self, w_init, t_max, step, store):\n",
" self.w = w_init\n",
" self.t = None\n",
" self.t_max = t_max\n",
" self.step = step\n",
" self.store = store\n",
" \n",
" # Keep record of all updates (optional).\n",
" if self.store:\n",
" self.wstore = np.zeros((self.w.size,t_max+1), dtype=np.float32)\n",
" self.wstore[:,0] = self.w.flatten()\n",
" else:\n",
" self.wstore = None\n",
" \n",
" def __iter__(self):\n",
" self.t = 0\n",
" return self\n",
" \n",
" def __next__(self):\n",
" if self.t >= self.t_max:\n",
" raise StopIteration\n",
" \n",
" def newdir(self, model, data):\n",
" return (-1) * np.mean(model.g_tr(w=self.w, data=data), axis=0, keepdims=True)\n",
"\n",
" def update(self, model, data):\n",
" \n",
" # Parameter update.\n",
" stepsize = self.step(self.t)\n",
" newdir = self.newdir(model=model, data=data)\n",
" self.w = self.w + stepsize * np.transpose(newdir)\n",
" \n",
" # Monitor update.\n",
" self.t += 1\n",
" \n",
" # Keep record of all updates (optional).\n",
" if self.store:\n",
" self.wstore[:,self.t] = self.w.flatten()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"念頭におくべき要点:\n",
"\n",
"- `w_init`が初期値の$w_{(0)}$で、アルゴリズムの初期化の際に必要である。\n",
"\n",
"- `newdir`は更新方向を返す。\n",
"\n",
"- `step`はコールバック関数で、たとえばステップ数の`t`に応じてステップサイズを変えるときなどに活用できる。\n",
"\n",
"- `update`はパラメータの更新を実際に行うメソッドである。そのほかに計算資源のコスト等を管理するメソッドを呼び出す。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"まずは試しに、乱数を発生させてデータを生成し、自前の学習機を動かしてみる。"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"X_tr: (100, 2) ... sumcheck = -9.94677583241537\n",
"y_tr: (100, 1) ... sumcheck = -40.20352170786684\n",
"X_te: (100, 2) ... sumcheck = -1.089890767972685\n",
"y_te: (100, 1) ... sumcheck = 16.95676631804479\n"
]
}
],
"source": [
"# Generate data.\n",
"data = DataSet()\n",
"n = 100\n",
"d = 2\n",
"w_star = math.pi * np.ones((d,1), dtype=np.float32)\n",
"X = np.random.standard_normal(size=(n*2,d))\n",
"epsilon = np.random.standard_t(df=3, size=(n*2,1))\n",
"y = X.dot(w_star) + epsilon\n",
"data.init_tr(X=X[0:n,:], y=y[0:n,:]) # former half for training\n",
"data.init_te(X=X[n:,:], y=y[n:,:]) # latter half for testing\n",
"\n",
"print(\"X_tr:\", data.X_tr.shape, \"... sumcheck =\", np.sum(data.X_tr))\n",
"print(\"y_tr:\", data.y_tr.shape, \"... sumcheck =\", np.sum(data.y_tr))\n",
"print(\"X_te:\", data.X_te.shape, \"... sumcheck =\", np.sum(data.X_te))\n",
"print(\"y_te:\", data.y_te.shape, \"... sumcheck =\", np.sum(data.y_te))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Initialize model.\n",
"mod = LinearL2()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Initialize learning algorithm.\n",
"\n",
"w_init = np.random.uniform(size=(d,1))\n",
"\n",
"def alpha_fixed(t, val):\n",
" return val\n",
"\n",
"def make_step(u):\n",
" def mystep(t):\n",
" return alpha_fixed(t=t, val=u)\n",
" return mystep\n",
"\n",
"al = Algo_SimpleGD(w_init=w_init,\n",
" t_max=15,\n",
" step=make_step(0.15),\n",
" store=True)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"State of algorithm after completion:\n",
"t = 15\n",
"w = [[2.97596384]\n",
" [2.74678402]]\n"
]
}
],
"source": [
"# Iterate the learning algorithm.\n",
"for onestep in al:\n",
" al.update(model=mod, data=data)\n",
"\n",
"print(\"State of algorithm after completion:\")\n",
"print(\"t =\", al.t)\n",
"print(\"w =\", al.w)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"比較対象として、線形モデルの下で二乗誤差の和を最小にする解が解析的に求まることはわかっているので、それも合わせて計算しておく(`w_OLS`と書く)。"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"w_OLS = np.linalg.lstsq(a=data.X_tr, b=data.y_tr, rcond=None)[0]"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAGrCAYAAABg2IjeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VIW9xvH3RwTCDgJlkVVQi4CNENxBFLfeti6Fi4obaKHuXmup1j5eFa29lVYpLiiCoAiCoKhFsV7RusM11CCbYlHUsJTFyiapEH73jzNJE5KQCZnMOTPz/TzPPEnOOTnznplk3jnLnGPuLgAAoqBO2AEAAChGKQEAIoNSAgBEBqUEAIgMSgkAEBmUEgAgMiglVMrMssxsh5l1CjtLoljgSTP7xszeq8X7mWRmt9bSvAvMbGBtzDuqzGygmS0POwdqH6WURmIFUnzba2a7Sv18UXXn5+5F7t7Y3b+sYa53zGx4TeaRQAMlnSypvbufUFt34u4/c/d7amv+B8rMmpnZODP7wsx2xr7ONrN+sfEHmZnHxu0ws81m9pqZ/WcSMxZn6FI8zN3/6u49k5UB4aGU0kisQBq7e2NJX0r6Salh0/ed3swOSn7K6jGzOmaWyL/TzpI+d/dvEzjPlGBm2ZLekPR9Sf8hqamkIyU9I+mH+0zeM/Z39H1JT0maYGa/ifN+Iv93hQhzd25peJO0RtJp+wy7W9IsSU9L2i5puKTjJS2U9I2k9ZLGS6obm/4gSS6pS+znbEn3SfpK0j8kPSwpu9T8fyopX9I2SX+XdIak30sqklQoaYekcbFpT5KUJ2mrpP+TdGyp+bwj6S5J70vaJenXkhbtsyw3S5pTybJ3kDRP0teSPpV0eWz4qFiOoliW2yr43cMUvHBvkbRZ0jRJzSq5nzqxx2tjbDk+knRkbNxTku6IfX9a7Pn4taRNktZJ+omkH8fyfS3pVxU8T7Njz1OepN6lxhdIGlgqw62SVsfyzpTUopK8V0paK6nBfv5uyjznpYZfEHsumlfyewWSRktaKum72LCekt6M/W0tlfSjfZ7j4aV+/pmkv8a+fy+WYWfseRpc/Bjuc3+/iM13q4K/6fqlxp8taUnsvt+R1Cvs/0lu8d1CD8Ctlp7Yykvpu9gLYh1JDST1k3Rs7MXoUEmrJF0bm37fUnpQ0lxJLRS8y35Z0l2xcSfEXgAGxebdUdIRsXH7vgC1ir2QXBi7j4sVlECLUtOvkdRDUl1JTWLzPqzUPJZKOqeSZX9X0gMKSrRP7MX65Ni4khe/Sn738Ngy1JP0vdi8/lDJtD9SUKjNYst8pKS2sXH7ltIeSb+JLc9VCorsKUmNJR2loCw7lXqedks6Lzb9LQpK/qDY+NKl9MtYxkNiyztZ0rRK8s6RNKmKv5vKSilb0l5Jp1fyewWSFit4Q9Ag9vh9LulXsWU4TUHBdK/kb6J0KZXLoIpLaaGktpJaKvi7/VlsXD8Fb5r6ScqSdLmC0q4X9v8lt6pvbL7LPO+4+5/dfa+773L3D9x9kbvvcffPJE1UsM+ljNgmtJ9J+i93/6e7b5P0OwXvoCXpCkmPufuC2Ly/cvdPKsnwE0nL3f3p2P0+JekzBS/yxR5395XuvtvdtytYa7g4liVHUjsFpbhvzq6SjpF0i7sXuvvfJE2RdEk8D467r4otw3fuvlHS/RU9HjG7FZTz92O/u8LdN1QybaGk/3H33QrWZlpLut/dd7j7R5I+UVBOxRa5+9zY9GNj99Ovgvn+XNKt7r7W3Qsl3SFpaCWbPFtJKslnZrmxAz62VXUQQWzeX0s6eD+T/cndC9x9l6QTFRTT2Nhz+Jqk+fr330sijHP3De6+RcGacU5s+ChJD8f+tovc/fHY8IoeP0QMpZR5vir9g5l938xeMrMNZrZN0hgFL177aiupvqQlsReybxS8EHwvNr6jgnej8Wgv6Yt9hn2h4N1+hTklPSGp+GCNiyXNir1gVzTvze6+cz/zrpSZtTWzZ8xsbezxmKqKHw+5+6uSHpE0QdI/zOwRM2tSyaw3u3tR7Ptdsa//KDV+l4K1pmIlyx/7vbWxZdtXJ0l/LvWcLFWwlvG9CqbdoqDMi+eb5+7NJQ1V8NxWKrY/6mAFxVSZ0s9Ze0lfunvpMz7H/TzEqfQbgG/178evs6Sbix+T2OPSLsH3jVpCKWWefU8L/6ikZQo2qzSV9N+SrILf+4eCTX9HuHvz2K2ZuzeLjf9KUrc473OdgheO0jopeOGt8Hfc/R1JMrMTFWz2m1bJfa2T1MrMGu1n3vvze0n/UrAPp6mC/W4VPR7Fuca5ex9JvRRsvvtFnPdTlY7F38TWeg5RsGz7KlCwSa15qVt2JWtsCySdZWYNDyDPuQoelw/2M03p52ydpI5mVvqxK/087JRUOkfbSuZzIL6SdOc+j0lDd3+mhvNFElBKaKJg/85OM+uhYHNQObF365MkjTOz1rHP+3QwszNik0yW9DMzOyV2xFwHMzsiNu4fCvZXFZsnqaeZnR87/HeYpO6qYHPcPqYpWCvZ6e4LK8n5uYIDA+4xs/qxTX0jJJU7+rASTRS8YG41s44K9tlUyMyOid0Oiv3OdwoOokiEY8zsHDOrG8uwXRUXwiMKlrVTLNP3zOzsSuY5RcH+tefMrGfsc2gNJOVWFsLMWprZJQr20f3O3b+JM/97Cvaj3WRmdc3sVAVH/BUXQ76kwWbWwMwOV7DfR1LJ39oWlf2bqY6Jkq4xs36xv9PGZvaTfd6oIKIoJdwk6TIFL3qPKjjqa3/TfqFg5/5WSa8qOFpN7v6epJEKjkbbquAItuJ3++MkXRjblHKfu29ScHTUzQpefG6U9GN339+mIUl6UsEaSWVrScXOj+XaoGDn/q3u/kYVv1PsdgX7pLZKelHSs/uZtrmCMv5GwYEZ6xXsg0qEuQo2U36tYHl+6u57KpjuPkmvSFpgZtsVlEGF+05i+3pOVrD/ar6CoyQ/lvQDld/Xs9zMdig4OnCEpOvcfUy84d39Xwr2HZ6joAjHSxrm7qtik/xBwRrRRkmPKzjoo7TbJc2I/c38NN77jd33IgUHk0yQ9E8FB0FcXJ15IDxWdpMv8G9mVk/BJptD3L2iTUfJztNIwYtYr9gaUVoys7sldXD34WFnAZKNNSXsTy8FO5A3hh0k5hpJ76ZzIQGZjk9eo0Jmdr6khxR8qLOizUbJzlOg4BDsc8LOAqD2sPkOABAZbL4DAERGrWy+a9WqlXfp0qU2Zg0ASEGLFy/e7O6tq5quVkqpS5cuysvLq41ZAwBSkJntexaXCrH5DgAQGZQSACAyKCUAQGTwOSUAqKHdu3eroKBAhYWFYUcJXXZ2tjp06KC6dese0O9TSgBQQwUFBWrSpIm6dOmisidGzyzuri1btqigoEBdu3Y9oHmw+Q4AaqiwsFAtW7bM6EKSJDNTy5Yta7TGSCkBQAJUu5B27pQuukj69tvaCRSSmhYzpQQAYXj/fWnGjOArSlBKABCG114r+zUCpk6dqnXrwr1KDaUEAGF46aXg67x54eYo5UBKac+exF5EgFICgNp2zjmSWdnbqthFeFetKj/unOpdoeXee+/V+PHjJUk33nijTj31VEnSggULdPHF5S+6W1RUpOHDh6tXr17q3bu37r//fs2ZM0d5eXm66KKLlJOTo127dmnMmDHq16+fevXqpVGjRqn4qhIDBw7UrbfeqpNPPll/+tOfavDAlEcpAUBtu+ceqVMnKTv738O++67sVykY37lzMH01DBgwQG+//bYkKS8vTzt27NDu3bv1zjvvqH///uWmz8/P19q1a7Vs2TItXbpUI0aM0JAhQ5Sbm6vp06crPz9fDRo00LXXXqsPPvhAy5Yt065duzSv1FrdN998ozfffFM33XRTtbJWhVICgNrWs6e0YoV09tlSw4YVT9OwYbCGtHx5MH019O3bV4sXL9b27dtVv359HX/88crLy9Pbb79dYSkdeuih+uyzz3TdddfplVdeUdOmTSuc7xtvvKFjjz1WvXv31uuvv67ly5eXjDv//POrlTFelBIAJEOjRtKsWdKYMVKDBmXHNWgQDJ85M5iumurWrasuXbpoypQpOuGEE9S/f3+98cYbWr16tXr06FFu+hYtWmjJkiUaOHCgHnroIf3sZz8rN01hYaGuvvpqzZkzR0uXLtXIkSPLfP6o0QHkjAelBADJtHq1VFQU7Dtq2DD4WlQkffZZjWY7YMAA/eEPf9CAAQPUv39/PfLII8rJyanwc0ObN2/W3r17NXjwYN11113629/+Jklq0qSJtm/fLkklBdSqVSvt2LFDc+bMqVG+eFFKAJAs69dLkyYF33fqJE2fLnXsGPz82GPShg0HPOv+/ftr/fr1Ov7449WmTRtlZ2dXuOlOktauXauBAwcqJydHw4cP1+9+9ztJ0vDhw3XllVcqJydH9evX18iRI9W7d2+de+656tev3wFnqw4rPpoikXJzc52L/AHIFCtXrqxwM1k5V18tTZggnX++NHlysKlu507p8sulZ54Jxj/0UO0HrmUVPR5mttjdc6v6XdaUACBZvv02KKPS+46K9zVNnhwUVIbjLOEAkCxTp1Y+7vLLg1uCHXvssfrXv/5VZti0adPUu3fvhN9XIlBKAJDGFi1aFHaEamHzHQAgMiglAEBkUEoAgMiglAAAkcGBDgCQJEc/erTyN+RXOj6nbY4+/PmHSUwUPawpAUCSHN/heNXLqlfhuHpZ9XRChxOSlmXr1q269NJL1a1bN3Xr1k2XXnqptm7dKklas2aNevXqVe53Fi5cqGOPPVY5OTnq0aOH7rjjjoTnopQAIEluG3Cb6ljFL7tZlqXbTr4taVmuuOIKHXrooVq9erVWr16trl27Vnhi1tIuu+wyTZw4Ufn5+Vq2bJmGDh2a8FyUEgAkSbsm7TQiZ0S5taV6WfU0ImeE2jZue0Dzre5F/v7+979r8eLFuu22f5fgf//3fysvL0+rV6+u9H42btyodu3aSZKysrJ05JFHHlDe/aGUACCJKlpbqulaUnUv8rdixQrl5OQoKyvr3xmyspSTk1Pmmkn7uvHGG3XEEUfovPPO06OPPlrmUhaJQikBQBLtu7ZU07UkqfoX+XP3Ci9pUdnwYsVrU2eccYZmzJihs84664AzV4ZSAoAkK722lIh9SdW9yF/Pnj314Ycfau/evSXD9u7dqyVLllR5tvNu3brpqquu0oIFC7RkyRJt2bKlRtn3RSkBQJIVry3VsTo1XksqVp2L/HXv3l1HH3207r777pJhd999t/r06aPu3btXeh8vvfSSii939OmnnyorK0vNmzevcfbSKCUACMFtA27TSZ1OStgRd9W5yJ8kTZ48WatWrVL37t3VrVs3rVq1SpMnTy4Z/8knn6hDhw4lt9mzZ2vatGk64ogjlJOTo0suuUTTp08vs18qEbjIHwDUUNwX+csQXOQPAJAWOM0QAKSxtLzIn5k1lzRJUi9JLulyd3+/NoMBAGou1S7yF++a0p8kveLuQ8ysnqSGtZgJAJChqiwlM2sqaYCk4ZLk7t9J+q52YwEAMlE8BzocKmmTpClm9qGZTTKzRvtOZGajzCzPzPI2bdqU8KAAgPQXTykdJKmPpAnufrSknZJu2Xcid5/o7rnuntu6desExwSA9LBjh3T77VLr1lKdOsHX228PhiO+UiqQVODuxXvL5igoKQBANezYIR13nHTvvdLmzZJ78PXee4PhFFMcpeTuGyR9ZWZHxAYNkrSiVlMBQBoaO1ZavVra9+TahYXB8LFjD2y+1b10RZTF++HZ6yRNN7OPJOVIuqf2IgFAenr44fKFVKywUJow4cDmW91LV0RZXIeEu3u+pCpPDwEAqFxVJ9Q+0BNu73vpij59+pRcuqJ4DSpVcEYHAEiSli2DfUj7G38g9r10xVFHHbXfS1dEGee+A4AkufpqKTu74nHZ2dJVVx34vKtz6Yooo5QAIElGj5a6dStfTNnZwfDRow983tW9dEVUsfkOAJKkcWNp4cLgKLsJE4J9SC1bBmtIo0cH4w/UoEGDtHv37pKfV61alYDEyUcpAUASNW4s3XlncEN5bL4DAEQGpQQACVAbV/FORTV9HCglAKih7OxsbdmyJeOLyd21ZcsWZVd2iGEc2KcEADXUoUMHFRQUiCskBAXdoUOHA/59SgkAaqhu3brq2rVr2DHSApvvAACRQSkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFBKQEAIoNSAgBEBqUEAIgMSgkAEBmUEgAgMiglAEBkUEoAgMiglAAAkUEpAQAig1ICAEQGpQQAiAxKCQAQGZQSACAyKCUAQGRQSgCAyKCUAACRQSkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFxUDwTmdkaSdslFUna4+65tRkKAJCZ4iqlmFPcfXOtJQEAZDw23wEAIiPeUnJJr5rZYjMbVdEEZjbKzPLMLG/Tpk2JSwgAyBjxltKJ7t5H0g8lXWNmA/adwN0nunuuu+e2bt06oSEBAJkhrlJy93WxrxslzZV0TG2GAgBkpipLycwamVmT4u8lnSFpWW0HAwBknniOvmsjaa6ZFU8/w91fqdVUAICMVGUpuftnkn6QhCwAgAzHIeEAgMiglAAAkUEpAQAig1ICAEQGpQQAiAxKCQAQGZQSACAyKCUAQGRQSgCAyKCUAACRQSkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFBKQEAIoNSAgBEBqUEAIgMSgkAEBmUEgAgMiglAEBkUEoAgMiglAAAkUEpAQAig1ICAEQGpQQAiAxKCQAQGZQSACAyKCUAQGRQSgCAyKCUAACRQSkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFBKQEAIoNSAgBERtylZGZZZvahmc2rzUAAgMxVnTWlGyStrK0gAADEVUpm1kHSjyRNqt04AIBMFu+a0jhJv5K0t7IJzGyUmeWZWd6mTZsSEg4AkFmqLCUz+7Gkje6+eH/TuftEd89199zWrVsnLCAAIHPEs6Z0oqSzzWyNpJmSTjWzp2o1FQAgI1VZSu7+a3fv4O5dJF0g6XV3v7jWkwEAMg6fUwIARMZB1ZnY3f8q6a+1kgQAkPFYUwIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFBKQEAIoNSAgBEBqUEAIgMSgkAEBmUEgAgMiglAEBkUEoAgMiglAAAkUEpAQAig1ICAEQGpQQAiAxKCQAQGZQSACAyKCUAQGRQSgCAyKCUAACRQSkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFBKQEAIoNSAgBEBqUE1JadO6WLLpK+/TbsJEDKoJSA2vL++9KMGcFXAHGhlIDa8tprZb8CqBKlBNSWl14Kvs6bF24OIIUcFHYAIB0cfXMz5TfcVnbgkNhNy5RzpenDR0uNO/ts6YUXkhcQSBGsKQEJcHyvH6renorH1dsjnfBV7IfsbKlzZ+mee5KWDUgllBKQALede7/q1M+ucFyWS7e9JalhQ+mcc6Tly6WePZMbEEgRlBKQAO2atNOIo0eoXla9MsPr7ZFGfCi1LWogjRkjzZwpNWoUUkog+iglIEFuG3Cb6ljZf6mStaSiIumzz8IJBqQQSglIkHZN2mnEYUNL9i3V2yONaHuW2h7cKRjw2GPShg3hBQRSQJWlZGbZZvZ/ZrbEzJab2Z3JCAakotvekup48H1WdrZuGz5FWrFCOvdcafdu6a67wg0IRFw8a0r/knSqu/9AUo6ks8zsuNqNBaSmdjtMI5oPVB2roxE5l6tt47bBPqRZs6TJk4NTDwGolLl7/BObNZT0jqSr3H1RZdPl5uZ6Xl5eAuIBqWf99vW64NkLNGvIrKCUAMjMFrt7blXTxbVPycyyzCxf0kZJ/1tRIZnZKDPLM7O8TZs2VT8xkCbaNWmnN4e/SSEBByCuUnL3InfPkdRB0jFm1quCaSa6e66757Zu3TrROQEAGaBaR9+5+zeS/irprFpJAwDIaPEcfdfazJrHvm8g6TRJH9d2MABA5onnhKztJD1hZlkKSuwZd+e0xwCAhKuylNz9I0lHJyELACDDcUYHAEBkUEoAgMiglAAAkUEpAdWwY4d0++1S69ZSnTrB19tvD4YDqDkuhw7EaccO6bjjpNWrpcLCYNjmzdK990rPPistXCg1bhxuRiDVsaYExGns2LKFVKywMBg+dmw4uYB0QikBcXr44fKFVKywUJowIbl5gHREKQFx2rKlZuMBVI1SAuLUsmXNxgOoGqUExOnqq6Xs7IqvP5adLV11VZIDAWmIUgLiNHLkN6pb9yvVr7+3zPDsbKlbN2n06JCCAWmEUgLi8MUXX2jQoGPVr991uvnmOmU+p/SrX3E4OJAofE4JiEOLFi104YUXqlu3brrkEunOO8NOBKQnSgmoQl5enq6++mq99dZbql+/fthxgLTG5jugEu6u++67TyeccIL69u2r7OxsmVnYsYC0RikBlTAzHXnkkWrbtq0uu+yysOMAGYHNd0AlZs2apRUrVmjx4sVq1apV2HGAjEApAfvYsmWLrrnmGs2aNUtLly5V69atw44EZAw23wGl7NmzR7fccotmzZqlk046Sb169Qo7EpBRWFMCSsnKytLZZ58td9epp54adhwg41BKQMyXX36pq6++WoMGDdKjjz6qPXv2hB0JyDhsvkPGKyoq0oMPPqiePXvqrbfe0hVXXKGsrCw+kwSEgFJCxps2bZpuueUW7dixQyNHjlTTpk3DjgRkLEoJGW/w4ME666yzdMghh+j6668POw6Q0dinhIzl7po9e7YefPBBLViwQMuXL1fnzp3DjgVkNNaUkJFWrlyp008/Xeeff76GDRumunXrKicnJ+xYQMajlJBRtm/frtGjR+uoo47SggUL1L59ew0fPjzsWABiKCVknEaNGpUc7j169GhlZ2eHnAhAMUoJGaVhw4bq3bu3brjhBrVv314jR44MOxKAUiglZIxXXnlFOTk5+vbbbzVu3Dj9+c9/VqNGjcKOBaAUSglp76OPPtKZZ56pH/7wh9q7d6+GDRsmSerTp0/IyQDsi1JC2lq3bp2uuOIK5eTk6NVXX5Uk3XXXXcrKygo5GYDKUEpIW++//76WLFkid5ck9e3bV+edd17IqQDsD6WEtDV48GD17dtX9erVkyT99re/5XLmQMRRSkhLX375pUaNGqULLrhA+fn5Ov3003XGGWeEHQtAFSglpJX169fr+uuv12GHHab169frlFNOUY8ePTRz5kzWkoAUwLnvkBa2bNmi3//+93rwwQe1a9cuZWVlaezYsSXjDz744BDTAYgXpYSUtnXrVt133326//77tX379pLhV155pb7//e+HmAzAgWDzHVLa119/LUnavXt3ybBmzZrp9ttvDysSgBqglJDSunbtqv/8z/8sU0q/+c1v1Lp16xBTAThQbL5DSluxYoWefvppzZgxQx999JFmzJih6667LuxYAA4QpYSUtHjxYt1zzz167rnnNHv2bA0ZMkSnnHKKjj76aM76DaQwSgkp5a233tI999yjv/zlL5KkM844Q4MHD5YktW7duuR7AKmpyn1KZtbRzN4ws5VmttzMbkhGMKCYu2v+/Pnq37+/Tj755JJCqlu3rh544AE+fwSkkXjWlPZIusnd/2ZmTSQtNrP/dfcVtZwN0Hvvvadrr71WH374Yblxo0eP1uGHHx5CKgC1pco1JXdf7+5/i32/XdJKSYfUdjBAko477jiNHTtW3bp1KzO8U6dOuvXWW0NKBaC2VOuQcDPrIuloSYsqGDfKzPLMLG/Tpk2JSYeMV6dOHR155JHasGFDmeHjxo3jAn1AGoq7lMyssaRnJf2Xu2/bd7y7T3T3XHfP5TMiSJQtW7boqaee0mOPPabrrrtOdevW1VlnnaVzzz037GgAakFcpWRmdRUU0nR3f652IwHSJ598oquuukodO3bUN998owsvvFDjx4/XoEGDNH78eA5uANKUFV8ArdIJgv/+JyR97e7/Fc9Mc3NzPS8vLwHxkEncXQsWLND999+vl19+WZJ0+OGHa8mSJSWfPVq/fr3atWsXZkwAB8DMFrt7blXTxXP03YmSLpG01MzyY8NudfeXaxIQKFZYWKinn35a999/v5YuXVpm3GOPPVbmw7AUEpDeqiwld39HEttKkHAbN27UhAkT9PDDD2vjxo3lxl955ZUaMGBACMkAhIUzOiAUBQUFGjJkiBYtKncgpySpffv2+p//+Z8kpwIQNs4SjlB06NBB77//vp588skKx0+YMEHNmjVLcioAYaOUEBoz0xtvvFFu+NChQ3X22WeHkAhA2CglhGbJkiU6+eST1adPH5155pmSpBYtWmj8+PEhJwMQFvYpIal27dql2bNna8KECVq+fLmWLl2qoUOHas+ePWrRooXuu+8+tWnTJuyYAEJCKSEpVq1apUcffVRTp04tuYT55MmT1blz55Jpfv3rX+uyyy4LKyKACKjyw7MHgg/PQpJ2796tF154QRMmTNDrr79eZtxPfvITvfDCC2XOzODunKkBSFOJ/PAsUC1ffvmlHnvsMU2aNKnciVQlqWXLlpo4cWK5AqKQAFBKSIiioiL95S9/0SOPPKKXXnpJe/furXTaRx55RG3btk1iOgCpgqPvkBDurm3btmnt2rX7LaSLLrpIQ4YMSWIyAKmEUkJCHHTQQbrggguUl5enYcOGVTjNIYccogceeCDJyQCkEkoJCbVw4ULNmDGjwnGTJ09WixYtkpwIQCqhlJBQBx98sK688kpJUumLPV511VUlH5AFgMpwoANqbOvWrZo1a5amTp2qpk2b6sUXX1SfPn3Uo0cP9e/fX926ddPYsWPDjgkgBVBKOCBFRUVasGCBpk6dqrlz56qwsFBt2rTRkiVLVK9ePY0cOVKFhYXKzs7WE088oUaNGoUdGUAKoJRQLatWrdITTzyhJ598UgUFBWXGPfnkk2VzgoyqAAANUklEQVROEZSdna1p06bpxBNPTHZMACmKUkKVtm7dqmeeeUZTp07Ve++9V+E0o0eP1hlnnFFuOId/A6gOSgmVevvtt/XII4/oueeeU2FhYaXT5ebm6u67705iMgDpiqPvUKnDDz9chxxyiLKysiqdpnHjxnr66adVr169JCYDkK4oJVSqTZs2uvfeezV58uRKp5kwYYK6d++exFQA0hmlhErt3btXY8aM0YUXXljh+EsuuUQXX3xxklMBSGfsU0KFtmzZomHDhunVV18tGdamTRv94x//kCR1795dDz30UFjxAKQp1pRQoXr16mnNmjUlPw8ePFjLly9XvXr1VLduXc2cOVNNmjQJLyCAtEQpoQx316JFizRmzBj98Y9/VLNmzTRu3DjNnj1bLVu2VK9evfS73/1Offv2DTsqgDTE5jvI3fXRRx9p5syZmjlzptasWaOxY8fqxz/+sdasWaPmzZuXTHvzzTfz2SMAtYZSymAff/yxZs2apZkzZ+rjjz8uGX7eeefppptukqQyhSRJQ4cOTWpGAJmFUsowa9asKSmi/Pz8cuO7d++uKVOmcGlyAKGglDLAunXrNHv2bM2cOVMLFy6sdLoGDRro2WefVbNmzZKYDgD+jVLKACtWrNC0adO0ePHi/U43YcIEHXXUUUlKBQDlcfRdBjjttNP0wQcf6N577610mpEjR+qyyy5LYioAKI9SyhAzZszQXXfdVeG4Pn36aPz48UlOBADlUUppbtu2bSWnA9q+fXu58S1atNCcOXOUnZ0dQjoAKItSSmMLFy5UTk6OnnrqqZJhubm5mjp1asnP06ZNU9euXUNIBwDlUUppqKioSL/97W910kkn6fPPP5ckmZluvvlmvfvuuzr77LMlSb/5zW/0ox/9KMyoAFAGR9+lob1792revHkqKiqSJLVv315PPvmkBg0aJCk4r93FF1+sO++8M8yYAFAOa0pp5vPPP9cDDzygm266SU2aNNE555yjJUuWlBRSsccff3y/F+8DgDBQSinO3ZWfn6877rhDP/jBD3TooYdq6dKlGjx4sBYvXqy5c+eqVatW5X6vbt26IaQFgP1j810KKioq0rvvvqu5c+fq+eefL3OJiWOPPVYTJkyQmemwww4LLyQAHABKKUXs2rVLr732mp5//nm9+OKL2rx5c7lp2rVrp+eee47DuwGkLEop4ubPn6/JkyfrlVde0c6dOyudrn79+nr++efVvn37JKYDgMRin1LEHXfccercubO+++67/U43ceJEHXPMMUlKBQC1g1KKuBYtWuiPf/yj7rnnnkqnufHGG3XppZcmMRUA1A4230XcP//5T1133XWaPn16heNPP/30/Z5oFQBSSZVrSmb2uJltNLNlyQiEf5s/f7569epVppDq169f8n23bt00c+ZMHXQQ7y0ApId4Nt9NlXRWLedAKdu2bdPIkSP1H//xH1q3bp2k4DRBv/zlL0vOY9e4cWO98MILOvjgg8OMCgAJVeVbbHd/y8y61H4USNIbb7yhESNG6IsvvigZ1q1bN02dOlUnnXSSVqxYIUmaPn26evbsGVZMAKgVCTvQwcxGmVmemeVt2rQpUbPNGN9++62uv/56nXrqqWUK6ZprrtGSJUt00kknSZK6dOmiu+++u+SkqgCQTszdq54oWFOa5+694plpbm6u5+Xl1SxZBsnPz9fQoUP16aeflgzr2LGjpkyZUu6cdVJwaiEzS2ZEAKgRM1vs7rlVTcch4RHQvHlzbdiwoeTnyy+/XEuXLq2wkCRRSADSFodthcTdtXTpUr388ssqKCjQ73//e40ZM0aTJk3iGkcAMlaVpWRmT0saKKmVmRVIut3dJ9d2sHS0detWvfbaa5o/f75eeeUVrV27Vm3atNGiRYvUqVMnDRs2TM2aNQs7JgCEJp6j7y5MRpB0VLw2NH/+fM2fP1/vvvuu9uzZUzK+QYMG+vOf/6zOnTtLEoUEIOOx+S7Btm3bVrI2NH/+fK1du7bC6cxM06dPV79+/ZKcEACii1JKoM2bN2vIkCF68803q5x27NixOu+885KQCgBSB0ffJVCrVq30+uuv64knnlDLli0rne7KK6/UL37xiyQmA4DUQCkl2HfffadPPvlEW7durXD8mWeeqQceeIDDugGgAmy+S6A333xTo0aN0qpVqyoc37t3bz3zzDOcQBUAKsGaUgJs3bpVP//5zzVw4MAyhXTuueeqTp3gIW7btq3mzZunpk2bhhUTACKPUqqh559/Xj169NDEiRNLhrVt21bPPvus5s6dqzZt2qhhw4aaN2+eOnXqFGJSAIg+SukArV+/XkOGDNF5552n9evXlwwfOXKkVq5cqZ/+9KeSpE6dOmnGjBnq27dvWFEBIGWwc6Oa3F2TJ0/WL3/5yzIHM3Tv3l0TJ07UKaecUmb6hx9+WH369El2TABISawpVcOnn36qU089VSNHjiwppKysLN1yyy366KOPyhWSJAoJAKqBNaU47N69W/fdd5/uuOMOFRYWlgzv27evJk2apJycnBDTAUD6oJSqsGPHDvXv31/5+fklwxo0aKC77rpLN9xwA4d3A0ACsfmuCo0bN1aPHj1Kfj7ttNO0bNky3XTTTRQSACQYpRSHcePG6bDDDtOUKVP06quv6tBDDw07EgCkJd7qx+F73/ueVq5cqaysrLCjAEBaY00pThQSANQ+SgkAEBmUEgAgMiglAEBkUEoAgMiglAAAkUEpAQAig1ICAEQGpQQAiAxKCQAQGZQSACAyKCUAQGRQSgCAyKCUAACRQSkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFBKQEAIoNSAgBEBqUEAIgMSgkAEBmUEgAgMiglAEBkUEoAgMiIq5TM7Cwz+8TM/m5mt9R2KABAZqqylMwsS9JDkn4o6UhJF5rZkbUdDACQeeJZUzpG0t/d/TN3/07STEnn1G4sAEAmiqeUDpH0VamfC2LDyjCzUWaWZ2Z5mzZtSlQ+AEAGiaeUrIJhXm6A+0R3z3X33NatW9c8GQAg48RTSgWSOpb6uYOkdbUTBwCQyeIppQ8kHWZmXc2snqQLJL1Yu7EAAJnooKomcPc9ZnatpL9IypL0uLsvr/VkAICMU2UpSZK7vyzp5VrOAgDIcJzRAQAQGZQSACAyKCUAQGRQSgCAyKCUAACRQSkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZFBKAIDIoJQAAJFBKQEAIoNSAgBEBqUEAIgMSgkAEBmUEgAgMiglAEBkUEoAgMgwd0/8TM02Sfoi4TMur5WkzUm4n2RjuVJHOi6TlJ7LlY7LJKXOcnV299ZVTVQrpZQsZpbn7rlh50g0lit1pOMySem5XOm4TFL6LReb7wAAkUEpAQAiI9VLaWLYAWoJy5U60nGZpPRcrnRcJinNliul9ykBANJLqq8pAQDSCKUEAIiMlC0lMzvLzD4xs7+b2S1h50kEM3vczDaa2bKwsySKmXU0szfMbKWZLTezG8LOlAhmlm1m/2dmS2LLdWfYmRLFzLLM7EMzmxd2lkQxszVmttTM8s0sL+w8iWJmzc1sjpl9HPsfOz7sTDWVkvuUzCxL0ipJp0sqkPSBpAvdfUWowWrIzAZI2iHpSXfvFXaeRDCzdpLaufvfzKyJpMWSzk2D58okNXL3HWZWV9I7km5w94UhR6sxM/uFpFxJTd39x2HnSQQzWyMp191T4UOmcTOzJyS97e6TzKyepIbu/k3YuWoiVdeUjpH0d3f/zN2/kzRT0jkhZ6oxd39L0tdh50gkd1/v7n+Lfb9d0kpJh4SbquY8sCP2Y93YLfXe4e3DzDpI+pGkSWFnwf6ZWVNJAyRNliR3/y7VC0lK3VI6RNJXpX4uUBq80KU7M+si6WhJi8JNkhixzVz5kjZK+l93T4flGifpV5L2hh0kwVzSq2a22MxGhR0mQQ6VtEnSlNjm1klm1ijsUDWVqqVkFQxL+Xep6czMGkt6VtJ/ufu2sPMkgrsXuXuOpA6SjjGzlN7kamY/lrTR3ReHnaUWnOjufST9UNI1sU3lqe4gSX0kTXD3oyXtlJTy+9dTtZQKJHUs9XMHSetCyoIqxPa5PCtpurs/F3aeRIttMvmrpLNCjlJTJ0o6O7b/ZaakU83sqXAjJYa7r4t93ShproJdAKmuQFJBqTX0OQpKKqWlail9IOkwM+sa27l3gaQXQ86ECsQOCJgsaaW73xd2nkQxs9Zm1jz2fQNJp0n6ONxUNePuv3b3Du7eRcH/1OvufnHIsWrMzBrFDrJRbPPWGZJS/ghXd98g6SszOyI2aJCklD6ASApW/1KOu+8xs2sl/UVSlqTH3X15yLFqzMyeljRQUiszK5B0u7tPDjdVjZ0o6RJJS2P7XyTpVnd/OcRMidBO0hOxI0HrSHrG3dPmEOo000bS3OD9kQ6SNMPdXwk3UsJcJ2l67M35Z5JGhJynxlLykHAAQHpK1c13AIA0RCkBACKDUgIARAalBACIDEoJABAZlBIAIDIoJQBAZPw/mxn1fdGPbloAAAAASUVORK5CYII=\n",
"text/plain": [
"