0_stan_quick_intro_sergey.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "0-stan-quick-intro-sergey",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "lLupYdeZ4gqg",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import pystan\n",
        "import numpy as np\n",
        "\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x1l_8BtK54ic",
        "colab_type": "text"
      },
      "source": [
        "# Regression with normal distribution.\n",
        "\n",
        "Data Types\n",
        "\n",
        "\n",
        "\n",
        "* Basic types: real, int, vector, row_vector, matrix\n",
        "* Constrained types: simplex, unit_vector, ordered, positive_ordered, corr_matrix, cov_matrix\n",
        "\n",
        "Bounded Variables\n",
        "* Can put upper and lower bounds on variables\n",
        "* Ex: real<lower=0> sigma;\n",
        "\n",
        "Program Blocks\n",
        "* data (input data we want to fit)\n",
        "* transformed data (functions of your data)\n",
        "* parameters (parameters of your model)\n",
        "* transformed parameters (functions of your parameters)\n",
        "* model (defines priors of your parameters and \"loss\" function)\n",
        "* generated quantities (some additional quantities you want to output)\n",
        "\n",
        "The only required program block is the model block."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5aAH55ze4n3v",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "script = \"\"\"\n",
        "data {\n",
        "    int<lower=1> N;\n",
        "    vector[N] y;\n",
        "    vector[N] x;\n",
        "}\n",
        "parameters {\n",
        "    real alpha;\n",
        "    real beta;\n",
        "    \n",
        "}\n",
        "transformed parameters {\n",
        "  vector[N] mu;\n",
        "\n",
        "  mu = alpha + x * beta;\n",
        "}\n",
        "model {\n",
        "  // the following defines priors \n",
        "  alpha ~ normal(0, 1000);\n",
        "  beta ~ normal(0, 1000);\n",
        "  // the following defines \"loss\"\n",
        "  y ~ normal(mu, 1);\n",
        "}\n",
        "\"\"\"\n"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GU9EJ2NW5W8h",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "0b6e9e1f-e6b8-4992-92c4-2e0855badc8d"
      },
      "source": [
        "sm = pystan.StanModel(model_code=script) # this may take up to 1 minute"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_821eaf683878c5e35d7a262b03476e19 NOW.\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "g_arM58AZRsm",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 405
        },
        "outputId": "b5d07054-54eb-4189-c903-ee5c4436f451"
      },
      "source": [
        "n_train = 20\n",
        "x_train = np.random.normal(size=n_train)\n",
        "y_train = 3*x_train + 1\n",
        "\n",
        "fit = sm.sampling(data=dict(y=y_train, N=len(y_train), x=x_train))\n",
        "s = fit.extract()\n",
        "s"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "OrderedDict([('alpha',\n",
              "              array([0.61508944, 0.95163289, 1.00954125, ..., 0.72609786, 0.71219839,\n",
              "                     1.02309446])),\n",
              "             ('beta',\n",
              "              array([2.12775344, 2.23503874, 3.16301888, ..., 2.62063453, 3.08662956,\n",
              "                     3.39397544])),\n",
              "             ('mu',\n",
              "              array([[ 0.67253152,  0.36923411,  2.70192768, ...,  2.50370631,\n",
              "                       2.05304057, -2.35013129],\n",
              "                     [ 1.01197131,  0.69338108,  3.14369342, ...,  2.93547735,\n",
              "                       2.4620882 , -2.16309981],\n",
              "                     [ 1.09493196,  0.64406421,  4.11173751, ...,  3.81707087,\n",
              "                       3.14713223, -3.39841739],\n",
              "                     ...,\n",
              "                     [ 0.79684604,  0.42329164,  3.29633943, ...,  3.05220129,\n",
              "                       2.49714158, -2.92599816],\n",
              "                     [ 0.79552685,  0.35554789,  3.73947424, ...,  3.45192403,\n",
              "                       2.79816492, -3.58930469],\n",
              "                     [ 1.11472021,  0.6309311 ,  4.35180615, ...,  4.03562361,\n",
              "                       3.31676755, -3.70672346]])),\n",
              "             ('lp__',\n",
              "              array([-4.02219774, -2.21815783, -0.10072627, ..., -1.18686993,\n",
              "                     -0.88409105, -0.58828393]))])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "o0kKiJlxhbax",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 152
        },
        "outputId": "d6d4856d-bfcc-4be9-e000-3a327d3f38c4"
      },
      "source": [
        "s['mu'].mean(0), y_train"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([ 1.08521483,  0.6574899 ,  3.94716915, -1.76399451,  1.40875988,\n",
              "         1.04567583, -1.98804073,  0.46557258, -0.44164262, -0.7919078 ,\n",
              "         0.66844152,  2.07012121,  2.88185076,  0.21292117, -0.17503042,\n",
              "         2.23931928,  1.77548164,  3.66762761,  3.03207661, -3.177493  ]),\n",
              " array([ 1.08098976,  0.65335928,  3.94231212, -1.76759043,  1.40446337,\n",
              "         1.04145949, -1.99158718,  0.46148434, -0.44553054, -0.79571837,\n",
              "         0.66430848,  2.06567866,  2.87722897,  0.20888871, -0.1789772 ,\n",
              "         2.23483937,  1.77110415,  3.66283231,  3.02742164, -3.1807768 ]))"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_7HPnOYoCoEt",
        "colab_type": "text"
      },
      "source": [
        "# Exercise\n",
        "\n",
        "Modify the previous script so that it has additional parameter on variance. Make a prior to that parameter from exponetial distribution (https://mc-stan.org/docs/2_20/functions-reference/exponential-distribution.html). What parameter sigma did it learn compared to the previous constant value. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_E8tZ8_0Cm8v",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "a674dd0b-a857-4a3b-e8d5-62b61970de57"
      },
      "source": [
        "script = \"\"\"\n",
        "data {\n",
        "    int<lower=1> N;\n",
        "    vector[N] y;\n",
        "    vector[N] x;\n",
        "}\n",
        "parameters {\n",
        "    real alpha;\n",
        "    real beta;\n",
        "    // replace next line with your code\n",
        "    // <define parameter sigma>\n",
        "    \n",
        "}\n",
        "transformed parameters {\n",
        "  vector[N] mu;\n",
        "\n",
        "  mu = alpha + x * beta;\n",
        "}\n",
        "model {\n",
        "  // the following defines priors \n",
        "  alpha ~ normal(0, 1000);\n",
        "  beta ~ normal(0, 1000);\n",
        "  // <define prior for parameter sigma>\n",
        "  // <define the loss with sigma>\n",
        "}\n",
        "\"\"\"\n",
        "\n",
        "sm = pystan.StanModel(model_code=script) # this may take up to 1 minute"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_3be6059edddab6ab4e44ab38c965d11e NOW.\n"
          ],
          "name": "stderr"
        }
      ]
    }
  ]
}