hmm-gibbs-sampling.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "NLPhw4.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lNPnXo6Vrw-z"
      },
      "source": [
        "# NLP HW 4: Hidden Markov Model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mSJPMjXWDVeU"
      },
      "source": [
        "## Setup\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WaO6nRdwZ5qW"
      },
      "source": [
        "Downloading data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Kj46rvwvoci3"
      },
      "source": [
        "# from https://medium.com/@Keshav31/colab-features-download-and-upload-e1ec537a83df\n",
        "from urllib.request import urlretrieve\n",
        "import os\n",
        "from zipfile import ZipFile\n",
        "\n",
        "url = 'https://ttic.uchicago.edu/~kgimpel/teaching/31210-s19/data/31210-s19-hw3.zip'\n",
        "file = '31210-s19-hw3.zip'\n",
        "\n",
        "if not os.path.isfile(file):\n",
        "    urlretrieve(url,file)\n",
        "\n",
        "with ZipFile(file) as zipf:\n",
        "    zipf.extractall()\n",
        "\n",
        "!rm -rf data/\n",
        "!mkdir data/\n",
        "!mv 31210-s19-hw3/* data/\n",
        "!rm -rf 31210-s19-hw3.zip 31210-s19-hw3/"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a6SG96QHZ07o"
      },
      "source": [
        "Loading training data and doing bigram counts"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UwButzrDrhp3",
        "outputId": "cae219f2-9501-4343-d8d4-07e8ecb30568",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 69
        }
      },
      "source": [
        "from collections import Counter\n",
        "    \n",
        "word_tag_counts = Counter()\n",
        "tag_bigram_counts = Counter()\n",
        "vocab = set()\n",
        "pos_tags = set()\n",
        "    \n",
        "start_token = ('<s>', '<s>')\n",
        "end_token = ('</s>', '</s>')\n",
        "    \n",
        "with open('data/en_ewt.train', 'r') as trainfile:\n",
        "    prev = start_token\n",
        "    for row in trainfile.readlines():\n",
        "        row = row.split()\n",
        "        if len(row) == 0:\n",
        "            row = end_token\n",
        "        elif prev[0] == '</s>': \n",
        "            prev = start_token\n",
        "        if row[0] not in ['<s>', '</s>']:\n",
        "            vocab.add(row[0])\n",
        "            pos_tags.add(row[1])\n",
        "        word_tag_counts[(row[0], row[1])] += 1\n",
        "        tag_bigram_counts[(prev[1], row[1])] += 1\n",
        "        prev = row\n",
        "\n",
        "print('data loaded')\n",
        "print('vocab size:', str(len(vocab)))\n",
        "print('num pos tags:', str(len(pos_tags)))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "data loaded\n",
            "vocab size: 19380\n",
            "num pos tags: 50\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v2g_Z04wfgIc"
      },
      "source": [
        "Compute probabilities"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pshhjt7tffOX"
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "log_p_tau = {}\n",
        "lambda_tau = 0.1\n",
        "\n",
        "for y_p in list(pos_tags):\n",
        "    denom = sum([tag_bigram_counts[(y_p, y_pp)] + lambda_tau\n",
        "                       for y_pp in list(pos_tags) + ['</s>']])\n",
        "    \n",
        "    for y in list(pos_tags) + ['</s>']:\n",
        "        numer = tag_bigram_counts[(y_p, y)]\n",
        "        log_p_tau[(y, y_p)] = np.log((numer + lambda_tau)/denom)\n",
        "\n",
        "for y in list(pos_tags):\n",
        "    denom = sum([tag_bigram_counts[('<s>', y_pp)] + lambda_tau\n",
        "                        for y_pp in list(pos_tags)])\n",
        "    numer = tag_bigram_counts[('<s>', y)]\n",
        "    log_p_tau[(y, '<s>')] = np.log((numer + lambda_tau)/denom)\n",
        "\n",
        "log_p_eta = {}\n",
        "lambda_eta = 0.001\n",
        "for y in list(pos_tags):\n",
        "    denom = sum([word_tag_counts[(x_p, y)] + lambda_eta\n",
        "                       for x_p in list(vocab)])\n",
        "    for x in list(vocab):\n",
        "        numer = word_tag_counts[(x, y)]\n",
        "        log_p_eta[(x, y)] = np.log((numer + lambda_eta)/denom)\n",
        "\n",
        "# make sure that all the zero probs are right\n",
        "log_p_tau[('</s>', '<s>')] = -float('inf')\n",
        "for y in list(pos_tags):\n",
        "    log_p_tau[('<s>', y)] = -float('inf')\n",
        "    log_p_tau[(y, '</s>')] = -float('inf')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "avKga_pnDQsv"
      },
      "source": [
        "Functions to calculate log-probability and accuracy of tagged sentence"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "65zvhF9sDxRj"
      },
      "source": [
        "# takes a list of (word, tag) tuples\n",
        "# should start with <s>\n",
        "def log_sentence_probability(s):\n",
        "    total = log_p_tau[('</s>', s[-1][1])]\n",
        "    for i in range(1, len(s)):\n",
        "        total += log_p_tau[(s[i][1], s[i-1][1])]\n",
        "        total += log_p_eta[(s[i][0], s[i][1])] \n",
        "    return total\n",
        "\n",
        "def calc_accuracy_rate(sentences, gold_standards):\n",
        "    total = 0\n",
        "    total_correct = 0\n",
        "    assert(len(sentences) == len(gold_standards))\n",
        "    for s_num in range(len(sentences)):\n",
        "        for w_num in range(len(sentences[s_num])):\n",
        "            total += 1\n",
        "            if sentences[s_num][w_num][1] == gold_standards[s_num][w_num][1]:\n",
        "                total_correct += 1\n",
        "    return total_correct / total"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IRK9vettmizq"
      },
      "source": [
        "Reusable function to score prediction functions"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k4UyKw7Kr3gA"
      },
      "source": [
        "import time, tqdm\n",
        "\n",
        "# start_work_func prepares all memoization\n",
        "#     takes no args, returns nothing\n",
        "# pred_func does predictions\n",
        "#     takes a sentence of [(word, gs_tag)]\n",
        "#     returns predictions [(word, pred_tag)]\n",
        "def score_prediction_function(start_work_func, pred_func, do_setup=True):\n",
        "    start_time = time.time()\n",
        "    if do_setup:\n",
        "        print('doing setup, building lookup tables...', flush=True)\n",
        "        start_work_func()\n",
        "    print('making predictions...', flush=True)\n",
        "    with open('data/en_ewt.dev', 'r') as trainfile:\n",
        "        current_sent = [start_token]\n",
        "        predictions = []\n",
        "        gold_standards = []\n",
        "        total_log_prob = 0\n",
        "        lines = trainfile.readlines()\n",
        "        for row in tqdm.tqdm(lines):\n",
        "            row = row.split()\n",
        "            if len(row) == 0:\n",
        "                gold_standards.append(current_sent)\n",
        "                pred = pred_func(current_sent)\n",
        "                log_prob = log_sentence_probability(pred)\n",
        "                predictions.append(pred)\n",
        "                total_log_prob += log_prob\n",
        "                current_sent = [start_token]\n",
        "            else:\n",
        "                current_sent.append(row)\n",
        "\n",
        "    end_time = time.time()\n",
        "    elapsed_time = end_time - start_time\n",
        "\n",
        "    print()\n",
        "    tagging_accuracy = calc_accuracy_rate(predictions, gold_standards)\n",
        "    print('tagging accuracy on DEV:', tagging_accuracy)\n",
        "    print('time required for predictions (s):', elapsed_time)\n",
        "    print('log-probability of predictions:', total_log_prob, flush=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MYp2zKI0sAAd"
      },
      "source": [
        "## 1. Gibbs Sampling for HMMs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7vMvf1qJsFM4"
      },
      "source": [
        "### a) Gibbs Sampling Derivation\n",
        "\n",
        "$P(Y_t=y|Y_{-t}=y_{-t}, \\textbf{X}=\\textbf{x})$\n",
        "\n",
        "by HMM independence rules and $A\\bot B \\rightarrow P(A|B)=P(A)$:\n",
        "\n",
        "$=P(Y_t=y|Y_{t-1}=y_{t-1}, Y_{t+1}=y_{t+1}, X_t=x_t)$\n",
        "\n",
        "by Bayes' theorem:\n",
        "\n",
        "$=\\frac{P(Y_{t-1}=y_{t-1}, Y_{t+1}=y_{t+1}, X_t=x_t|Y_t=y) P(Y_t=y)}{P(Y_{t-1}=y_{t-1}, Y_{t+1}=y_{t+1}, X_t=x_T)}$\n",
        "\n",
        "by HMM independence rules, cancelling constant W.R.T. $Y_t$, and rearranging:\n",
        "\n",
        "$\\propto (P(Y_t=y)P(Y_{t-1}=y_{t-1}|Y_t=y))P(X_t=x_t|Y_t=y)P(Y_{t+1}=y_{t+1}|Y_t=y_t)$\n",
        "\n",
        "one more Bayes' theorem:\n",
        "\n",
        "$=(P(Y_{t-1}=y_{t-1})P(Y_t=y|Y_{t-1}=y_{t-1}))P(X_t=x_t|Y_t=y)P(Y_{t+1}=y_{t+1}|Y_t=y_t)$\n",
        "\n",
        "cancel one more constant, and done!\n",
        "\n",
        "$\\propto P(Y_t=y|Y_{t-1}=y_{t-1})P(X_t=x_t|Y_t=y)P(Y_{t+1}=y_{t+1}|Y_t=y_t)$\n",
        "\n",
        "To normalize with log-probabilities, plug in every possible value of $y$ to the formula above, then `np.logaddexp` them all together. Add the three terms above for a given $y$ then subtract the normalization factor to get the normalized log-probability."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "czV6NmLEuTh3"
      },
      "source": [
        "### b) Gibbs Sampling Special Cases\n",
        "\n",
        "I'm assuming you just want us to get these into a form where we can use the lookup tables for $p_\\tau$ and $p_\\eta$, so that's all the simplification that I bothered with. Normalization is done the same way as in the last part.\n",
        "\n",
        "---\n",
        "\n",
        "$P(Y_1=y | Y_{-1}=y_{-1}, \\textbf{X}=\\textbf{x})$\n",
        "\n",
        "Let `<s>` stand in for $y_0$. $P(Y_0=\\texttt{<s>})=1$ regardless of any other $\\textbf{X}, \\textbf{Y}$, so $Y_0 \\bot Y_1$, so we can just add the condition of $Y_0$'s value. Plug in to the previous formula.\n",
        "\n",
        "$\\propto P(Y_1=y|Y_0=\\texttt{<s>})P(X_1=x_1|Y_1=y)P(Y_2=y_2|Y_1=y_1)$\n",
        "\n",
        "$\\square$\n",
        "\n",
        "---\n",
        "\n",
        "$P(Y_T=y | Y_{-T}=y_{-T}, \\textbf{X}=\\textbf{x})$\n",
        "\n",
        "Let `</s>` stand in for $y_{t+1}$. $P(Y_{t+1}=\\texttt{</s>})=1$ regardless of any other $\\textbf{X}, \\textbf{Y}$, so $Y_T \\bot Y_{T+1}$, so we can just add the condition of $Y_{T+1}$'s value. Plug in to the previous formula.\n",
        "\n",
        "$\\propto P(Y_T=y|Y_{T-1}=y_{T-1})P(X_T=x_T|Y_T=y)P(Y_{T+1}=\\texttt{</s>}|Y_T=y_T)$\n",
        "\n",
        "$\\square$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N4bAx21J1bWZ"
      },
      "source": [
        "### c) Implement Gibbs sampling"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vDAQuuW5Wc9N"
      },
      "source": [
        "import functools, random\n",
        "\n",
        "# make sure things stay ordered properly for prob dists\n",
        "l_pos_tags = list(pos_tags)\n",
        "\n",
        "# computes log gibbs probability for values of variables above\n",
        "log_gibbs_prob = lambda yt, ytm, ytp, xt: \\\n",
        "    (log_p_tau[(yt, ytm)] + log_p_eta[(xt, yt)] + log_p_tau[(ytp, yt)])\n",
        "\n",
        "log_gibbs_probs = lambda *args:\\\n",
        "    [log_gibbs_prob(yt, *args) for yt in l_pos_tags]\n",
        "\n",
        "# computes log normalization factor to be subtracted from above\n",
        "log_gibbs_probs_normed = lambda probs:\\\n",
        "    probs - functools.reduce(np.logaddexp, probs)\n",
        "\n",
        "probs_memo = {}\n",
        "\n",
        "gibbs_setup = lambda: None\n",
        "\n",
        "def gibbs_iteration(Y_T, sentence):\n",
        "    global probs_memo, beta, beta_update\n",
        "    beta_update()\n",
        "    Y_T1 = ['<s>']\n",
        "    for ix in range(1, len(sentence) - 1):\n",
        "        dist_args = (Y_T[ix-1], Y_T[ix+1], sentence[ix][0])\n",
        "        if (dist_args, beta) in probs_memo:\n",
        "            probs = probs_memo[(dist_args, beta)]\n",
        "        elif dist_args in probs_memo:\n",
        "            probs = probs_memo[dist_args]\n",
        "            probs = np.exp(log_gibbs_probs_normed([p*beta for p in probs]))\n",
        "            #probs_memo[(dist_args, beta)] = probs\n",
        "        else: \n",
        "            probs = log_gibbs_probs(*dist_args)\n",
        "            probs_memo[dist_args] = probs\n",
        "            probs = np.exp(log_gibbs_probs_normed([p*beta for p in probs]))\n",
        "            #probs_memo[(dist_args, beta)] = probs\n",
        "        \n",
        "        #Y_T1.append(np.random.choice(l_pos_tags, p=probs))\n",
        "        Y_T1.append(random.choices(l_pos_tags, weights=probs)[0])\n",
        "    Y_T1.append('</s>')\n",
        "    return Y_T1\n",
        "\n",
        "def gibbs_sampling(sentence):\n",
        "    # because of how this is called in the test function\n",
        "    # K can't be an argument, has to be set externally\n",
        "    # beta_update is a function set globally that has a parameter\n",
        "    # reset=True means a new sentence, =False means a new iteration\n",
        "    # beta is set globally\n",
        "    global K, beta_update\n",
        "    beta_update(reset=True)\n",
        "    Y_T = list(np.random.choice(list(pos_tags), len(sentence), replace=True))\n",
        "    Y_T = [start_token[0]] + Y_T + [end_token[0]]\n",
        "    sentence = sentence + [end_token]\n",
        "    for _ in range(K):\n",
        "        Y_T = gibbs_iteration(Y_T, sentence)\n",
        "    return list(zip([s[0] for s in sentence], list(Y_T)))[:-1]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9-l6M-yEWeJ2"
      },
      "source": [
        "### d) Test Gibbs Sampling"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q7_gXnxa1aYl",
        "outputId": "9f14823e-b9f5-44b2-c327-77f28b005e2f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 418
        }
      },
      "source": [
        "def beta_update(reset=False):\n",
        "    pass\n",
        "\n",
        "beta = 1.\n",
        "\n",
        "def test_all_k():\n",
        "    global K\n",
        "    for K in [2, 5, 10, 50, 100, 500, 1000]:\n",
        "        print('----- K = ' + str(K) + ' -----')\n",
        "        score_prediction_function(gibbs_setup, gibbs_sampling)\n",
        "        if K < 1000: print()\n",
        "\n",
        "test_all_k()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "----- K = 2 -----\n",
            "doing setup, building lookup tables...\n",
            "making predictions...\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            " 43%|████▎     | 11741/27152 [00:04<00:06, 2338.22it/s]\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "error",
          "ename": "KeyboardInterrupt",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-7-21f4ab4f30f2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     11\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mK\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m1000\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m \u001b[0mtest_all_k\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
            "\u001b[0;32m<ipython-input-7-21f4ab4f30f2>\u001b[0m in \u001b[0;36mtest_all_k\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mK\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m50\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m500\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1000\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'----- K = '\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m' -----'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m         \u001b[0mscore_prediction_function\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgibbs_setup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgibbs_sampling\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mK\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m1000\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-5-233a726878c3>\u001b[0m in \u001b[0;36mscore_prediction_function\u001b[0;34m(start_work_func, pred_func, do_setup)\u001b[0m\n\u001b[1;32m     22\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m                 \u001b[0mgold_standards\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcurrent_sent\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 24\u001b[0;31m                 \u001b[0mpred\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpred_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcurrent_sent\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     25\u001b[0m                 \u001b[0mlog_prob\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlog_sentence_probability\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpred\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m                 \u001b[0mpredictions\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpred\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-6-999c9af1abec>\u001b[0m in \u001b[0;36mgibbs_sampling\u001b[0;34m(sentence)\u001b[0m\n\u001b[1;32m     51\u001b[0m     \u001b[0msentence\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msentence\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mend_token\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     52\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 53\u001b[0;31m         \u001b[0mY_T\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgibbs_iteration\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY_T\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msentence\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     54\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ms\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msentence\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY_T\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-6-999c9af1abec>\u001b[0m in \u001b[0;36mgibbs_iteration\u001b[0;34m(Y_T, sentence)\u001b[0m\n\u001b[1;32m     31\u001b[0m             \u001b[0mprobs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlog_gibbs_probs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdist_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     32\u001b[0m             \u001b[0mprobs_memo\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdist_args\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 33\u001b[0;31m             \u001b[0mprobs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlog_gibbs_probs_normed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mbeta\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     34\u001b[0m             \u001b[0;31m#probs_memo[(dist_args, beta)] = probs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-6-999c9af1abec>\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(probs)\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;31m# computes log normalization factor to be subtracted from above\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mlog_gibbs_probs_normed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m:\u001b[0m    \u001b[0mprobs\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mfunctools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlogaddexp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0mprobs_memo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vfhIzaqDTuwK"
      },
      "source": [
        "summary of results:\n",
        "\n",
        "| K    | Runtime (s) | Tag Accuracy | Log-Prob of Preds |\n",
        "|------|-------------|--------------|-------------------|\n",
        "| 2    | 9.59        | 83.70%       | -177110           |\n",
        "| 5    | 11.44       | 86.89%       | -169121           |\n",
        "| 10   | 13.56       | 87.30%       | -168110           |\n",
        "| 50   | 43.80       | 87.53%       | -167661           |\n",
        "| 100  | 79.96       | 87.54%       | -167485           |\n",
        "| 500  | 374.08      | 87.77%       | -167251           |\n",
        "| 1000 | 740.01      | 87.55%       | -167289           |"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wmqKv1tzWpyM"
      },
      "source": [
        "### e) Setting $\\beta$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8kByuvxKW_vt"
      },
      "source": [
        "for beta in [0.5, 2., 5.]:\n",
        "    print('----- beta = ' + str(beta) ' -----')\n",
        "    test_all_k()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s24G33J0Z6JJ"
      },
      "source": [
        "| K    | $\\beta$ | Runtime (s) | Tag Accuracy | Log-Prob of Preds |\n",
        "|------|------|-------------|--------------|-------------------|\n",
        "| 2    | 0.5     |   9.83   |    70.16%    |   -219167          |\n",
        "|      |   2.0   |    9.54         |      86.23%        |      -171951             |\n",
        "|      | 5.0     |     9.51        |      86.86%        |      -170494             |\n",
        "| 5    |   0.5   |  14.05      |    77.18%    |      -195628     |\n",
        "|      |    2.0  |       10.72      |        88.19%      |       -166613            |\n",
        "|      |    5.0  |    10.68         |     88.48%         |       -166334            |\n",
        "| 10    |   0.5   |   18.21     |    77.75%    |      -193608     |\n",
        "|      |    2.0  |      13.29       |       88.52%       |     -166022              |\n",
        "|      |    5.0  |     13.18        |       88.74%       |      -166093             |\n",
        "| 50   |   0.5   |   57.85     |   78.20%     |        -193383   |\n",
        "|      |    2.0  |      43.09       |      89.02%        |      -165271             |\n",
        "|      |    5.0  |      42.30       |      89.22%        |      -165342             |\n",
        "| 100    |   0.5   |   98.01     |   78.06%     |      -193377     |\n",
        "|      |    2.0  |       79.18      |        89.01%      |        -164980           |\n",
        "|      |    5.0  |  78.82           |      89.29%        |       -165162            |\n",
        "| 500    |   0.5   |    427.28    |    78.32%    |    -192331       |\n",
        "|      |    2.0  |    369.96         |   89.28%           |      -164813             |\n",
        "|      |    5.0  |      369.26       |        89.26%      |   -164846                |\n",
        "| 1000    |   0.5   |   798.24     |   77.90%     |      -193385     |\n",
        "|      |    2.0  |      733.53       |       89.26%       |     -164669              |\n",
        "|      |    5.0  |      732.23       |      89.30%        |       -164732            |\n",
        "\n",
        "I see a trend where the tag accuracy increases sharply when moving from $\\beta=0.5$ to $\\beta=2.0$, and a small increase from $\\beta=2.0$ to $\\beta=5.0$. The trend with $K$ is the same as last time, where higher $K$ means better performance, but with diminishing returns as $K$ grows.\n",
        "\n",
        "Interestingly enough, when $\\beta=0.5$, everything runs noticably slower than when $\\beta > 1$, but I'm exponentiating in log space, so either way the same amount of multiplication is happening. I think that's going on is that lower $\\beta$ means that the sampler is more likely to pick more unlikely tags, so there is more variety of combinations of $(Y_{t-1}, Y_{t+1}, X_t)$ that show up as it runs, so there are a lot more cache misses."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PUgPH4CTdOmu"
      },
      "source": [
        "### f) Annealing $\\beta$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JpMWVM_-dSMa"
      },
      "source": [
        "# annealing schedule from the assignment\n",
        "def beta_update(reset=False):\n",
        "    global beta\n",
        "    if reset: beta = 0.1\n",
        "    else: beta += 0.1\n",
        "            \n",
        "test_all_k()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "384bNAfqdlno"
      },
      "source": [
        "#### $\\beta_0 = 0.1$, $\\beta_t = 0.1+\\beta_{t-1}$ \n",
        "\n",
        "This is a little better than constant $\\beta$ for high $K$. The tag accuracies are pretty terrible for low $K$ as expected (because $\\beta$ doesn't reach 1 until $K=10$). However, above $K=10$, it gets a whole fraction of a percent advantage. \n",
        "\n",
        "Unfortunately, especially for large $K$, this takes a very significant amount longer to compute because there are so many more possibilities to memoize, and so many fewer iterations per beta. Basically the values that *are* memoized are spread more sparsely across the whole space, because there are $K$ options for beta as well. I tried to memoize the intermediate step of the distribution before exponentiating to beta, but this didn't help much, because exponentiation and normalization are still pretty expensive :\\(.\n",
        "\n",
        "Things are speeding up a bit as it runs, but since as it continued all the cache misses became on really rare probabilities, so the speed sort of leveled out. I want more cores but I'm too lazy to set up a Jupyter server on Midway2, SSH tunnel it into my laptop, port forward it, then connect Colaboratory to it... such is life. I am going to go read Zen and the Art of Motorcycle Maintenance, hopefully it will teach me the patience to cope with staring at this glacially slow-moving TQDM bar.\n",
        "\n",
        "It has now finished running. I just spent 1700 seconds of my life staring at a progress bar for 0.01% improvement. :''(\n",
        "\n",
        "| K    | Runtime (s) | Tag Accuracy | Log-Prob of Preds |\n",
        "|------|-------------|--------------|-------------------|\n",
        "| 2    | 11.24       | 40.50%       | -339231           |\n",
        "| 5    | 23.56       | 79.19%       | -189529           |\n",
        "| 10   |   38.56     | 87.54%    | -167990           |\n",
        "| 50   |   152.74     |  89.44%     |      -164495     |\n",
        "| 100  |  219.67     |    89.63%    |    -164340       |\n",
        "| 500  |  473.38*   |   89.77%     |   -164454       |\n",
        "| 1000 |    1706.84**  |   89.78%     |     -164365      |\n",
        "\n",
        "\\* artifically low, I started this run and let it get about halfway through, then killed it to add a progress bar because I thought it might be hanging, then restarted it- so there were a lot of things that got memoized the first time around that probably drastically sped up the second run\n",
        "\n",
        "\\*\\* this column is honestly meaningless at this point, I got so fed up with $K=1000$ that I went absolutely crazy with the cacheing and micro-optimizations but was too tired to re-run the rest.... sorry :/"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dWErJApNePck"
      },
      "source": [
        "def test_all_k():\n",
        "    global K\n",
        "    for K in [1000]:\n",
        "        print('----- K = ' + str(K) + ' -----')\n",
        "        score_prediction_function(gibbs_setup, gibbs_sampling)\n",
        "        if K < 1000: print()\n",
        "\n",
        "# my own idea- exponential annealing!\n",
        "def beta_update(reset=False):\n",
        "    global beta\n",
        "    if reset: beta = 0.01\n",
        "    else: beta *= 1.01\n",
        "        \n",
        "test_all_k()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S4yC3oppe7ts"
      },
      "source": [
        "#### $\\beta_0 = 0.01$, $\\beta_t = 1.01 * \\beta_{t-1}$ \n",
        "\n",
        "I picked these constants because they were nice factors of 10 where $\\beta$ didn't end up larger than the number of atoms in the universe after 1000 iterations :P. $\\beta$ ends up around 200, which is only about double the max $\\beta$ from the last one, but it has much more spread because it starts so much lower and initially grows so much slower.\n",
        "\n",
        "The way that we're running this simulation basically is doing stochastic gradient descent in a space of tag states, with a \"metric\" of transition probabilities. This choice of $\\beta$ means the simulation has a longer period at first to get out of local \"islands\" of high-probability states, then gets more and more \"confined\" to the highest-probability states it finds as $\\beta$ grows quickly.\n",
        "\n",
        "Unsurprisingly at lower $K$ this does very poorly, because $\\beta$ stays very small the whole time, so decisions are made pretty much at random. I could probably tune constants to make it better for certain $K$.\n",
        "\n",
        "At high $K$, specifically at 1000, this has the 3rd-best performance I've seen in any implementation including MBR, well outperforming constant $\\beta$ or linear annealing. This is because high $\\beta$ means that the tags are unlikely to change between iterations because the distribution is so imbalanced, so not many iterations are needed at that high $\\beta$- the simulation will stay within whatever small island of high-probability states it settled on in the earlier rounds. At lower $\\beta$ there's likely to be a lot of movement, which is good, because it means the tag state will be able to jump across islands of local high-probability states, then eventually settle on the global area of high-probability states, and then let a high $\\beta$ do the rest of the work to make sure the final stable state is the most probable one in the local area.\n",
        "\n",
        "| K    | Runtime (s) | Tag Accuracy | Log-Prob of Preds |\n",
        "|------|-------------|--------------|-------------------|\n",
        "| 2    |  7.93   |     9.48%  |          -506580|\n",
        "| 5    |    19.32    |   9.38%     |  -506477       |\n",
        "| 10   |  32.89      | 9.51%  |      -506464   |\n",
        "| 50   |   162.79     |     9.58%  |     -505395     |\n",
        "| 100  |   347.86    |   9.93%     |     -501930      |\n",
        "| 500  |  1606.49*   |    88.65%    |    -165397      |\n",
        "| 1000 |    2899.95*  |   90.05%     |    -163752       |\n",
        "\n",
        "\\* ran out of RAM, had to turn off part of the caching, so this is slower than it would have been with more memory :/"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oi6CTqKqfNrJ"
      },
      "source": [
        "## 2. Gibbs Sampling for Minimum Bayes Risk Inference"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7Ps7HlOKfcyz"
      },
      "source": [
        "### a) MBR and MAP\n",
        "\n",
        "$\\hat{\\textbf{y}}= \\arg\\min_\\textbf{y} \\sum_{\\textbf{y'}} P(\\textbf{Y}=\\textbf{y'} | \\textbf{X} = \\textbf{x}) \\;\\text{cost}(\\textbf{y}, \\textbf{y}')$\n",
        "\n",
        "$= \\arg\\min_\\textbf{y} \\left(P(\\textbf{Y}=\\textbf{y} | \\textbf{X} = \\textbf{x}) \\;\\text{cost}(\\textbf{y}, \\textbf{y})+ \\sum_{\\textbf{y'}\\neq y} P(\\textbf{Y}=\\textbf{y'} | \\textbf{X} = \\textbf{x}) \\;\\text{cost}(\\textbf{y}, \\textbf{y}')\\right)$\n",
        "\n",
        "$= \\arg\\min_\\textbf{y} \\left(P(\\textbf{Y}=\\textbf{y} | \\textbf{X} = \\textbf{x}) \\cdot 0+ \\sum_{\\textbf{y'}\\neq y} P(\\textbf{Y}=\\textbf{y'} | \\textbf{X} = \\textbf{x}) \\cdot 1\\right)$\n",
        "\n",
        "$= \\arg\\min_\\textbf{y} \\sum_{\\textbf{y'}\\neq y} P(\\textbf{Y}=\\textbf{y'} | \\textbf{X} = \\textbf{x})$\n",
        "\n",
        "$= \\arg\\min_\\textbf{y} P(\\textbf{Y}\\neq \\textbf{y}|\\textbf{X}=\\textbf{x})$\n",
        "\n",
        "$= \\arg\\min_\\textbf{y} 1- P(\\textbf{Y}=\\textbf{y}|\\textbf{X}=\\textbf{x})$\n",
        "\n",
        "$= \\arg\\max_\\textbf{y}  P(\\textbf{Y}=\\textbf{y}|\\textbf{X}=\\textbf{x})$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hKgZnIFthaUl"
      },
      "source": [
        "### b) Approximation of Hamming cost MBR\n",
        "\n",
        "Same idea as 0-1 cost, instead of looking at the probability over all possible samples that $Y_t=y | \\textbf{X}=\\textbf{x}$, just average its indicator over a bunch of samples for that value of $\\textbf{x}$.\n",
        "\n",
        "$P(Y_t=y|\\textbf{X}=\\textbf{x}) = \\sum_{y_{-t}} P(Y_t=y, Y_{-t}=y_{-t}|\\bf{X}=\\bf{x})$\n",
        "\n",
        "$\\simeq \\frac{1}{K} \\sum_{i=1}^K \\mathbb{I}[\\tilde{y_t}^{(i)} = y]$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vWvqv6fHuOV9"
      },
      "source": [
        "### c) Implementing MBR inference"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OCu3ddQruejS"
      },
      "source": [
        "My plan here is to keep a count of how many times a tag is predicted at a certain index, then just take the maximum tag at each index."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZgW14O1OuN1v"
      },
      "source": [
        "from collections import Counter\n",
        "\n",
        "MBR_setup = lambda: None\n",
        "\n",
        "def table_score_prediction_function(start_work_func, pred_func, do_setup=True):\n",
        "    start_time = time.time()\n",
        "    if do_setup:\n",
        "        start_work_func()\n",
        "    with open('data/en_ewt.dev', 'r') as trainfile:\n",
        "        current_sent = [start_token]\n",
        "        predictions = []\n",
        "        gold_standards = []\n",
        "        total_log_prob = 0\n",
        "        lines = trainfile.readlines()\n",
        "        for row in lines:\n",
        "            row = row.split()\n",
        "            if len(row) == 0:\n",
        "                gold_standards.append(current_sent)\n",
        "                pred = pred_func(current_sent)\n",
        "                log_prob = log_sentence_probability(pred)\n",
        "                predictions.append(pred)\n",
        "                total_log_prob += log_prob\n",
        "                current_sent = [start_token]\n",
        "            else:\n",
        "                current_sent.append(row)\n",
        "\n",
        "    end_time = time.time()\n",
        "    elapsed_time = end_time - start_time\n",
        "\n",
        "    tagging_accuracy = calc_accuracy_rate(predictions, gold_standards)\n",
        "    return f' {elapsed_time} | {tagging_accuracy * 100}% |' \n",
        "def MBR_inference(sentence):\n",
        "    # because of how this is called in the test function\n",
        "    # K can't be an argument, has to be set externally\n",
        "    # beta_update is a function set globally that has a parameter\n",
        "    # reset=True means a new sentence, =False means a new iteration\n",
        "    # beta is set globally\n",
        "    global K, beta_update, beta\n",
        "    beta_update(reset=True)\n",
        "    Y_T = list(np.random.choice(list(pos_tags), len(sentence), replace=True))\n",
        "    Y_T = [start_token[0]] + Y_T + [end_token[0]]\n",
        "    sentence = sentence + [end_token]\n",
        "    counters = [Counter() for i in range(len(sentence))]\n",
        "    for _ in range(K):\n",
        "        Y_T = gibbs_iteration(Y_T, sentence)\n",
        "        for ix, y_t in enumerate(Y_T):\n",
        "            counters[ix][y_t] += 1\n",
        "    return list(zip([s[0] for s in sentence], \n",
        "                    [y_t.most_common(1)[0][0] for y_t in counters]))[:-1]\n",
        "\n",
        "def test_all_k_MBR():\n",
        "    global K, beta\n",
        "    for K in [2, 5, 10, 50, 100, 500, 1000]:\n",
        "        result_str = '|'\n",
        "        if K == 2:\n",
        "            result_str += ' ' + str(beta) + ' '\n",
        "        result_str += '| ' + str(K) + ' | '\n",
        "        result_str += table_score_prediction_function(MBR_setup, \n",
        "                                                MBR_inference)\n",
        "        print(result_str)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J7NP03MNxItS",
        "outputId": "a4f477c3-6cd1-49db-de13-90585dea8bb0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 365
        }
      },
      "source": [
        "def beta_update(reset=False):\n",
        "    pass\n",
        "\n",
        "# this just outputs the rows of the table, i was being lazy\n",
        "for beta in [0.5, 1., 2., 5.]:\n",
        "    test_all_k_MBR()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "| 0.5 | 2 |  9.330647945404053 | 55.60179728933412% |\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "error",
          "ename": "KeyboardInterrupt",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-9-b4cfa2debb55>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m# this just outputs the rows of the table, i was being lazy\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mbeta\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m5.\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m     \u001b[0mtest_all_k_MBR\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
            "\u001b[0;32m<ipython-input-8-0cba1c7a8aaf>\u001b[0m in \u001b[0;36mtest_all_k_MBR\u001b[0;34m()\u001b[0m\n\u001b[1;32m     57\u001b[0m         \u001b[0mresult_str\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m'| '\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m' | '\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     58\u001b[0m         result_str += table_score_prediction_function(MBR_setup, \n\u001b[0;32m---> 59\u001b[0;31m                                                 MBR_inference)\n\u001b[0m\u001b[1;32m     60\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_str\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-8-0cba1c7a8aaf>\u001b[0m in \u001b[0;36mtable_score_prediction_function\u001b[0;34m(start_work_func, pred_func, do_setup)\u001b[0m\n\u001b[1;32m     17\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrow\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m                 \u001b[0mgold_standards\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcurrent_sent\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m                 \u001b[0mpred\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpred_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcurrent_sent\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m                 \u001b[0mlog_prob\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlog_sentence_probability\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpred\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m                 \u001b[0mpredictions\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpred\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-8-0cba1c7a8aaf>\u001b[0m in \u001b[0;36mMBR_inference\u001b[0;34m(sentence)\u001b[0m\n\u001b[1;32m     43\u001b[0m     \u001b[0mcounters\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mCounter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msentence\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 45\u001b[0;31m         \u001b[0mY_T\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgibbs_iteration\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY_T\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msentence\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     46\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_t\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY_T\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m             \u001b[0mcounters\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mix\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0my_t\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-6-999c9af1abec>\u001b[0m in \u001b[0;36mgibbs_iteration\u001b[0;34m(Y_T, sentence)\u001b[0m\n\u001b[1;32m     31\u001b[0m             \u001b[0mprobs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlog_gibbs_probs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdist_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     32\u001b[0m             \u001b[0mprobs_memo\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdist_args\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 33\u001b[0;31m             \u001b[0mprobs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlog_gibbs_probs_normed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mbeta\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     34\u001b[0m             \u001b[0;31m#probs_memo[(dist_args, beta)] = probs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-6-999c9af1abec>\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(probs)\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;31m# computes log normalization factor to be subtracted from above\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mlog_gibbs_probs_normed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m:\u001b[0m    \u001b[0mprobs\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mfunctools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlogaddexp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0mprobs_memo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i7TTbspxxo9A"
      },
      "source": [
        "For low $K$, accuracy increases with $\\beta$. This is probably because there are so few iterations that quickly pulling the tags into a local minima (with very high-probability tags from whatever state it start in) ends up being the best on average.\n",
        "\n",
        "With high $K$, $\\beta=1$ actually has the best performance, then $0.5$, then $2$, then $5$. I think this is because with more iterations for the tag state to wander around, and since (assuming the model is right) regression to the mean means that if these states are moving around enough they'll end up doing a random-ish walk centered around the \"correct\" tagging. \n",
        "\n",
        "Lower betas will move around enough to escape any local minima and explore around the \"correct\" answer, so that their average ends up being closest to the truth. \n",
        "\n",
        "Higher $\\beta$ probably has a better chance of being stuck in a local minima during all the iterations, because it has a very high probability of picking between only a few highly-locally-likely tags each time and having an average around that local minima rather than around the correct answer.\n",
        "\n",
        "| $\\beta$ | $K$ | Runtime (s) | Tag Accuracy |\n",
        "|------|------|-------------|--------------|\n",
        "| 0.5 | 2 |  3.44 | 55.59% |\n",
        "|| 5 |  4.53 | 82.44% |\n",
        "|| 10 |  6.23 | 86.79% |\n",
        "|| 50 |  18.68 | 89.31% |\n",
        "|| 100 |  34.08 | 89.56% |\n",
        "|| 500 |  157.98 | 89.92% |\n",
        "|| 1000 |  303.68 | 89.91% |\n",
        "| 1.0 | 2 |  4.38 | 69.62% |\n",
        "|| 5 |  5.31 | 87.35% |\n",
        "|| 10 |  6.68 | 88.75% |\n",
        "|| 50 |  19.76 | 89.65% |\n",
        "|| 100 |  31.71 | 89.97% |\n",
        "|| 500 |  139.05 | 90.15% |\n",
        "|| 1000 |  272.64 | 90.18% |\n",
        "| 2.0 | 2 |  3.78 | 72.55% |\n",
        "|| 5 |  4.63 | 88.33% |\n",
        "|| 10 |  6.00 | 89.09% |\n",
        "|| 50 |  16.77 | 89.70% |\n",
        "|| 100 |  30.23 | 89.59% |\n",
        "|| 500 |  138.66 | 89.89% |\n",
        "|| 1000 |  272.35 | 89.85% |\n",
        "| 5.0 | 2 |  5.98 | 73.27% |\n",
        "|| 5 |  6.21 | 88.49% |\n",
        "|| 10 |  6.69 | 88.87% |\n",
        "|| 50 |  17.31 | 89.23% |\n",
        "|| 100 |  30.42 | 89.25% |\n",
        "|| 500 |  136.71 | 89.49% |\n",
        "|| 1000 |  268.00 | 89.53% |"
      ]
    }
  ]
}