Estimate_and_Inference_without_pyro.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Untitled",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/hellocybernetics/Tensorflow-Probability-Tutorials/blob/master/Estimate_and_Inference_without_pyro.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "metadata": {
        "id": "hA-qlJZ9Nw4v",
        "colab_type": "code",
        "outputId": "0bfa0739-8e5c-471c-8cb1-65b7811f6dee",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 197
        }
      },
      "cell_type": "code",
      "source": [
        "!pip install pyro-ppl"
      ],
      "execution_count": 179,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: pyro-ppl in /usr/local/lib/python3.6/dist-packages (0.3.0)\n",
            "Requirement already satisfied: torch>=1.0.0 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (1.0.0)\n",
            "Requirement already satisfied: graphviz>=0.8 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (0.10.1)\n",
            "Requirement already satisfied: six>=1.10.0 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (1.11.0)\n",
            "Requirement already satisfied: opt-einsum>=2.3.0 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (2.3.2)\n",
            "Requirement already satisfied: numpy>=1.7 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (1.14.6)\n",
            "Requirement already satisfied: networkx>=2.2 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (2.2)\n",
            "Requirement already satisfied: contextlib2 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (0.5.5)\n",
            "Requirement already satisfied: tqdm>=4.28 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (4.28.1)\n",
            "Requirement already satisfied: decorator>=4.3.0 in /usr/local/lib/python3.6/dist-packages (from networkx>=2.2->pyro-ppl) (4.3.2)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "459cvjGfOGyN",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Pyroの基本的な役割\n",
        "PyroはPyTorchを計算のバックエンドに構えた確率プログラミング言語です。確率プログラミング言語の主な役割は、様々な確率分布からのサンプリングや、同時分布・条件付き分布・周辺分布の取扱を容易にすることです。**まずはTorchを使ってみて、その確率プログラミング言語の必要性を体感してみましょう。**"
      ]
    },
    {
      "metadata": {
        "id": "77PKI6F3N6Ec",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import torch\n",
        "import torch.distributions as torchdist\n",
        "\n",
        "import pyro\n",
        "import pyro.distributions as dist\n",
        "import pyro.poutine as poutine"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "LSELqDR0PTtu",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### 分布の記述\n",
        "まず、肩慣らしに標準正規分布 $\\rm {Normal} (0, 1)$ を書いてみましょう。`torch.distributions`  モジュールを使うことで下記のように記述することができます。"
      ]
    },
    {
      "metadata": {
        "id": "bPIXw1EnPDU2",
        "colab_type": "code",
        "outputId": "64b2974b-3c4d-4019-ae09-3b81f3aed621",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "cell_type": "code",
      "source": [
        "normal_dist = torchdist.Normal(loc=torch.tensor(0.), scale=torch.tensor(1.))\n",
        "normal_dist"
      ],
      "execution_count": 181,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Normal(loc: 0.0, scale: 1.0)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 181
        }
      ]
    },
    {
      "metadata": {
        "id": "k93zrklnPuis",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### サンプリング\n",
        "続いて、記述した分布から\n",
        "$$\n",
        "x \\sim {\\rm Normal}(0, 1)\n",
        "$$\n",
        "とサンプリングを得るには下記のようにします。"
      ]
    },
    {
      "metadata": {
        "id": "9Apgnf6sQFe0",
        "colab_type": "code",
        "outputId": "9fe6f5ac-a5a6-4e6c-82d0-730da0b3a2fd",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "cell_type": "code",
      "source": [
        "x = normal_dist.sample()\n",
        "x"
      ],
      "execution_count": 182,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor(0.0850)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 182
        }
      ]
    },
    {
      "metadata": {
        "id": "mq-6TcUIQPqS",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "更に同じ分布から独立に複数のサンプルを得たい場合には下記のように、`sample()` メソッドに引数を渡すことで実施します。このとき引数はリストあるいはタプルで渡すようにします。"
      ]
    },
    {
      "metadata": {
        "id": "kc0-utrYQnVK",
        "colab_type": "code",
        "outputId": "98a8a80d-220e-43a6-ad3f-161eb0781769",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "cell_type": "code",
      "source": [
        "X1 = normal_dist.sample([3,])\n",
        "X1"
      ],
      "execution_count": 183,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor([ 2.1259,  1.4169, -2.0609])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 183
        }
      ]
    },
    {
      "metadata": {
        "id": "Puq8R9cqRihc",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "上記のサンプルは各成分が各々独立に正規分布から生成されています。数式で表現するならば下記のようになります。\n",
        "\n",
        "$$\n",
        "X_1[i] \\sim \\rm{Normal}(0, 1)\n",
        "$$\n",
        "\n",
        "同じ理屈でもっと多次元の配列を作ることもできます。"
      ]
    },
    {
      "metadata": {
        "id": "VzE_FNBCQqKa",
        "colab_type": "code",
        "outputId": "138c605a-8b57-43a4-c749-20afaef1a5c7",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 197
        }
      },
      "cell_type": "code",
      "source": [
        "X2 = normal_dist.sample([10, 2])\n",
        "X2"
      ],
      "execution_count": 184,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor([[-0.8648,  0.5027],\n",
              "        [-0.8110, -0.5906],\n",
              "        [-0.5288,  0.2666],\n",
              "        [-1.3417,  0.3022],\n",
              "        [-1.7919, -3.0321],\n",
              "        [ 1.1864, -0.3182],\n",
              "        [-0.0318, -0.5959],\n",
              "        [ 1.1256, -0.6396],\n",
              "        [-1.1621,  1.4708],\n",
              "        [ 0.7577, -0.8019]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 184
        }
      ]
    },
    {
      "metadata": {
        "id": "yfcsEdNdRvic",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "上記の配列も同様に各成分が独立に正規分布から生成されています。\n",
        "\n",
        "$$\n",
        "X_2[i, j] = \\rm{Normal}(0, 1)\n",
        "$$\n",
        "\n",
        "配列として取り出すと、もはや元々これが何の分布の話だったのか、情報を持ち合わせていないので注意が必要です。例えば、下記のような2次元正規分布を考えましょう。2次元正規分布からの1つのサンプルは当然2つの要素を持っています。"
      ]
    },
    {
      "metadata": {
        "id": "K5GrSGa7RBAK",
        "colab_type": "code",
        "outputId": "5850447a-39ff-4047-f019-e4dfc476c091",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "cell_type": "code",
      "source": [
        "multi_normal_dist = torchdist.MultivariateNormal(\n",
        "    loc=torch.tensor([1., -2.]),\n",
        "    covariance_matrix=torch.tensor([[1., -1.],\n",
        "                                    [-1., 2.]])\n",
        ")\n",
        "\n",
        "multi_normal_dist.sample()"
      ],
      "execution_count": 189,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor([ 0.5684, -3.4546])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 189
        }
      ]
    },
    {
      "metadata": {
        "id": "8gKtZFYpVB0C",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "この2次元正規分布から10個のサンプルを得ると、`(10, 2)` の `tensor`が得られますが、この `tensor` の各成分は独立ではありません。各行は独立ですが、各列は2次元の正規分布からの一つのサンプルであり、相関を持っていることに注意する必要があります。\n",
        "\n"
      ]
    },
    {
      "metadata": {
        "id": "M1BGa6oXUILa",
        "colab_type": "code",
        "outputId": "a2af3be1-96d4-4d58-e9fc-c989d60a0a65",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 197
        }
      },
      "cell_type": "code",
      "source": [
        "multi_normal_dist.sample([10])"
      ],
      "execution_count": 190,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor([[ 0.5396, -2.3828],\n",
              "        [ 2.3775, -4.0210],\n",
              "        [ 0.3301, -2.0491],\n",
              "        [ 0.8666, -1.0754],\n",
              "        [ 0.7196, -0.7638],\n",
              "        [ 1.1182,  0.0462],\n",
              "        [ 1.2652, -0.1359],\n",
              "        [ 1.9287, -3.0320],\n",
              "        [ 1.2091, -2.2953],\n",
              "        [ 0.2908, -1.6284]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 190
        }
      ]
    },
    {
      "metadata": {
        "id": "sHtUrjA7VszA",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "長いコードを書いているうちに、ある`tensor` がどういうサンプルのされ方をしたものであったのかを見失うことも起こりうるので注意が必要です。"
      ]
    },
    {
      "metadata": {
        "id": "ED47SK9QVuBv",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### 尤度の計算\n",
        "機械学習ではしばしば確率分布 $p(x)$ の尤度の計算が必要になります。あるサンプル $x^*$ の $p(x)$ における尤度は単純に $p(x^*)$ の値になります。では $\\{x_1, x_2, x_3,...,\\}$ の尤度はどのように計算するのかというと、\n",
        "\n",
        "$$\n",
        "p(x_1)p(x_2)p(x_3)...\n",
        "$$\n",
        "\n",
        "と各確率値の積によって表されます。確率値は $[0, 1]$ の値でありサンプルの数が多くなると、$1$ 未満の数の積が連なり非常に小さな値となってしまいます。そこで、コンピュータでは通常、対数尤度というものを計算することにします。\n",
        "\n",
        "$$\n",
        "{\\rm log}\\{p(x_1)p(x_2)p(x_3)...\\} = {\\rm log}p(x_1) + {\\rm log}p(x_2) + {\\rm log}p(x_3) + ...\n",
        "$$\n",
        "\n",
        "と負の値の足し算になり数値的に安定します。また、$\\rm log$ は単調増加関数であるので、微分にとっても符号には影響が無く（絶対値には影響するが）、積が和となっているため都合が良いです。サンプル $x^*$ の対数尤度は簡単に下記のコードで求まります。"
      ]
    },
    {
      "metadata": {
        "id": "1LnrmXtHWQs8",
        "colab_type": "code",
        "outputId": "fb49e9b4-02cf-41c2-cf3f-82d5a15e619c",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 89
        }
      },
      "cell_type": "code",
      "source": [
        "sample = multi_normal_dist.sample()\n",
        "print(\"*** a sample from 2dim normal: \\n\", sample)\n",
        "log_likelihood = multi_normal_dist.log_prob(sample)\n",
        "print(\"*** log likelihood of the sample \\n\", log_likelihood)"
      ],
      "execution_count": 191,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "*** a sample from 2dim normal: \n",
            " tensor([0.2084, 0.0300])\n",
            "*** log likelihood of the sample \n",
            " tensor(-2.9180)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "qZb5n79oWebI",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "また複数のサンプルに関しては足し算する前の各々の対数尤度が返ってきます。"
      ]
    },
    {
      "metadata": {
        "id": "KRHth-wsYfF-",
        "colab_type": "code",
        "outputId": "2b61f3f4-b424-45fd-86db-d4c568ca2697",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 89
        }
      },
      "cell_type": "code",
      "source": [
        "samples = multi_normal_dist.sample([5,])\n",
        "print(\"*** 5 samples from 2dim normal: \\n\", sample)\n",
        "log_likelihoods = multi_normal_dist.log_prob(samples)\n",
        "print(\"*** log likelihoods of the samples: \\n\", log_likelihoods)"
      ],
      "execution_count": 192,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "*** 5 samples from 2dim normal: \n",
            " tensor([0.2084, 0.0300])\n",
            "*** log likelihoods of the samples: \n",
            " tensor([-2.5051, -4.1052, -2.2041, -2.2630, -2.2343])\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "wy4ONNQrYkgi",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### MAP推定\n",
        "ここでMAP推定を実施してみます。上記の関数と自動微分機能を駆使すれば難しくありません。\n",
        "\n",
        "#### 訓練データ\n",
        "人工訓練データを作りますが、パラメータは知らない体で推定をします。"
      ]
    },
    {
      "metadata": {
        "id": "M1mu_XCBZ-8e",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def toy_data():\n",
        "    true_mu = torch.tensor([2., 3.])\n",
        "    true_cov = torch.tensor([[2., -3.],[-3, 5]])\n",
        "    \n",
        "    X = torchdist.MultivariateNormal(loc=true_mu, \n",
        "                                     covariance_matrix=true_cov).sample([100,])\n",
        "    return X\n",
        "\n",
        "X_train = toy_data()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "RePK_G1Kac-j",
        "colab_type": "code",
        "outputId": "1d45d7c7-6ea5-48b2-ad51-35342320deed",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 369
        }
      },
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.plot(X_train.numpy()[:,0], X_train.numpy()[:,1], \"o\")"
      ],
      "execution_count": 194,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7fdf76f8f0b8>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 194
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFOCAYAAADHOhe+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH85JREFUeJzt3W1wVeXd7/FfSAwaEiKETSRGwbG9\nK0nHEVtN4wO1LbdOjZMZ7VCwt9qZWoszhWnHYmuxox0dnQH7Qout0IrVGVplwljlHDpitWT0dGK8\nRY8doQ9gDw8BjDsQSUgIN0n2edEmkGSvnbXXXuu61sP384qsvbP2tcD2t//XY1Emk8kIAAAYM8V2\nAwAASBrCFwAAwwhfAAAMI3wBADCM8AUAwDDCFwAAw0q8/mJfX59+9KMf6dixYzp16pS++93v6tpr\nr/WzbQAAxJLn8P3973+viy66SD/4wQ/U2dmpb37zm3rllVf8bBsAALHkudt5xowZ+uSTTyRJPT09\nmjFjhm+NAgAgzooK2eHqzjvv1P79+9XT06P169frsssu87NtAADEkufK9+WXX1ZNTY3++Mc/6rnn\nntNDDz2U8/2Dg0NePwoAgFjxPOb77rvv6pprrpEkXXLJJfr44481NDSk4uLirO/v7u73+lG+SaUq\nlE732m6GUTxzMvDMyZDEZ5ai+9ypVIXja54r37lz5+r999+XJB08eFDTpk1zDF4AAHCa58p3yZIl\nWrVqlW677TYNDg7qpz/9qY/NAgAgvjyH77Rp0/TEE0/42RYAABKBHa4AADCM8AUAwDDCFwAAwwhf\nAAAMI3wBADDM82znOGjf1amtbXt1qKtfNbPK1NQ4Tw111babBQCIucRWvu27OrV+y051pPs0nMmo\nI92n9Vt2auUv/qz2XZ22mwcAiLHEhu/Wtr1Zrx/tPan1W3YSwACAwCQ2fA915d5remvbPkMtAQAk\nTWLHfGtmlakj3ef4+uEjzq+FAePVABBdia18mxrn5Xx9TtU0Mw3xwGm8mq5yAIiGxIZvQ121ljXX\na+b0qVlfb2qca7hF7jmNV9NVDgDRkNhuZ+lfAdxQV/3vLtx9OnykT3OqpqmpcW6ou3CdxqvD3lUO\nAPiXRIfviJEQjgqn8eowd5UDAE5LbLdzlDmNV4e5qxwAcBqVbwSNVOlR6ioHAJxG+EZU1LrKAQCn\n0e0MAIBhkat82VwCABB1kQrfkc0lRoxsLiGJAAYAREakup3ZXAIAEAeRCl82lwAAxEGkwrdmVlnW\n62wuAQCIkkiFL5tLAADiIFITrthcAgAQB5EKX4nNJQAA0RepbmcAAOIgcpVvodikAwBgW6LC9433\nOtikAwBgXaK6nVte3531Opt0AABMSlTlu7+zN+v1OG7SMdq9fqRfNVV0rwNAmCSq8r2wuiLr9bht\n0jGyB3ZHuk/Dw5nR7vX2XZ22mwYAUMLCd/FXPp31etw26WAPbAAIt0R1Oy9cUKuenoHYb9LBHtgA\nEG6JCl8p/006org0qWZWmTrSE4M2bt3rABBViep2zteYsdNMdMZO2QMbAMItcZVvPnKNnYa5+mUP\nbAAIN8I3hyiPnY50r6dSFUqnsy+xAgDYQbdzDpwfDAAIAuGbA2OnAIAg0O2cA2OnAIAgEL6T4Pxg\nAIDfCF8EJoprpAHABMIXruQbpCNrpEdwfCMAnEb4xkDQFaaXII3qGmkAMKGg2c5btmxRc3Ozbrnl\nFrW2tvrUJOTDxC5cXg5qiPIaaQAImufw7e7u1i9+8Qv97ne/07p16/T666/72S64ZOIEIy9Byhpp\nAHDmOXzb2trU2Nio8vJyzZ49Ww8//LCf7YJLJipML0HKGmkAcOY5fDs6OjQwMKC7775b3/jGN9TW\n1uZnu+CSiQrTS5A21FVrWXO9alPlKp5SpNpUuZY11zPeCwCSijKZTMbLL/7qV7/Su+++qyeffFKH\nDh3SHXfcoe3bt6uoqCjr+wcHh1RSUlxQYzHRG+916LGNOyZcv/e2z2nhglpfP6fl9d060NmrC6or\ntPgrn/b1/gCQJJ5nO1dVVWnBggUqKSnRhRdeqGnTpuno0aOqqqrK+v7u7uzdoyZF5ZCBfGYvz6+t\n1LLm+gm7cM2vrVQ63evbM8+vrdQD3/z8mGuF3DfIGdpunzlO65Cj8t+2n3jm5Ijqc6dSFY6veQ7f\na665Rvfdd5/uuusuHTt2TP39/ZoxY4bX2+HfvCzridouXGFYAxyGNgBILs/hW11drRtuuEFf//rX\nJUk/+clPNGUK5zQUysv62HwqONvVXvuuTj2zdVfW10yuAWYdMgCbCtpkY+nSpVq6dKlfbYHyn72c\nTwVnu9ob//njmVwDzDpkADZRqoZMvrOXnSq49Vt26oEN7XrjvY5J3+vnmuBcnD5/RGV5qZF2SKxD\nBmAX4Rsy+S7rcargpH9Vto9t3DG625Xtai9XWyVJnubde8M6ZAA2Eb4hk+/6WKcK7kwjla3tam+y\nth7r+x8j7ZBYhwzALg5WCKF8Zi83Nc7LOY4qna5snd5rqtqbrK2mu3yjNkscQHwQvhE3Eh5b2/ap\nI30863tGQu3M9565JthUAI18TkvrHh3tOTnhdbp8ASQF4RsDIxWc02ziM0PNdrV3ZlttfQkAANsI\n3xjJVtneesNnNL+20nLLJrL9JQAAbCJ8Y2Z8qEV1WzYAiDNmOwMAYBiVb4TY3hrSq6i2GwCCQvhG\nhO2tIb2KarsBIEh0O0eE7a0hvYpquwEgSIRvRNjeGtKrqLYbAIJEt3NE1MwqU0d6YmAFvStUoeO1\nttoNAGFG5RsRNg4CGBmv7Uj3aTiTGR2vHTmowQ0OMACAiah8I8LG1pB+HDhve0tLAAgjwjdCTO8K\n5dd4LbtZAcBYhC8cOY3XVk4r1QMb2lm3CwAeMeYLR07jtUd7TxY0DgwASUf4wlG2A+dnVkzN+l7W\n7QKAe3Q7I6fx47XfXr096/tsrNtl20oAUUXli7zUzCrLet30ul0/lkEBgC1UvshLU+O8MXs1n77u\nbd2u1+rVaRlUS+sebfvvA9r/US/VMIDQInyRFz/X7RZy6ILTMqijPSd1tOdk3vcDAJMIX+TNr3W7\nhWzi4bQMKpuW7XsIXwChQvhiUkFNbCpkEw+n7u9sjvaeVPuuTldtZhIXABMIX+SUq2tYUkFBdW55\nqY72npxw3c3krWzd3/0Dp7Leb+R9k7WNs4cBmEL4IifHiU3b94wJunyDqn1Xp2NQup28Nb77e3x4\nnslNNe3HXtYA4AZLjZCT48SmHBWmG05BN7Niquega6irdtwExE01zdnDAEyh8kVO+UxsktwHlVPQ\nHev7H0nex14Xf+lTnpdCcfYwAFOofJGT0/7OhVSYUu7NOgrZQKOhrlr33va5MVtiLmuudxXcnD0M\nwBQqX+TktK5XUkGbbeTarKPQsdeFC2o1v7bSVTvOxNnDAEwhfDGpXOt6vQZVrqD79f/alfV3TIy9\ncvYwABMIX3hWaFA5/b6bsVfW4wKIMsIXoTPZ/tFhXY/LFwIAbhG+CJ3Jxl7DuB43rF8IAIQT4YtQ\nytWlHcb1uGH8QgAgvFhqhMgJy5nCZwrjFwIA4UXli8jx+0xhPzhNEqssL9UDG9oZBwYwBuGLyPGy\nHjfoyVBOXwiO9pzUUXG+MICxCF9EUj7LnExMhsrnlKVntv7V188GED2EL2LP1GSo8V8Ivr16e9b3\nnRoapgIGEo4JV4g9W5OhnCaGjXB7AhSA+CF8EXu2Zkc7HdQwgpnQQHIRvog9W6cVNdRVa1lzvc4q\nzv4/M44qBJKroPAdGBjQokWL9OKLL/rVHsB3IyHo5ZhBPz77W03zs77GUYVAchU04eqpp55SZWX+\nR7cBptk8rYijCgGM5zl8P/zwQ+3Zs0fXXXedj80B4omjCgGcyXO38+rVq3Xffff52RYAABLBU+X7\n0ksv6bLLLtMFF1zg+ndmzChTSUmxl4/zVSpVYbsJxvHMdrzxXodaXt+t/Z29urC6Qou/8mktXFAb\n2OeF4ZlN45mTI27P7Sl8W1tbdeDAAbW2tuqjjz5SaWmpzjvvPF111VWOv9PdnX2tpUmpVIXS6V7b\nzTCKZ7Zj/K5aew/36LGNO9TTMxBI93MYntk0njk5ovrcub4weArfxx9/fPTPa9eu1fnnn58zeIGk\n4YhBALmwvSQQALe7ark98CHogyFMitOzAF4VHL4rVqzwox1ArDgdMXjmxhpuD3wwcTCEKXF6FqAQ\n7HAFBMDNrlq5uqbH/uzufVEQp2cBCkG3M2LPVjfnzIqpo0cKzpw+VYuv+9SYz3XbNW3rYIggxOlZ\ngEJQ+SLWRro5O9J9Gs5kRrs523d1Bv6ZZ57le7Rn4rm+bg98sHUwRBDi9CxAIQhfxJqNbk63n+n2\nwAdbB0MEIU7PAhSCbmfEmo1uTref6XbP5zjtDe32WZgRjbgjfBFrbmYd2/xMt3s+x2lv6MmehRnR\nSAK6nRFrNro56VotDDOikQRUvog1G122ceomtoEZ0UgCwhexZ6PLNk7dxKbZGCoATKPbGUCo0G2P\nJKDyBQxhBq87dNsjCQhfwICwzOCNyhcAuu0Rd4QvYECQRwy27+rUtv9+R/s/6p30ZKQwfAEAQPgC\nRvh9xOCZ73cbqJwxDIQHE64AA9zsaexlH+p81sSyhAcIDypfIAu/x0abGueNqVBPX3d3xKDTZ+cT\nqIUs4YnKWDEQFVS+wDhBnITUUFetZc31qk2Vq3hKkWpT5VrWXO/piMEz5XNKkNclPDZOhgLijsoX\nGCeosdHJZvB6qUzdVNTS6cq1SFJJ8RQNDQ+rZla5qyU8jBUD/qPyBcaxNTbqpTIdqajnzZnuWFGf\nWblmJJ0aGtZwRq7XzjJWDPiPyhcYZ7IKNKjxz3w2lxjfhltvuETzayuz3rfQypXtHgH/Eb7AOLm6\ncoNeK+tmc4lsbXhs444JFe+IQitXt13bANwjfIFxclWgD2xoz/o7+Y5/FlI951vJFlq5st0j4D/C\nF8jCqQJ1U0VOFqyFVs/5VrJ+VK5s9wj4iwlXQB4mW9rjZllOoYfF57O8SHK3zAmAWVS+QB4mqyLd\ndAmbHIMdX4V/+6Y6V6HLphpAsAhfIA+TjX+6CdYgxmBvveEzE2Y7e+3e5gAGIHiEL2Lnjfc69Py2\nvwVWteUa/3QTrEGMwaZSFUqne8e8x+3ErPFVbv/AoKvfA+Ad4YtYsV21uQlWU7OH3U4OG//35YRN\nNQD/EL6IFdtbIboNVhOzh91U4U5/X9lMKSpS+65Oql/AB4QvYiUMWyGGZVmOmyrc6e8rm1NDw4z9\nAj4hfBErbIV4mlMVLkkPbGjXoa5+FU+Rhofyuy9jv0DhCF9EymRLYOK4FWIhy37GV+Hjx3jzDV7J\nXC8Cy50QZ4QvIsPNZKqGumpNn362nt/291hshej3BDKnMd6ziqdocHhYmczk9zDRi2B74hwQNMIX\nkeF2MtXCBbWOJ/xEjd8TyJzGeIczGW340Zf/XW3+q5u6clqpjvaenPBeE70ItifOAUEjfBEZYZhM\nZZrfzzzZmHi2bmobByok8d8ayUL4IjKSOJnK72fOd0w86JnbTuO6Sfy3RrJwsAIio6lxnsP16E6m\nmozfzxymQxZyHUKRxH9rJAuVLyIjiefKBvHMflezbmcl57ON5UN3Xjn656T8WyNZCF9ESlg2sDAp\nzM/sdlayl20sw/zcQKEIXyAhglg363ZWcj7bWDKuiyQgfIEECGrdrNtZyflsY8m4LpKA8AViJluF\nG9S6Wbezkp3eN3P6VJVNPYtxXSQO4QvEiFOFW+Tw/kLXzbpduuT0vsXXfYqwRSIRvkCMOFW4JcVT\ndGpoeML1QsdX8zlC0c37gKQgfIEYcRpbHRqeGLySP+Orbmclm5y9zKEMCLuCwnfNmjXasWOHBgcH\ntWzZMl1//fV+tQuAB05jqzWzytXUODcRlSeHMiAKPIfvW2+9pd27d2vTpk3q7u7WzTffTPgCluUa\ng03KulkOZUAUeA7fK664Qpdeeqkkafr06Tpx4oSGhoZUXFzsW+MA5IexVQ5lQDR4Dt/i4mKVlZVJ\nkjZv3qyFCxcSvEAImKhwgx5T9bplZVPjPA5lQCQUZTJujs929tprr2n9+vV65plnVFFR4fi+wcEh\nlZQQzkDUvfFehx7buGPC9Xtv+5wWLqgN/P5vvNehltd3a99HPcr2/143XXOR/vf/+X+u2zdyv/2d\nvbqwukKLv/JpX54DyKWgCVdvvvmm1q1bp6effjpn8EpSd7f7HW6CkkpVKJ3utd0Mo3jmZDD5zM9v\n+5vD9b9rfm1loPfv6RnIOqZ9pv/797SWNddP6HqfX1s54e9o/OSsvYd79NjGHerpGQhlV30S/9uW\novvcqZRzLnoO397eXq1Zs0bPPvuszj33XK+3ARAxQY+p5rq/mz2iDx/pc931zuQs2OI5fP/whz+o\nu7tb3//+90evrV69WjU1Nb40DEA4BT2mmuv+h7omD/h82sHkLNjiOXyXLFmiJUuW+NkWABHgdktJ\nN7JNmMp1/61te3MeR5hvO5icBVum2G4AgGhpqKvWsuZ61abKVTylSLWpci1rrs+7m3ZkvLUj3afh\nTGbMZhhO929qnJf1XlOK5KkdTvfjZCUEje0lAeTNj+VMucZbH7rzyqz393sdM+uiYQvhC8AKp/HW\ng13Hc/7eSPD7NQM2KTt/IVwIXwBWOI23ZjLSb//4D/3Xf/6H53tzsALCjjFfAFY4jbdK0us7OtS+\nq9PTfZ3Gkr3eDwgC4QvAioa6ahUVOb++tW2fp/vmGksGwoJuZwDWnD9rmuPSIa9rbU2s3aVbG4Ui\nfIEE8Dss/Lqf05peyfta26DX7nJeMPxA+AIx53dY+Hm/hrpq7Tl4TK/v6Jjw2mRrbd94r0PPb/vb\nhC8Afm4Ckk1QW1JSTScL4QvEnN9h4ff9/us//0OfOr8yr7W2br4ABLV2N4hubarp5CF8gZjzOyyC\nCJ9819pO9gUgyLW7QXRrc8BD8jDbGYi5mlllWa8XMqbq5/28sHkgQhBbUnLAQ/IQvkDM+R0WYdgP\n2eYXAL/2tj5TGL7QwCy6nYGYi+N+yEFPqpqM393atp8H5hG+QAL4HRZ+3s/LLN+GumpNn362nt/2\n91gciBCGLzQwi/AFYE0hs3wXLqjV/NrKQNtnEgc8JAtjvgCsYStIJBXhC8AaZvkiqQhfANYwyxdJ\nRfgCsCYMy5YAG5hwBcAaZvnmNn4m+K03XBKrSWZJRvgCsIpZvtllmwn+2MYdBW/ogXAgfAHAoyBP\nImK/53gjfAHAg6BPImImeLwx4QoAPAh6jTIzweON8AUAD4KuTJkJHm90OwMInSDHUv0SxLm+Z8o2\nE/zWGz7DbOeYIHwBhEq+Y6m2gtrESUTjZ4KnUhVKp3t9uz/sIXwBhEo+s3yDnvSUC2uUUQjCF0Co\n5DOWans5DmuU4RUTrgCESj6zfFmOg6gifAGESj6zfFmOg6gifAGESkNdtZY116s2Va7iKUWqTZU7\nbqnIchxEFWO+AELH7Vgqk54QVYQvgEhj0hOiiG5nAAAMo/IFAMOisIMXgkX4AoBBNjcGQXgQvgDg\nIIgK1fbGIAgHwhcAsgiqQmVjEEhMuAKArII6r5eNQSBR+QJIiHy7kIOqUE2choTwI3wBxJ6XLuSg\nzutlYxBIhC+ABPAyySnICpWNQUD4Aog9L13IVKgIkufwffTRR/X++++rqKhIq1at0qWXXupnuwDA\nN167kKlQERRPs53ffvtt7du3T5s2bdIjjzyiRx55xO92AYBvOP0IYeOp8m1ra9OiRYskSRdffLGO\nHTum48ePq7y83NfGAYAf6EJG2HgK366uLtXX14/+PHPmTKXTacIXQGjRhYww8WXCVSaTmfQ9M2aU\nqaSk2I+PK0gqVWG7CcbxzMnAMydDEp9Zit9zewrf2bNnq6ura/Tnjz/+WKlUKufvdHdnn21oUipV\noXS613YzjOKZk4FnToYkPrMU3efO9YXB04Srq6++Wtu2bZMk7dy5U7Nnz6bLGQAAlzxVvpdffrnq\n6+u1dOlSFRUV6cEHH/S7XQAAxJbnMd+VK1f62Q4AABKDU40AADCM8AUAwDDCFwAAwwhfAAAMI3wB\nADCM8AUAwDDCFwAAwwhfAAAMI3wBADCM8AUAwDDCFwAAwwhfAAAMI3wBADDM86lGAAB/tO/q1Na2\nvTrU1a+aWWVqapynhrpq281CgAhfALCofVen1m/ZOfpzR7pv9GcCOL7odgYAi7a27XW4vs9oO2AW\n4QsAFh3q6s96/fCRPsMtgUmELwBYVDOrLOv1OVXTDLcEJjHmCwAWNTXOGzPme/r6XKPtYNKXWYQv\nAFg0EnBb2/bp8JE+zamapqbGuUaDj0lf5tHtDACWNdRVq6lxruZUlelQV5+2tu1V+65OY5/PpC/z\nqHwBwDLblSeTvsyj8gUAy2xXnkz6Mo/wBQDLbFeeTY3zHK6bnfSVJHQ7A4BlNbPK1JGeGLSmKs8w\nTPpKGsIXACwLw3KjhrpqwtYgwhcALKPyTB7CFwBCIKyVJ5tvBIPwBQBkZWIJVFLDndnOAICsgl4C\nNRLuHek+DWcyo+FucoMRWwhfAEBWQS+Bsr2+2SbCFwCQVdCbb9he32wTY74AkABexlaDXgJle32z\nTVS+ABBzXsdWG+qqtay5XrWpchVPKVJtqlzLmut9mxCV5J21qHwBIOZyja1OFqRBLoFK8vpmwhcA\nYi7MY6thXd8cNLqdASDmOLUofKh8ASDC3EykCsPe0RiL8AWAiHK7A1WSx1bDivAFgIjKZyJVUsdW\nw4oxXwCIqDBPpEJuhC8ARBQTqaKL8AWAiEryJhVRx5gvAEQUE6mii/AFgAhjIlU0eQrfwcFB3X//\n/dq/f7+Ghob0wx/+UJ///Of9bhsAALHkKXxffvllnXPOOXr++ee1e/du/fjHP9bmzZv9bhsAALHk\nKXybm5t10003SZJmzpypTz75xNdGAQAQZ57C96yzzhr983PPPTcaxAAAYHJFmUwmk+sNLS0tamlp\nGXNtxYoVuvbaa/Xb3/5Wf/rTn7Ru3boxgZzN4OCQSkqKC28xAAARN2n4OmlpadErr7yiX/7yl5o6\ndeqk70+ne718jK9SqYpQtMMknjkZeOZkSOIzS9F97lSqwvE1T93OBw4c0AsvvKCNGze6Cl4AAHCa\np/BtaWnRJ598ou985zuj1zZs2KDS0lLfGgYAQFx5Ct977rlH99xzj99tAQAgEdjbGQAAwwhfAAAM\nI3wBADCM8AUAwDDCFwAAwwhfAAAMI3wBADDM0zpfAADion1Xp7a27dWhrn7VzCpTU+M8NdRVB/qZ\nhC8AILHad3Vq/Zadoz93pPtGfw4ygOl2BgAk1ta2vQ7X9wX6uYQvACCxDnX1Z71++EhfoJ9L+AIA\nEqtmVlnW63OqpgX6uYQvACCxmhrnOVyfG+jnMuEKAJBYI5Oqtrbt0+EjfZpTNU1NjXOZ7QwAQJAa\n6qoDD9vx6HYGAMAwwhcAAMMIXwAADCN8AQAwjPAFAMAwwhcAAMMIXwAADCN8AQAwjPAFAMCwokwm\nk7HdCAAAkoTKFwAAwwhfAAAMI3wBADCM8AUAwDDCFwAAwwhfAAAMS1z4vv3222psbNT27dttNyVw\njz76qJYsWaKlS5fqL3/5i+3mGPGPf/xDixYt0saNG203xZg1a9ZoyZIl+trXvqZXX33VdnMCd+LE\nCX3ve9/TbbfdpsWLFyfif8sjBgYGtGjRIr344ou2mxK49vZ2feELX9Dtt9+u22+/XQ8//LDtJvmq\nxHYDTNq/f79+85vf6PLLL7fdlMC9/fbb2rdvnzZt2qQPP/xQq1at0qZNm2w3K1D9/f16+OGH1djY\naLspxrz11lvavXu3Nm3apO7ubt188826/vrrbTcrUNu3b9dnP/tZ3XXXXTp48KC+9a1v6Utf+pLt\nZhnx1FNPqbKy0nYzjLnyyiv185//3HYzApGoyjeVSunJJ59URUWF7aYErq2tTYsWLZIkXXzxxTp2\n7JiOHz9uuVXBKi0t1a9//WvNnj3bdlOMueKKK/TEE09IkqZPn64TJ05oaGjIcquCdeONN+quu+6S\nJB0+fFjV1dWWW2TGhx9+qD179ui6666z3RT4IFHhe84556i4uNh2M4zo6urSjBkzRn+eOXOm0um0\nxRYFr6SkRGeffbbtZhhVXFyssrIySdLmzZu1cOHCxPw3vnTpUq1cuVKrVq2y3RQjVq9erfvuu892\nM4zas2eP7r77bt16663685//bLs5voptt3NLS4taWlrGXFuxYoWuvfZaSy2yi11E4+21117T5s2b\n9cwzz9huijEvvPCC/vrXv+ree+/Vli1bVFRUZLtJgXnppZd02WWX6YILLrDdFGPmzZun5cuX66tf\n/aoOHDigO+64Q6+++qpKS0ttN80XsQ3fxYsXa/HixbabYc3s2bPV1dU1+vPHH3+sVCplsUUIyptv\nvql169bp6aefTsSQygcffKCqqirNmTNH8+fP19DQkI4ePaqqqirbTQtMa2urDhw4oNbWVn300Ucq\nLS3Veeedp6uuusp20wJTXV2tG2+8UZJ04YUXatasWers7IzNF5DYhm/SXX311Vq7dq2WLl2qnTt3\navbs2SovL7fdLPist7dXa9as0bPPPqtzzz3XdnOMeOedd3Tw4EHdf//96urqUn9//5ghljh6/PHH\nR/+8du1anX/++bEOXknasmWL0um07rzzTqXTaR05ciRW4/uJOtWotbVVGzZs0D//+U/NnDlTqVQq\n1t10P/vZz/TOO++oqKhIDz74oC655BLbTQrUBx98oNWrV+vgwYMqKSlRdXW11q5dG+tQ2rRpk9au\nXauLLrpo9Nrq1atVU1NjsVXBGhgY0P3336/Dhw9rYGBAy5cv15e//GXbzTJmJHxvueUW200J1PHj\nx7Vy5Ur19PTo1KlTWr58ub74xS/abpZvEhW+AACEQaJmOwMAEAaELwAAhhG+AAAYRvgCAGAY4QsA\ngGGELwAAhhG+AAAYRvgCAGDY/wegKGGHRvHvfgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 576x396 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "2C3CoLtMbFP-",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "#### 推定のためのモデル\n",
        "データを可視化して、コレは2次元正規分布に従っていると睨んだとしましょう。すると私達が推定したいパラメータは、平均ベクトルと共分散行列になります。これらをまず初期化しておきます。"
      ]
    },
    {
      "metadata": {
        "id": "O5WouLd4beCO",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "mu_param = torch.nn.Parameter(torch.tensor([0., 0.]))\n",
        "cov_tril_param = torch.nn.Parameter(torch.tensor([[1., 0.],\n",
        "                                                  [0., 1.]]))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "IdOxDc-Nc0gP",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "#### 対数尤度関数（最適化の目的関数）\n",
        "データ $X = \\{x_1, \\cdots, x_N\\}$ に対する対数尤度関数は\n",
        "\n",
        "$$\n",
        "{\\rm LogLikelihood}(X, \\mu, \\Sigma) = \\frac{1}{N}\\sum_{i=1}^N {\\rm log}\\{{\\rm Normal}(x_i \\mid \\mu, \\Sigma)\\}\n",
        "$$\n",
        "\n",
        "であり、データ $X$ は既に手元にあるので、パラメータ $\\mu, \\Sigma$ に関して最大化するのが最尤推定法になります。続いて事前分布 $p(\\mu)$ と $p(\\Sigma)$ を目的関数に考慮するのがMAP推定になります。具体的には下記のように修正されます。\n",
        "\n",
        "$$\n",
        "{\\rm LogJointProb}(X, \\mu, \\Sigma) = \\frac{1}{N}\\sum_{i=1}^N {\\rm log}\\{{\\rm Normal}(x_i \\mid \\mu, \\Sigma)\\} + {\\rm log}p(\\mu) +{\\rm log}p(\\Sigma)\n",
        "$$\n",
        "\n",
        "今回は適当に分散の大きな正規分布を仮定し、無情報事前分布に近い状態にしておきます。\n",
        "\n",
        "#### 補足\n",
        "ちなみにベイズの定理よりパラメータの事後分布は\n",
        "\n",
        "$$\n",
        "p(\\mu, \\Sigma \\mid X) = \\frac{p(X, \\mu, \\Sigma)}{p(X)}\n",
        "$$\n",
        "\n",
        "と表され、この分母は定数になっています。したがって分子の同時分布について最大化すれば良いというのかMAP推定であり、同時分布を\n",
        "\n",
        "$$\n",
        "p(X,\\mu,\\Sigma)=p(X\\mid \\mu, \\Sigma) p(\\mu) p(\\Sigma)\n",
        "$$\n",
        "\n",
        "と表現しつつ、対数を取ることで上記の目的関数が導出されます。"
      ]
    },
    {
      "metadata": {
        "id": "Iya1lNR9b2Eu",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def log_joint_prob(mu_param, cov_param, X):\n",
        "    \n",
        "    prior_mu = torchdist.Normal(loc=torch.zeros_like(mu_param),\n",
        "                                 scale=100*torch.ones_like(mu_param))\n",
        "    prior_cov = torchdist.Normal(loc=torch.zeros_like(cov_param),\n",
        "                                 scale=100*torch.ones_like(cov_param))\n",
        "    \n",
        "    model = torchdist.MultivariateNormal(\n",
        "        loc=mu_param, \n",
        "        scale_tril=cov_param\n",
        "    )\n",
        "    \n",
        "    return (\n",
        "        model.log_prob(X).mean() \n",
        "        + prior_mu.log_prob(mu_param).mean()\n",
        "        + prior_cov.log_prob(cov_param).mean()\n",
        "    )"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Mkk01MKRd2UZ",
        "colab_type": "code",
        "outputId": "bee4a99e-d14e-472a-ad70-f9ac2a46141c",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 53
        }
      },
      "cell_type": "code",
      "source": [
        "print(\"*** initial log joint prob\")\n",
        "log_joint_prob(mu_param, cov_param, X_train)"
      ],
      "execution_count": 197,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "*** initial log joint prob\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor(-50.5738, grad_fn=<AddBackward0>)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 197
        }
      ]
    },
    {
      "metadata": {
        "id": "_bOHbf89kWuo",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "optimizer = torch.optim.Adam(params=[mu_param, cov_param], lr=1e-3)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "FyyuYqyEkw3b",
        "colab_type": "code",
        "outputId": "ac4d8582-da88-4ff6-bc67-e3bac5a2bb75",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 215
        }
      },
      "cell_type": "code",
      "source": [
        "for i in range(10000):\n",
        "    optimizer.zero_grad()\n",
        "    log_joint_prob_value = log_joint_prob(mu_param, cov_param, X_train)\n",
        "    loss_value = - log_joint_prob_value\n",
        "    loss_value.backward()\n",
        "    optimizer.step()\n",
        "    \n",
        "    if (i+1) % 1000 == 0 or (i==0):\n",
        "        print(loss_value.detach().numpy())"
      ],
      "execution_count": 199,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "50.573753\n",
            "19.344479\n",
            "16.497795\n",
            "15.430585\n",
            "14.788506\n",
            "14.359862\n",
            "14.2038765\n",
            "14.019022\n",
            "13.922331\n",
            "13.9182\n",
            "13.917614\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "w0AAZ-oPlakF",
        "colab_type": "code",
        "outputId": "fb8ab641-403c-4a8c-8ee3-a5f3edc0bfa6",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 107
        }
      },
      "cell_type": "code",
      "source": [
        "\"\"\"\n",
        "def toy_data():\n",
        "    true_mu = torch.tensor([2., 3.])\n",
        "    true_cov = torch.tensor([[2., -3.],\n",
        "                             [-3., 5.]])\n",
        "    \n",
        "    X = torchdist.MultivariateNormal(loc=true_mu, \n",
        "                                     covariance_matrix=true_cov).sample([100,])\n",
        "    return X\n",
        "\"\"\"\n",
        "\n",
        "print(\"mu map estimated: \\n\", mu_param.data)\n",
        "print(\"cov map estimated \\n\", cov_param.mm(cov_param.t()).data)"
      ],
      "execution_count": 200,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "mu map estimated: \n",
            " tensor([1.8759, 3.1613])\n",
            "cov map estimated \n",
            " tensor([[ 2.0109, -2.9852],\n",
            "        [-2.9852,  4.9609]])\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "5qLM4A32m_Ko",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### 回帰モデル\n",
        "次に同じようにTorchで回帰モデルを実装していきます。スカラーの入力 $x$ とスカラーの出力 $y$ の関係性が下記の用に多項式で表されると仮定しましょう。\n",
        "\n",
        "$$\n",
        "y = -3 + 4x + x^2 + \\epsilon\n",
        "$$\n",
        "\n",
        "ただし、ここで $\\epsilon \\sim {\\rm Normal}(0, 1)$ の正規乱数としておきます。"
      ]
    },
    {
      "metadata": {
        "id": "oRiDD1Ux1vP4",
        "colab_type": "code",
        "outputId": "367dd027-ce3f-4e48-8826-9a4c3e15d79e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 365
        }
      },
      "cell_type": "code",
      "source": [
        "def toy_poly():\n",
        "    \n",
        "    x = 5 * torch.rand(100, 1) \n",
        "    linear_op = -3 - 4*x + 1*x**2 \n",
        "    y = torchdist.Normal(linear_op, 1).sample()\n",
        "    return x, y\n",
        "\n",
        "x_train, y_train = toy_poly()\n",
        "\n",
        "plt.plot(x_train.numpy(), y_train.numpy(), \"o\")"
      ],
      "execution_count": 201,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7fdf76f67358>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 201
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFKCAYAAABcq1WoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XFsVfX9//FXaUVX2mKBS6UW6G/O\nTagh4kTs/Fm3yTSuDYnZGuo2913EjS3ZnH+4RSFBo4kb+Eu26ebYDzEuwUBTXBwJZhgZBLNgVWY0\ngHOVRGgR6y0UKSAq7f39wa+Flntu7z33c875fM55Pv6ih7b3w4fevs7nc96fz6ckk8lkBAAAQjMh\n6gYAAJA0hC8AACEjfAEACBnhCwBAyAhfAABCRvgCABCysrBeKJ0eMPr9qqvL1d9/yuj3TBr6sHj0\noRn0Y/How+KZ7sNUqtLz75wd+ZaVlUbdBOfRh8WjD82gH4tHHxYvzD50NnwBAHAV4QsAQMgIXwAA\nQkb4AgAQMsIXAICQEb4AAISM8AUAIGSELwAAIQtthysAAGzQua9XW3a9rw/6Tql2WrmaG+u1cG5N\nqG0gfAEAidG5r1d/2bx35OOe9MmRj1tu9t4O0jSmnQEAibFl1/se1w+E2g7CFwCQGB/0ZT844fCR\nk6G2g/AFACRG7bTyrNdnTJ0UajsIXwBAYjQ31ntcnx1qOyi4AgAkxnBV85ZdB3T4yEnNmDpJzY2z\nqXYGACBIC+fWhB62YxUVvqtXr9bu3bt15swZLVu2TLfeequpdgEAEFu+w/fVV19VV1eX2tvb1d/f\nrzvuuIPwBQAgD77Dd8GCBZo3b54kqaqqSp988okGBwdVWlpqrHEAAMRRSSaTyRT7Tdrb2/XGG2/o\n8ccf9/ycM2cGVVZGMAMAUHTB1csvv6xNmzbpmWeeyfl5/f3ZFzb7lUpVKp0eMPo9k4Y+LB59aAb9\nWDz6sHim+zCV8t6usqjwfeWVV7RmzRo9/fTTqqwMb09MAABc5jt8BwYGtHr1aj377LO69NJLTbYJ\nAIBY8x2+L774ovr7+3XfffeNXFu1apVqa2uNNAwAgLjyHb5LlizRkiVLTLYFAIBEYIcrAEDsdO7r\n1ZZd7+uDvlOqnVau5sb6yHe1Oh/hCwCIlc59vfrL5r0jH/ekT458bEsAc6oRACBWtux63+P6gVDb\nkQvhCwCIlQ/6su8rcfjIyZBb4o3wBQDESu208qzXJ5SUqHNfb8ityY7wBQDESnNjfdbrnw8O6S+b\n91oRwBRcAQBiZbio6pkt7+jzwaEL/n7LrgOjCq9GKqOPnFLt1HAqoxn5AgBiZ+HcGg0OZT836Pxn\nv8OV0T3pkxoayoxURgc9OiZ8AQCx5PXsd8bUSSN/jqoymvAFAMSS17Pf5sbZI3+OqjKaZ74AgFga\nfm67ZdcBHT5yUjOmTlJz4+xRz3Nrp5WrJ31h0J4/Og4C4QsAiK2Fc2tyFk81N9aP2g3r3PXZWT7b\nHMIXAJBY+YyOg0D4AgCsFcYBCcOj41SqUun0gNHv7YXwBQBYyYUDEvyi2hkAYCUXDkjwi/AFAFjJ\nhQMS/CJ8AQBWymeTDFcRvgAAK+WzSYarKLgCAFgpqmVAYSB8AQDWGm+TDFcx7QwAQMgIXwAAQkb4\nAgAQMsIXAICQEb4AAISM8AUAIGQsNQIAGBXGSUSuI3wBAMbE+SQikwhfAIAxuU4isiV8bRiZE74A\nAGNsP4ko18i85ebK0NpBwRUAwBjbTyKy5YxgRr4AAGOaG+tHjSzPXZ9txXSvLSNzwhcAYIzXSUSS\nrCjEqp1Wrp70hUEb9sic8AUAGJXtJKKV6zqzfm7YhVi5RuZhInwBAIGzZbrXljOCCV8AQOBsme6V\n7DgjmPAFAAQuqOleG4q4/CB8AQCBC2K61+XdtHyH72OPPaa33npLJSUlWr58uebNm2eyXQCAmDE9\n3evCblpefIXva6+9pgMHDqi9vV379+/X8uXL1d7ebrptAAB4sqWIyw9f4btr1y4tWrRIknTFFVfo\n448/1okTJ1RRUWG0cQCA+DD9fNamIq5C+dpesq+vT9XV1SMfT5kyRel02lijAADxMvx8tid9UkOZ\nzMjz2c59vb6/Z3Njvcf1cNfs+mGk4CqTyYz7OdXV5SorKzXxciNSqfA2wY4r+rB49KEZ9GPxbO7D\nra+/4XG9Wy03f8nX92y5uVJVVZeoY1uXunsHNLOmUq23XKmm+XW+2xlWH/oK3+nTp6uvr2/k448+\n+kipVCrn1/T3Z5+b9yuVqlQ6PWD0eyYNfVg8+tAM+rF4tvfhwQ+zt627d6Cods+pm6yV/3PdqGt+\nv5/pPswV5L6mnW+88UZt3bpVkrR3715Nnz6d570AAE+2n3YUNl8j32uvvVYNDQ1qa2tTSUmJHnro\nIdPtAgDEiC17KtvC9zPf+++/32Q7AAAxZsueyrZghysAQChs2FPZFr6e+QIAAP8IXwAAQkb4AgAQ\nMsIXAICQUXAFALCSq2f15oPwBQBYx+WzevPBtDMAwDq5zuqNA8IXAGAdl8/qzYdz084jzwCOnFLt\n1Hg9AwAAnOXyWb35cGrkO+o8yCEz50ECAOzj8lm9+XBq5JvrGQCjXwCwWyHVy3HfC9qp8I37MwAA\niCs/1ctx3gvaqWlnzoMEADfFvXq5UE6Fb9yfAQBAXDFzOZpT085xfwYAAHEV9+rlQjkVvtK5ZwCp\nVKXS6YGomwMAyENzY/2oZ77nridz5tK58AUAuIeZy9EIXwBAKOJcvVwopwquAACIA8IXAICQEb4A\nAISMZ74AgEAUsp1k0hC+AADj/GwnmSRMOwMAjGM7ydwIXwCAcWwnmRvhCwAwjoNwciN8AQDGcRBO\nbhRcAQCMYzvJ3AhfALCQrct0CmkX20l6I3wBwDJey3Q6tr+nYyc+iyyMWT5kTizC19Y7RADww2uZ\nztGBTyVFF3q5lg/xO7cwzhdcDd+J9aRPaiiTGfmh7NzXG3XTAMAXr2U6Y4W9ZpblQ+Y4H74s5AYQ\nN17LdMYKO/RYPmSO8+HLnRiAuPFapjNW2KHH8iFznH/mWzutXD3pC4OWOzEArhq7TGdyxUQdPf7p\nBZ8XduixfMgc58O3ubF+VPXduevciQFw19hlOmcLS6MPPZYPmeF8+HInBiAJCL14cT58JX4oAQBu\n8RW+Z86c0YoVK3Tw4EENDg7q17/+ta677jrTbQsUa4MBAFHxFb5///vf9YUvfEEbNmxQV1eXHnzw\nQW3atMl02wLDLi0AgCj5Ct/FixerpaVFkjRlyhQdO3bMaKOCxi4tAOKMmT37+Qrfiy66aOTPf/3r\nX0eC2BWsDQYQV8zsuWHc8O3o6FBHR8eoa7/4xS9000036bnnntPevXu1Zs2acV+ourpcZWWl/lua\nRSpV6evrZl1WqfcPH7/g+syaSt/f01VJ+/cGgT40g34sXipVqa2vv5H177a+3q2Wm78UcovcE9bP\n4bjh29raqtbW1guud3R06J///KeeeuqpUSNhL/39+e1Vmq9UqlLp9ICvr71twcysa4NvWzDT9/d0\nUTF9iLPoQzPox+IN9+HBD7P3Y3fvAH08DtM/h7mC3Ne0c3d3tzZu3Kj169fr4osv9t2wqLA2GEBc\nseufG3yFb0dHh44dO6af/OQnI9fWrVuniRMnGmtY0FgbDCCO2PXPDSWZTCYTxguZnu54p+djbdj6\nH6r5isBUX/HoQzPox+Kd34e2bEXpGuunnaNGNR8AeGNmz35OHinIGb4AAJc5OfJlnS4A+GPrBhy2\ntisoToYv1XwAUDhbH9nZ2q4gOTnt3NxY73Gdaj4A8GLrIztb2xUkJ0e+C+fWqKrqEm3Y+i7VfACQ\nJ1sf2dnariA5Gb6S1DS/TnPqJkfdDABwhq2P7GxtV5CcnHYGAJx9VrpyXafuWbVdv/g/29W5rzfn\n59v6yM7WdgXJ2ZEvACTZ2CKl9w8fH7dIydatdW1tV5AIXwBwkN9zyW3dgMPWdgWFaWcAcFASi5Ti\nhPAFAAfVTivPej3ORUpxQvgCgIO8ipS+MuvScBsCXwhfAHDQwrk1uuWrdRdc37a7Z9yqZ0SP8AUA\nR717sD/r9TjvDBUXhC8AOIqiK3ex1CggSTuhA0D4krgzVFww8g3A8OL3nvRJDWUyIyd08BwGgElJ\n3BkqLhj5BsDv4ncAKEQSd4aKC8I3ADyHARCW4Z2hUqlKpdMDUTcHeSJ8A8BzGCA+qN/Ijf7xh2e+\nAeA5DBAP1G/kRv/4x8jXp1x3ezyHAeKB+o3c6B//CF8fxh7lNXy3J2lUAPPDB7gtKfUbfqeOk9I/\nQWDa2Ydcd3sA4iMJhxcUM3WchP4JCuHrA3d7QDIkoX6jmMFEEvonKEw7+0A1M5AMSajfKGYwkYT+\nCQrh60NzY/2oZ77nrnO3B8RN3Os3ih1MxL1/gsK0sw8L59Zo2eIG1aUqVDqhRHWpCi1b3MAPIADn\nMHUcDUa+PnG3ByAOmDqOhtPhG+bOKsOvdSh9UmWlJTozlNHl0yaxmwsAo6LYMYrBRPicDd+db/aM\nu9bWlLHrej8fzAT+mgCSJ589BLzsfLNHG7b+h20eHeFs+HZs68p6fbydVfzcVXqV4uf7mgCQD787\nRhUT2oiGswVXB3uzn96Rqzze72Jyr1L8fF4TAPLld9kPG/+4x9nwnVVTmfV6rvJ4vz+gXru45POa\nAJAvvztGsfGPe5wN39Zbrsx6PVd5vN8fUK9S/HxeEwDy5XfZD9s8usfZZ75N8+t0/Pjpgsrj/S4m\nP78U/1DfCZVNmKDBoSHVTqugJB+AMX6X/bDxj3ucDV+p8PL4Yn5AKcUHEAY/v2sWzq1RVdUl2rD1\nXdbqOsLp8C0Ui8kBuMDPqoym+XWaUzc5nAaiaIkKXym8EWwUC+UBuI9lQ8lQVMFVX1+fFixYoM7O\nTlPtiYVizscE4I7Ofb1aua5T96zarpXrOo28x1k2lAxFjXxXr16tmTNnmmpLbHi9eTp2vGflaJhR\nOlC4oEaoXqsyetIndM+q7bxHY8L3yHfXrl2aNGmSvvzlL5tsTyx4vXmOHv/UutEwo3TAn6BGqLn2\nFeA9Gh++Rr6fffaZ/vSnP+mpp57SY489ltfXVFeXq6ys1M/LeUqlsm+04cfON3vUsa1LB3sHNKum\ncmQd8dhrTfPrxv1esy6r1PuHj+f1ultf71bLzV8qqu3F2Pp6t+f1KNvlEpM/h0nmWj9+cMR734Bi\n/i133naVHl+/e9zPy/Yeda0PbRRWH44bvh0dHero6Bh1rampSa2traqqqsr7hfr7c2/RWKhUqlLp\ndPYtJgs1dvro/cPHL/jhH752/Pjpcad7blswM+uSpmy6eweM/TsKlUpV6uCH2V87yna5xOTPYZK5\n2I+1U733DSjm3zKnbrKWLW4YWZUxOJTJ+nlj36Mu9qFtTPdhriAfN3xbW1vV2to66lpbW5uGhob0\n3HPP6eDBg3r77bf1hz/8QVdemX3XKZt17uvVM1v25f35+RyikG1J06nTn+vowKcXfG7UO9D43XgE\nSLogN7Y4f1XGynWdvEdjyNe088aNG0f+/MADD+iOO+5wNnjzHaEOy3ev1LFLmrxeK+odaNgZB/Bn\nvH0DTBUy8h6Np8St8z3feEcFZuP3btPWDT5sbRfgAq99A0xWQvMejaeiw/e3v/2tiXZEYryjArMp\n5m7T1i0qbW0X4Cq/5/J64T0aP4ke+Xo977yodILubp4jibtNAIXjiD+MJ9Hh6/Us5e7mOSMhS9gC\nKBSFjBhPosM36Gcp7BwFxNN4722KpDCeRIevFNyzFDZHj59sv3BbbmZTg6TJ571NkRTGk/jwDYrp\nggtEy+sXblXVJRzjljD5vrcpkkIuRZ1qBG8UXMSL52EZ27pCbQeix3sbJhC+AfHaHJ2CCzd5/cLt\n7mU7v6ThvQ0TCF+fxjvHs7mxPuvXUXDhJq9fuDNreOabNLy3YQLPfH2g4CJ5vKpXh0+/QnLw3oYJ\nhK8PFFwkj9cv3Kb5dZwkk0C8t1EswtcHCi6SiV+4AEwhfH0IY/caNugAgPgifH0IevcaNugAzuFG\nFHFE+OYh25t/2eKGwAou2KADOIsbUcRVrMPXxB2z15t/2eIGPbL0esMtPotnyv4xSooXbkQRV7Fd\n5zscmj3pkxrKZEZCc+x63PHkevMHhUX8/pj6P4c9uBFFXMU2fE2FZhRvfhbx+xPFjRKCxY0o4iq2\n4WsqNKN48y+cW6NlixtUl6pQ6YQS1aUqtGxxA9Ns42CUFD/ciCKuYvvM19RyoKjO5WRNaeE4wDx+\n2E0KcRXb8DUVmrz53cEB5nYqtgiOG1HEUWzD12Ro8uZ3AzdK9mGpEJBdbMNXIjSTiP9zu9i4VIjl\naLBBrMPXVTb8crChDXCfbUVwjMRhC8LXMqZ/OfgJUX5BwRTbiuBsHIkjmQhfy5j85eA3RPkFBVNs\nK4ILYyTOrBHyQfhaxuQvB78hattUIdxlWxFc0CNxZo2QL8LXMiZ/OfgNUdumCuG2XEVwI6PEI6dU\nOzX4UWLQI3FmjZCv2O5w5SqTO/r43Z2LXYUQhlF7cQ+Fsxd30LvHMWuEfDHyNczEhgKSmWm6XHf5\nnft6tfX1N3Tww4EL2mnbVCHiKapRYpDL0Zg1Qr4IX4NMPe8x9cvBK0QlZW3n/928V5enJo0EcZzC\nliIY+8RxlGhbgRnsRfgaZOPznmwhunJdZ9bPzSieBSIuFcEM3yQc6jupsgklOjOYGXVDFCdxHCUy\na4R8Eb4GuXIn79XO88WpQMTGm6Jsxt4kfD6YkWT3zUIx4jpKjNusEYJB+Brkyp28VzvPZ9sNQzFc\nuSnyukk49/f53Sy4MsXOKBFJRvga5MqdvFc7z2fbDUMxXLkpGm9GIp+bBZem2KVzo8RUqlLp9EDU\nzQFCw1Ijg4JexmDKcDvrZ1RpQkn2z7HthqEYriyd8loaNiyfm4VcU+wA7MHI1zBXnvcsnFujlpu/\npHR64P9PU8Z36s+V6c3xZiTyuVlwZYodSDrCF87cMBTDhX/j+TcJH/SdUOmECTozNKTLp1XkfbPg\nyhQ7kHSEL2CRYm8SXKk7AJKO8AUiZLoy2ZUpdiDpfIfvunXrtHnzZpWVlemhhx7SvHnzTLYLiL2g\nKpNdmGIHks5XtXNXV5e2bNmi559/Xo888oh27NhhuFlA/FGZDCSXr5Hv9u3bdfvtt6usrEwNDQ1q\naGgw3S4g9qhMBpLLV/geOnRIpaWlWrp0qc6cOaMHH3xQV111lem2AZ5c2cUpFyqTgeQqyWQymVyf\n0NHRoY6OjlHX+vr6dNNNN+nhhx/W7t279Zvf/EbPP/98zhc6c2ZQZWWlxbcYztj5Zo86tnXpYO+A\nZtVUqvWWK9U0v87I9318/e4Lrv/qB1818v3DEpd/B4DCjRu+2TzxxBP64he/qJaWFknSDTfcoFdf\nfTXn15jeOo7t6IoXZB+OLSYaZmLHr5XrOrOOGOtSFXpk6fVFfe9CFduHLm9wYnL2gfdz8ejD4pnu\nw1Sq0vPvfE07NzU1aePGjWppadH+/fs1Y8YM341DPAV5klCcnpW6Wpns2h7SgG18he8111yjnTt3\nasmSJZKklStXGm0U8mPzc88gA9Lvs1Kb+8s1rhzTCNjK9zrfe++9V/fee6/JtqAAto88giwm8rOL\nU9D9lSvY4xj6cZp9AKLADleOsn3kEeQ2h352cTLVX2OD9M7brtLx46c9g12S1TdJflGpDRSH8HWU\n7SOPoLc5LPRZqYn+yjZ6fnz9bk2pvDjr55/dLCN7PaMtN0nnK2SEzh7SQHEIX0e5MPKwqZjIRH95\njZ6PDnya9frhIyfltZbAlpukYYVOy7OHNFAcwtdRjDwKY6K/DvUVFpiTKyaq/OIy62+SJH/T8jbd\nXAGuIXwdxcijMCb6q2xCiT4fvHAoWzqhRINDF14/evxTzf9qKmv42naT5DUt35M+oXtWbY9NoRhg\nC8LXYYw8ClNsf53JErySNJTJaErlxVmnn989eEzLFjdYf5PkNS0vnf33xaVQDLAF4QvkcH4RUllp\n9pHv5dMq9IHHlPThIyeduEnympYfy8ZCMcBFvo4UBJJguAipJ31SQ5lM1uCVzk4h104rz/p3tj3b\n9bJwbo2WLW5QXapCpRNKPD/PtkIxwFWMfAEPXkVIF5VO0FAmoxlTJ+nO276iOXWTJcn5ArjzR+he\n+2e7cjMB2I7wBTx4FSENZTJa++tvSDq3EXvcCuCopgeCRfgipzhujZivQtcG5/Ns15X+jNvNhF+u\n/H/BPYQvPNm+f3TQTI/+XOtPFwrFguTa/xfcQsEVPOXaeCEJxhYh1aUqijqPOOn96Rr+vxAkRr7w\nZPv+0WEwOfqjP93C/xeCxMgXnlxfPmMb+tMt/H8hSIx84SmKilebC1yyta3l5sq8v76Y/rS5X+KK\nim8EifCFJ6+KV+nsOlDTQWBzgYtX26qqLhlZ5zsevxXENvdLnFHxjSARvshp7DPPIIPA1IH3QfBq\nW8e2Lq38n+vy/j5+niH77RdGy8VLesU3gsMzXxQkyApQmwtcvNrW3TsQ2Wvn6pexW2MO3yR17usN\nqpkACsDIFwUJMiBNHHgfFK+2zazJ/czXxOjTT7/YPIsAgJEvChRkBWhzY73H9egLXLza1nrLlZ5f\nY2r06adfbJ5FAED4okBBBqTpTS1M8mpb0/w6z68xNUXvp19YJgPYjWlnFCToClCbC1wKbZvJ0Weh\nr80yGcBuhC8KZnNA2iTKZ9gskwHsRvgCAYl69MlNEmAvwhcICKNPAF4IXyRa0BtRMPoEkA3hi8Ri\n20YAUWGpERKL81oBRIXwRWKxEQWAqBC+SCw2ogAQFcIXiWXzdpYA4o2CKyQWS4EARIXwRaKxFAhA\nFAhfoAg73+zRhq3/4cB6AAUhfAGfWCcMwC8KrgCfWCcMwC9GvoBPrq8TDnprTQDeCF/ApyiPDCwW\nU+ZAtHxNO/f29mrp0qW666679P3vf1979uwx3S7Aei6vE2bKHIiWr5Hvs88+q29961tqa2vTv//9\nb/3ud7/TunXrTLcNsNrCuTWqqrpEG7a+69w6YdenzAHX+Qrf6upqHTt2TJJ0/PhxVVdXG20U4Iqm\n+XWaUzc56mYUzOUpcyAOfIXvj370I333u9/VCy+8oBMnTmjDhg2m2wXAoLHFVV+ZVZ01fF2YMgfi\noCSTyWRyfUJHR4c6OjpGXWtqalJpaal+9rOfafv27Xr++ef1xz/+MecLnTkzqLKy0uJbjFjb+WaP\nOrZ16WDvgGbVVKr1livVNL8u6mY5beebPXp8/e4Lrrf87/+lPfuPqLt3QDPpayBU44ZvNvfcc4/u\nu+8+XX311frss8906623aseOHTm/Jp0e8NvGrFKpSuPfM2ls68OxFbjDli1usPY5qm19mM3KdZ1Z\nR7l1qQo9svR639/X5FIlF/rRdvRh8Uz3YSpV6fl3vqqdZ8+erbfeekuS9Pbbb2v2bKaqUDwqcIMR\nRHHV8I1ST/qkhjKZkaVKnft6fX9PIEl8PfNdtmyZVqxYoX/84x+SpBUrVhhtFJKJCtxgBFFcletG\nydZZCsAmvsJ3+vTpWrt2rem2IOGowA1Gc2N91un8YoqruFECisPezrCGy5tW2Gzh3BotW9ygulSF\nSieUqC5VUfRz9Npp5Vmvc6ME5IftJRG58wt3plReLJVIH5/4zKlNK2xn+tziIEbTQJIQvojU2Arn\nowOfSrK7whnn9n/esuuAc7t7ATYgfBEpCnfcZXo0DSQJz3wRKQp3ACQR4YtIUbgDIIkIX0SKCmcA\nScQzX0SKwh0ASUT4InIU7gBIGqadAQAIGeELAEDICF8AAEJG+AIAEDLCFwCAkFHtjEQ4//CG2mnl\nam6sp8IaQGQIX8Te2MMbetInRz4mgAFEgWlnxF6uwxsAIAqEL2KPwxsA2IbwRexxeAMA2xC+iD0O\nbwBgGwquEHsc3gDANoQvEoHDGwDYhGlnAABCRvgCABAywhcAgJARvgAAhIzwBQAgZIQvAAAhI3wB\nAAgZ4QsAQMgIXwAAQlaSyWQyUTcCAIAkYeQLAEDICF8AAEJG+AIAEDLCFwCAkBG+AACEjPAFACBk\nzoXvY489piVLlqitrU1vv/121M1x1n//+18tWrRI69evj7opzlq9erWWLFmi73znO3rppZeibo5z\nPvnkE/3yl7/UD37wA7W2tmr79u1RN8lZp0+f1qJFi/S3v/0t6qY4qbOzUzfccIPuuusu3XXXXXr0\n0UcDf82ywF/BoNdee00HDhxQe3u79u/fr+XLl6u9vT3qZjnn1KlTevTRR9XY2Bh1U5z16quvqqur\nS+3t7erv79cdd9yhW2+9NepmOWX79u26+uqr9eMf/1iHDh3S3XffrW984xtRN8tJf/7znzV58uSo\nm+G066+/Xk888URor+dU+O7atUuLFi2SJF1xxRX6+OOPdeLECVVUVETcMrdMnDhRa9eu1dq1a6Nu\nirMWLFigefPmSZKqqqr0ySefaHBwUKWlpRG3zB3f/va3R/58+PBh1dTURNgad+3fv1/vvfeevv71\nr0fdFBTAqWnnvr4+VVdXj3w8ZcoUpdPpCFvkprKyMl1yySVRN8NppaWlKi8vlyRt2rRJTU1NBK9P\nbW1tuv/++7V8+fKom+KkVatW6YEHHoi6Gc5777339NOf/lR33nmn/vWvfwX+ek6NfMdiZ0xE7eWX\nX9amTZv0zDPPRN0UZ23cuFHvvPOOfvWrX2nz5s0qKSmJuknOeOGFF3TNNddo5syZUTfFafX19fr5\nz3+u22+/Xd3d3frhD3+ol156SRMnTgzsNZ0K3+nTp6uvr2/k448++kipVCrCFiHJXnnlFa1Zs0ZP\nP/20Kisro26Oc/bs2aOpU6dqxowZmjNnjgYHB3X06FFNnTo16qY5Y8eOHeru7taOHTv04YcfauLE\nibrsssv0ta99LeqmOaWmpmbkMcisWbM0bdo09fb2BnpT41T43njjjXryySfV1tamvXv3avr06Tzv\nRSQGBga0evVqPfvss7r00kujbo6T3njjDR06dEgrVqxQX1+fTp06NeqxEsb3+9//fuTPTz75pC6/\n/HKC14fNmzcrnU5r6dKlSqfy12+wAAAAvUlEQVTTOnLkSOA1CE6F77XXXquGhga1tbWppKREDz30\nUNRNctKePXu0atUqHTp0SGVlZdq6dauefPJJQqQAL774ovr7+3XfffeNXFu1apVqa2sjbJVb2tra\ntGLFCn3ve9/T6dOntXLlSk2Y4FQZCmLim9/8pu6//35t27ZNn3/+uR5++OFAp5wljhQEACB03GYC\nABAywhcAgJARvgAAhIzwBQAgZIQvAAAhI3wBAAgZ4QsAQMgIXwAAQvb/AEloGoE6zvHWAAAAAElF\nTkSuQmCC\n",
            "text/plain": [
              "<Figure size 576x396 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "G4wncEBU2n3x",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "#### モデル\n",
        "さて、上記のデータを見て、二次関数で表現できると睨んだとしましょう（出来レースですが）。\n",
        "\n",
        "$$\n",
        "y = w_0 + w_1x + w_2x^2 + \\epsilon\n",
        "$$\n",
        "\n",
        "すると、簡単な式変形から下記のような確率モデルを使うことが考えられます（$\\epsilon$の分散は既知としているが、未知にすることも当然必要であれば考えられる）。 \n",
        "\n",
        "\n",
        "$$\n",
        "y - (w_0 + w_1x + w_2x^2) \\sim {\\rm Normal}(0,1)\n",
        "$$\n",
        "\n",
        "あるいは、分布の平均の方にパラメータを持つ項を吸収させて\n",
        "\n",
        "$$\n",
        "y \\sim {\\rm Normal}(w_0 + w_1x + w_2x^2 ,1)\n",
        "$$\n",
        "\n",
        "とできます。事前分布を適当に置いてしまえば、MAP推定に必要な同時分布 $p(w_0, w_1, w_2, x, y)$ が一先ず書き下せますので、その対数値を最大化することで回帰問題を解くことが可能になります。今回も各パラメータに対してはそれぞれ適当な正規分布を仮定してしまいましょう。\n",
        "\n",
        "$$\n",
        "{\\rm LogJointProb}(w_0, w_1, w_2, X, Y) = \\frac{1}{N} \\sum_{i=1}^N {\\rm Normal}(y_i\\mid w_0 + w_1x_i + w_2x_i^2 ,1) +{\\rm log} p(w_0) +{\\rm log} p(w_1) +{\\rm log} p(w_2)\n",
        "$$"
      ]
    },
    {
      "metadata": {
        "id": "LmDgr1854-ZS",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "w0 = torch.nn.Parameter(torch.tensor(1.))\n",
        "w1 = torch.nn.Parameter(torch.tensor(1.))\n",
        "w2 = torch.nn.Parameter(torch.tensor(1.))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "yM3nPHcf5J25",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def log_joint_prob(w0, w1, w2, x, y):\n",
        "    \n",
        "    prior_w0 = torchdist.Normal(torch.tensor(0.), 10*torch.tensor(1.))\n",
        "    prior_w1 = torchdist.Normal(torch.tensor(0.), 10*torch.tensor(1.))\n",
        "    prior_w2 = torchdist.Normal(torch.tensor(0.), 10*torch.tensor(1.))\n",
        "\n",
        "    linear = w0 + w1*x + w2*x**2\n",
        "    likelihood = torchdist.Normal(linear, torch.ones_like(linear))\n",
        "    \n",
        "    return (\n",
        "        prior_w0.log_prob(w0).mean() +\n",
        "        prior_w1.log_prob(w1).mean() +\n",
        "        prior_w2.log_prob(w2).mean() +\n",
        "        likelihood.log_prob(y).mean()\n",
        "    )    \n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "DcnPsxO-T9Ie",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 575
        },
        "outputId": "a998b685-37ae-430d-8290-1370fe6222cd"
      },
      "cell_type": "code",
      "source": [
        "optimizer = torch.optim.Adam(params=[w0, w1, w2], lr=1e-3)\n",
        "\n",
        "for i in range(30000):\n",
        "    optimizer.zero_grad()\n",
        "    log_joint_prob_value = log_joint_prob(w0, w1, w2, x_train, y_train)\n",
        "    loss_value = - log_joint_prob_value\n",
        "    loss_value.backward()\n",
        "    optimizer.step()\n",
        "    \n",
        "    if (i+1) % 1000 == 0 or (i==0):\n",
        "        print(loss_value.detach().numpy())"
      ],
      "execution_count": 256,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "11.186682\n",
            "11.186682\n",
            "11.186681\n",
            "11.186681\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186681\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n",
            "11.186682\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "vfoligfvUrqV",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 71
        },
        "outputId": "7157e7c9-4598-4758-9fa1-39feae34d626"
      },
      "cell_type": "code",
      "source": [
        "\"\"\"\n",
        "def toy_poly():\n",
        "    \n",
        "    x = 5 * torch.rand(100, 1) \n",
        "    linear_op = -3 - 4*x + 1*x**2 \n",
        "    y = torchdist.Normal(linear_op, 1).sample()\n",
        "    return x, y\n",
        "\n",
        "\"\"\"\n",
        "print(w0.data)\n",
        "print(w1.data)\n",
        "print(w2.data)"
      ],
      "execution_count": 257,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "tensor(-2.5727)\n",
            "tensor(-3.9757)\n",
            "tensor(0.9658)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "Hh15BgdT-v3M",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### 変分推論\n",
        "ベイズ推論は事後分布 \n",
        "\n",
        "$$\n",
        "p(\\theta \\mid D) = \\frac{p(D\\mid \\theta)p(\\theta)}{p(D)}\n",
        "$$\n",
        "\n",
        "において、同時分布の最大化に甘んじず（すなわち関数の一番山となる部分だけを探すのではなく）、分布の形状全体を把握しようという試みになります。その試みは一般に困難を極めます。都合の良い尤度関数と都合の良い事前分布を選ばない限りは形状全体を上手に求めることはできません。\n",
        "\n",
        "したがって、形状全体を知りたいのだが、ある程度簡略化した形状で一番近いものを探せればいいというのが変分推論です。変分モデル $q(\\theta; \\eta)$ を仮定し（$\\eta$は変分パラメータと呼ばれる最適化すべきパラメータである）、$\\eta$ の調整で分布の形状を変えること $p(\\theta |D)$ に最も近い $q(\\theta; \\eta)$ を決定します。近いというのはKLダイバージェンスの意味であり\n",
        "\n",
        "$$\n",
        "{\\rm KL}[q(\\theta; \\eta) : p(\\theta\\mid D)] = {\\mathbb E}_{q(\\theta;\\eta)}[{\\rm log}q(\\theta; \\eta)] - {\\mathbb E}_{q(\\theta;\\eta)}[{\\rm log}p(\\theta)] - {\\mathbb E}_{q(\\theta;\\eta)}[{\\rm log}p(D\\mid\\theta)]\n",
        "$$\n",
        "\n",
        "を最小化するような $\\eta$ を求めます。\n",
        "\n",
        "こちらも最適化問題ではありますが、求めているものはパラメータ $\\theta $ の値ではなくパラメータ $\\theta$ が取りうる値の分布を網羅的に把握するために $\\eta$ を最適化していることに注意しましょう。最適化された $\\eta$ によって分布 $q(\\theta ; \\eta)$ が定まり、この分布からサンプリングをしたりすることで、単に点推定で $\\theta$ を決めてしまうよりも多くの情報を利用することができるというわけです。\n",
        "\n",
        "実際の推論では期待値計算の代わりに現在の $\\eta$ の値を用いて \n",
        "\n",
        "$$\n",
        "\\theta^* \\sim q(\\theta; \\eta)\n",
        "$$\n",
        "\n",
        "とサンプリングし、サンプリングされた$\\theta^*$ で現在のKLダイバージェンスを計算するということにします（そんなのいい加減すぎる！と思うのであれば、$q(\\theta; \\eta)$ は現在の $\\eta$ を使うとして重点サンプリングなどをしてもいいだろうし、大げさにもMCMCを使ってもいいだろう。単に計算量の問題である）。\n",
        "\n",
        "$$\n",
        "{\\rm KL}[q(\\theta; \\eta) : p(\\theta\\mid D)] \\simeq {\\rm log}q(\\theta^*; \\eta) - {\\rm log}p(\\theta^*) - {\\rm log}p(D\\mid\\theta^*)\n",
        "$$\n",
        "\n",
        "\n",
        "#### 変分モデル\n",
        "回帰問題の例に戻って、パラメータ $w_0, w_1, w_2$ に対してそれぞれ変分モデル\n",
        "\n",
        "$$\n",
        "q(w_i ; \\eta_i) = {\\rm Normal}(\\mu_i, \\sigma_i)\n",
        "$$\n",
        "\n",
        "を仮定しましょう。すなわち各 $w_i$ に対して正規分布を仮定して、あとはそれぞれの平均分散を変分パラメータとして最適化して $w_i$ の分布を得てしまおうということにしたのです。"
      ]
    },
    {
      "metadata": {
        "id": "n4xTF9EUCD-D",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "variational_params = {\n",
        "    \"w0_loc\": torch.nn.Parameter(torch.tensor(0.)),\n",
        "    \"w0_scale\": torch.nn.Parameter(torch.tensor(1.)),\n",
        "    \"w1_loc\": torch.nn.Parameter(torch.tensor(0.)),\n",
        "    \"w1_scale\": torch.nn.Parameter(torch.tensor(1.)),\n",
        "    \"w2_loc\": torch.nn.Parameter(torch.tensor(0.)),\n",
        "    \"w2_scale\": torch.nn.Parameter(torch.tensor(1.)),\n",
        "}\n",
        "\n",
        "def variational_model(variational_params):\n",
        "    \"\"\"\n",
        "    Variational model q(w; eta)\n",
        "    arg: variational parameters \"eta\"\n",
        "    return: w ~ q(w; eta)\n",
        "    \"\"\"\n",
        "    w0_q = torchdist.Normal(\n",
        "        variational_params[\"w0_loc\"],\n",
        "        torch.exp(variational_params[\"w0_scale\"]),\n",
        "    )\n",
        "    \n",
        "    w1_q = torchdist.Normal(\n",
        "        variational_params[\"w1_loc\"],\n",
        "        torch.exp(variational_params[\"w1_scale\"]),\n",
        "    )\n",
        "    \n",
        "    w2_q = torchdist.Normal(\n",
        "        variational_params[\"w2_loc\"],\n",
        "        torch.exp(variational_params[\"w2_scale\"]),\n",
        "    )\n",
        "    \n",
        "    return w0_q, w1_q, w2_q"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "huYHejDTJOSy",
        "colab_type": "code",
        "outputId": "5802db92-0378-4031-9a2e-bb268d66e554",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 89
        }
      },
      "cell_type": "code",
      "source": [
        "print(\"*** w0, w1, w2 :variational model with initial variatinal params\")\n",
        "variational_model(variational_params)"
      ],
      "execution_count": 259,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "*** w0, w1, w2 :variational model with initial variatinal params\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(Normal(loc: 0.0, scale: 2.7182817459106445),\n",
              " Normal(loc: 0.0, scale: 2.7182817459106445),\n",
              " Normal(loc: 0.0, scale: 2.7182817459106445))"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 259
        }
      ]
    },
    {
      "metadata": {
        "id": "UrU74P8QJdKv",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### KLダイバージェンス\n",
        "変分モデルを書き終えたので、後は変分モデルから出てくるサンプルを引数としたKLダイバージェンスを書き下せばいいです。KLダイバージェンスの近似は実は下記の通り、MAP推定でも利用してきたLog Joint Prob（データ $D$とパラメータ $\\theta$の対数同時確率）\n",
        "と新たに出てきたVariational model （勝手に仮定したパラメータ $\\theta $ の事後分布）の対数確率値で構成されています。\n",
        "\n",
        "$$\n",
        "\\begin{align}\n",
        "{\\rm KL}[q(\\theta; \\eta) : p(\\theta\\mid D)] &\\simeq {\\rm log}q(\\theta^*; \\eta) - {\\rm log}p(\\theta^*) - {\\rm log}p(D\\mid\\theta^*) \\\\\\ \n",
        "& = {\\rm log}q(\\theta^*; \\eta) - {\\rm LogJointProb}(\\theta^*, D)\n",
        "\\end{align}\n",
        "$$\n",
        "\n",
        "ですので、実装上はMAP推定に使っていた目的関数をそのまま流用できます。\n",
        "\n",
        "すなわち実装は、$\\eta$ の関数である変分モデルを作り、パラメータ $\\theta$ をサンプリングできるようにして、そのサンプリング値$\\theta^*$をMAP推定で使っていた目的関数に代入してやれば良いのです。更に、サンプリングされた $\\theta^*$の$q(\\theta; \\eta)$ におけるの対数確率値も追加してやれば、目的関数の作成は終了です。"
      ]
    },
    {
      "metadata": {
        "id": "fkMzu9zJNNPp",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def kl_divergence(variational_params, x, y):\n",
        "    w0_q, w1_q, w2_q = variatinal_model(variational_params)\n",
        "    \n",
        "    w0_sample = w0_q.sample()\n",
        "    w1_sample = w1_q.sample()    \n",
        "    w2_sample = w2_q.sample()\n",
        "    \n",
        "    log_joint_prob_value = log_joint_prob(w0_sample, w1_sample, w2_sample, x, y)\n",
        "    log_variational_prob_value = (\n",
        "        w0_q.log_prob(w0_sample) +\n",
        "        w1_q.log_prob(w1_sample) +\n",
        "        w2_q.log_prob(w2_sample)\n",
        "    )\n",
        "    \n",
        "    return log_variational_prob_value - log_joint_prob_value"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "-eCSuQ4cP2ua",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 575
        },
        "outputId": "e695af16-c1cd-4e15-86c1-8a7102430eb0"
      },
      "cell_type": "code",
      "source": [
        "optimizer = torch.optim.Adam(params=variational_params.values(), lr=1e-6)\n",
        "\n",
        "for i in range(30000):\n",
        "    optimizer.zero_grad()\n",
        "    loss_value =kl_divergence(variational_params, x_train, y_train)\n",
        "    loss_value.backward()\n",
        "    optimizer.step()\n",
        "    \n",
        "    if (i+1) % 1000 == 0 or (i==0):\n",
        "        print(loss_value.detach().numpy())"
      ],
      "execution_count": 273,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "35.93248\n",
            "91.161766\n",
            "17.718224\n",
            "18.213715\n",
            "19.38914\n",
            "15.932632\n",
            "14.892603\n",
            "16.744282\n",
            "22.030283\n",
            "14.744864\n",
            "18.736307\n",
            "15.626644\n",
            "20.092701\n",
            "21.637543\n",
            "16.463593\n",
            "38.320213\n",
            "15.833952\n",
            "33.40632\n",
            "32.47555\n",
            "62.660324\n",
            "18.766312\n",
            "34.779316\n",
            "13.467247\n",
            "38.219856\n",
            "17.71047\n",
            "36.875412\n",
            "18.367641\n",
            "21.203026\n",
            "24.813643\n",
            "49.102386\n",
            "23.804152\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "-z6gnrVnWD8C",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 125
        },
        "outputId": "f630fa1b-4d17-4001-ec8a-470a38aaa107"
      },
      "cell_type": "code",
      "source": [
        "for k, v in variational_params.items():\n",
        "    print(\"{} : {}\".format(k, v))"
      ],
      "execution_count": 274,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "w0_loc : -0.009247524663805962\n",
            "w0_scale : 0.9921383857727051\n",
            "w1_loc : -0.007853412069380283\n",
            "w1_scale : 1.0310311317443848\n",
            "w2_loc : 0.0171018298715353\n",
            "w2_scale : 1.0298856496810913\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "VBfVH_TZYnIS",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}