LassiterGoodman2013.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats, special, integrate\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Full model:\n",
    "\n",
    "\\begin{align}\n",
    "P_{L_1} ( A, \\theta_{tall}, \\theta_{short} \\mid u ) & \\propto P_{S_1} ( u \\mid A, \\theta_{tall}, \\theta_{short} ) P_{L_1} (A) p(\\theta_{tall}) p(\\theta_{short}) \\\\\n",
    "P_{S_1} ( u \\mid A, \\theta_{tall}, \\theta_{short} ) & = \\frac{e^{ \\alpha (\\log ( P_{L_0} (A \\mid u, \\theta_{tall}, \\theta_{short} ) ) - C(u) ) }}{\\sum_u e^{ \\alpha ( \\log ( P_{L_0} (A \\mid u, \\theta_{tall}, \\theta_{short} ) ) - C(u) ) }}\n",
    "\\\\\n",
    "P_{L_0} ( A \\mid u, \\theta_{tall}, \\theta_{short} ) & = P_{L_0} ( A \\mid [u]^{\\theta_{tall}, \\theta_{short} } = 1 )\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simplification:\n",
    "\n",
    "$$\n",
    "e^{ \\alpha \\log ( P_{L_0} (A \\mid u, \\theta_{tall}, \\theta_{short} ) ) - \\alpha C(u) } \\\\\n",
    "= e^{ \\log ( P_{L_0} (A \\mid u, \\theta_{tall}, \\theta_{short} )^{\\alpha} ) } e^{- \\alpha C(u)} \\\\\n",
    "= P_{L_0} (A \\mid u, \\theta_{tall}, \\theta_{short} )^{\\alpha}  e^{- \\alpha C(u)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- $P_{L_0} ( A \\mid [u]^{\\theta_{tall}} = 1 )$ is the integral from $\\theta_{tall}$ to $\\infty$ of $P(A)$\n",
    "- $P_{L_0} ( A \\mid [u]^{\\theta_{short}} = 1 )$ is the integral from $-\\infty$ to $\\theta_{tall}$ of $P(A)$\n",
    "- $\\alpha = 4$\n",
    "- Alternative utterances: {$A_{pos}$, $A_{neg}$, silence}\n",
    "- $P(\\theta_A)$ is uniform in [0, 1] for any adjective A\n",
    "- $C(u) = 2/3 \\times length(u)$. This means: C(silence) = 0, C(tall, short) = 2\n",
    "- For relative standard adjectives: $P(A) = N(0.5, 0.15) $ truncated in [0,1]\n",
    "- $p(\\theta_{tall})$ and $p(\\theta_{short})$ are uniform in the unit interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pragmatic_listener(degree, theta_short, theta_tall, prior_A, alpha=4, index_signal=1):\n",
    "    \"\"\"\n",
    "    Return a tuple with the (unnormalized) probability of that combination of input parameters \n",
    "    given sent signals [short, tall].\n",
    "    \"\"\"\n",
    "    \n",
    "    cost_silence = 0\n",
    "    cost_tall = 4/3\n",
    "    cost_short = 4/3\n",
    "    \n",
    "    # probability of degree for literal listener if thetas are as input\n",
    "    p_L0_short = prior_A.pdf(degree) / prior_A.cdf(theta_short) if degree < theta_short else 0.\n",
    "    p_L0_tall = prior_A.pdf(degree) / prior_A.sf(theta_tall) if degree > theta_tall else 0.\n",
    "    p_L0_silence = prior_A.pdf(degree)\n",
    "    \n",
    "    unnorm_p_s1_short = (p_L0_short ** alpha) * np.exp(-alpha*cost_short)\n",
    "    unnorm_p_s1_tall = (p_L0_tall ** alpha) * np.exp(-alpha*cost_tall)\n",
    "    unnorm_p_s1_silence = (p_L0_silence ** alpha) * np.exp(-alpha*cost_silence)\n",
    "    \n",
    "    s1_normalization_constant = unnorm_p_s1_short + unnorm_p_s1_tall + unnorm_p_s1_silence\n",
    "    \n",
    "    # probability for listener of that combination of degree, theta_tall, and theta_short\n",
    "    # given the two possible signals \"short\" and \"tall\"\n",
    "    l_1 = np.array(\n",
    "        [unnorm_p_s1_short, unnorm_p_s1_tall, unnorm_p_s1_silence]) / s1_normalization_constant * prior_A.pdf(degree)\n",
    "\n",
    "    return l_1[index_signal]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rough and ready MCMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MCMC(prior_property, index_signal):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    index_signal: int\n",
    "        0 for \"short\" and 1 for \"tall\"\n",
    "    \"\"\"\n",
    "    \n",
    "    point = np.array([0.5, 0.5, 0.5])\n",
    "    point_prob = pragmatic_listener(*point, prior_property, alpha=4, index_signal=index_signal)\n",
    "    \n",
    "    samples = []\n",
    "    for i in range(500000):\n",
    "        \n",
    "        while True:\n",
    "            proposal = np.random.normal(loc=point, scale=0.01)\n",
    "            if np.all((proposal > 0) & (proposal < 1)):\n",
    "                break\n",
    "        \n",
    "        proposal_prob = pragmatic_listener(*proposal, prior_property, alpha=4, index_signal=index_signal)\n",
    "        \n",
    "        acceptance_ratio = proposal_prob / point_prob\n",
    "        \n",
    "        if np.random.uniform() <= acceptance_ratio:\n",
    "            point = proposal\n",
    "            point_prob = proposal_prob\n",
    "            \n",
    "        samples.append(point)\n",
    "            \n",
    "    return np.array(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 327,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw(results, adj_name, index_signal_to_plot, ax, prior_property, burn_in=1000):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    results: array\n",
    "        Array containing the results of the MCMC run from function above.\n",
    "    adj_name: str\n",
    "        String with the name of the signal. E.g. \"tall\" or \"short\"\n",
    "    index_signal_to_plot: int\n",
    "        1 for negative adjective and 2 for positive adjective\n",
    "    ax: plt ax object\n",
    "        Ax on which to plot results\n",
    "    prior_property: scipy dist object\n",
    "        As long as it has a pdf method, it's all good\n",
    "    burn_in: int\n",
    "        Number of initial MCMC samples to discard\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    None\n",
    "    \"\"\"\n",
    "    \n",
    "    sample_degree, to_plot = results.T[(0, index_signal_to_plot),:]\n",
    "\n",
    "    sns.kdeplot(\n",
    "        sample_degree[burn_in:], \n",
    "        color=\"blue\", \n",
    "        clip=(0,1), \n",
    "        label=\"$p($ degree $\\mid$ \\\"{}\\\" $)$\".format(adj_name), \n",
    "        ax=ax\n",
    "    )\n",
    "    # CAREFUL TO PLOT THE SIGNAL THAT YOU RUN THE MCMC FOR!\n",
    "    sns.kdeplot(\n",
    "        to_plot[burn_in:], \n",
    "        color=\"green\", \n",
    "        clip=(0,1), \n",
    "        label=\"$p($ transition $\\mid$ \\\"{}\\\" $)$\".format(adj_name), \n",
    "        ax=ax\n",
    "    )\n",
    "\n",
    "    xs = np.linspace(0, 1, 1000)\n",
    "    ys = prior_property.pdf(xs)\n",
    "    ax.plot(xs, ys, color=\"red\", label=\"$p($ degree $)$\")\n",
    "\n",
    "    legend = ax.legend(fontsize=\"small\")\n",
    "    legend_frame = legend.get_frame()\n",
    "    legend_frame.set_facecolor(\"white\")\n",
    "    ax.set_xlabel(\"Degrees\")\n",
    "    ax.set_ylabel(\"Density\")\n",
    "\n",
    "    for text in [\"top\", \"left\", \"right\"]:\n",
    "        ax.spines[text].set_visible(False)\n",
    "    ax.set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure 1. & left part of Figure 2 (for thesis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 301,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### RELATIVE ADJECTIVE'S DISTRIBUTION (FIG. 1)\n",
    "myclip_a = 0\n",
    "myclip_b = 1\n",
    "my_mean = 0.5\n",
    "my_std = 0.15\n",
    "a, b = (myclip_a - my_mean) / my_std, (myclip_b - my_mean) / my_std\n",
    "prior_property_tall = stats.truncnorm(a=a, b=b, loc=my_mean, scale=my_std)\n",
    "\n",
    "results_relative_tall = MCMC(prior_property, index_signal=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### HALF-OPEN SCALE ADJECTIVE'S DISTRIBUTION (FIG. 2)\n",
    "prior_property_dangerous = stats.beta(a=1, b=9)\n",
    "\n",
    "results_absolute_dangerous = MCMC(prior_property, index_signal=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 331,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAADnCAYAAABSbO4uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd3hUZfbHv28y6Qmkk5BGSQKhmNCrUkUUKS4qWFYRy8qiKOradsWIjVUsiCsi/AC7olIUUZQuRQQCISEESEiDBNIrqTPn98ebCZlMuzOZOyV5P88zzyT3vve975175sy55z3nvIyIIBAIBAL5cbL1AAQCgaCzIBSuQCAQWAmhcAUCgcBKCIUrEAgEVkIoXIFAILASQuEKBAKBlVDI1K+INRO0F2aDcwq5FbQXg3IrLFyBQCCwEkLhCgQCgZUQClcgEAishFw+XItRXl6O4uJiNDY22nooAplwcXFBYGAgfH19LdpvUlLSLIVC8TwRhUhp7+rqatHzCzo25sit3SvcgoIC9OjRA+7u7mDMFvMoAjkhItTV1SE7O9viCtfZ2fmd3r1713p6ehZLlJ0oiw5A0GExV27tXuECgIeHh62HIJAJxpic99fZy8urTq7OBZ0Xc+VW+HAFAoHASgiFKxAIBFZCKFyBQCCwEkLhSuTnn3/GZ599pnd/eno65s2bZ70BtaKqqgpFRUV69+/duxeJiYkYP348ACA7Oxu//fabwfbPPPNMy7uuPtTk5eVh8eLF7b6GjkxHkp3W2HLcUrE3uRUKVyJr167F3LlzZelbpVK16/jjx4/jhx9+kNzemMI1hYiICFy+fBnl5eUW6a8j0pFkxxa09xp1YSu5dYgoBTVPPgmcPNm+PhISgPffN9wmPj4ekZGRyM3Nxfz583H//ffj6tWrWnGaTU1NmDt3LsrLyxEXFweAh4ssWrQIqampUCgUWL9+PUJCQjTaVVVVYd68eVi+fDmcnJwwe/Zs3HfffVrHhYWFaW0LDw83+ZqHDx+OAQMGYP78+QCAVatW4dChQzh27Bg+/fRTPPjgg2hoaEBwcDC+/fZbSX205vrrr8eOHTswZ84ck8dmDebPR0RqKjyNtfM00kLIju5xX7lyBffcc4+GDP3xxx9YtmwZ3NzckJ2djS+++AJxcXFa17Jhwwad152RkaFxjTt27MClS5fg5eWFL774AqdOncK2bduwfPlypKenY9myZfjHP/6BJ554Ah4eHhg3bhyWLl1qd3IrLNw2lJeXo7S0FBs2bMCRI0fw6aef4ty5c4iMjNRqu2XLFsTGxmLnzp247rrrAPDHRz8/P+zZswfLli3DsmXLsHnz5pZ28fHxLcdXVlZi69atuP/++3Uep2ubOXh6eiIwMLDlGhYsWIA5c+Zg7969CA0Nxa+//or9+/cjMjISu3fvltRHa3r16oW0tDSzxtaR6Ayyo2vcfn5+OmWosbERW7duxfLly7F+/Xq916JvrOpr9PT0RGRkJPbt24c5c+Zg5cqVOse6fft2LFmyBPv27UNiYqLO8bfGFnLrUBauMevCEqSkpOCee+5BQEAAVCoVAgICAADu7u5abTMyMjBkyBAA3BI4ePAg0tLSsHnzZuzfvx9EhIiICGRmZra0GzZsGA4cOAAAGDp0aEsyh67jdG2zNKWlpXj00UdRVlaGgoICxMfHIyrKtPh/e1+IdN065ElsOqQ95+kMsqNr3PpkKCEhAQB/fC8rK9N7LfrGqr7GzMxMDBs2DAAwYsQI7Ny5U8Mnq5a/hQsX4s0338TGjRsxd+5c3HLLLQavxRZyKyzcNqSkpLQI8qeffoqpU6ciNjYWWVlZWm2jo6Nx4sQJAMCxY8cAAH379sWdd96JvXv3Yt++fVi/fj169+7d0u748eMtxzs5Xfv4dR2na5slcHFxgVKpBAB8+eWXmDJlCvbt24dbb73VLCHMyspqebzszHQG2dE1bn0y1Dq7j4j0Xou+saqvMTo6GkePHgUAHDlyBDExMfDz80NeXp5GX127dsWKFSuwdu1aPPfcc0avxRZyKxRuG1JSUlBTU4OJEydiz549WLRoEXx9feHk5IS6Os2kpVmzZiE9PR2TJk1CUlISAGD69OkoKSnBhAkTMGHCBHz22WeYNWsW0tLSMGnSJBw5cgQuLi5a59V1nK5tlmDgwIE4fvw4br/9dgwfPhyrVq3CzJkzcfnyZbP6279/P2666SaLjM2R6Qyyo2vckyZNkiRD+q7F2FhnzZqF3Nxc3HDDDfjmm2/w2GOPYeDAgairq8PkyZNbLOXVq1fjhhtuwKhRoyRFT9hEbolIjpfFSEtLs2R3RrnppptIpVJpbd+2bRutX7/e7H4bGhqIiGjNmjW0bNkys/vRxZ49e2jVqlUW7VMqubm59MQTT7S7Hx33ud1yePLkyWwiOmbCyzBlZUQFBXp3C9kxjpzXYgq2kluH8uFag8bGRp1FcqZNm9aufmfOnInq6mq4ubnpjQRwRCIiIvC+NZzr9kBFBVBWBoToLj4mZMc49nIttpJbRvI4ji3W6ZkzZ4R/sBOg4z63uzRccnJydnx8fLEJhxieNMvP56/BgwEn4Y0TmC63QmoEAqmoY2lFbWaBmQiFKxBIRT1hJRSuwEyEwhUIpKJWuA0Nth2HwGERClcgkIqwcAXtRChcgUAqCgXAmFC4ArMRClcgkApj3MoVCldgJkLhCgSm4OIifLgCsxEKVyK6ikhbsq6smpMnT2LVqlVafau3txdDBafNLVTe+m97K/hscVxdTbZwO6PstMaRC5UDlpVdoXAloquItL4vTXsKJickJGDBggVafau3txdTCk6LQuU6MMOl0BllxxbIUagcsKzsOlRq75O/PomTl9tXgTwhJAHvTzWc0ie1iHTrQt6LFy/GmjVrWgomf/755xoFmZ2dnbF3716toszV1dUaRZMnTpyIbdu2IS8vr6XvrVu34sSJE9i2bRuWLVuG++67T6sYc9t+Bw4caPJn01ELlc/fOj8itTDVeAFyF8NNEkIS8H78c4BSyV/OzlpthOy0v1D5wIEDtY6vqqrC+vXrTS5U7u/vj7179xotVv7888/rlVvAcrIrLNw2mFJEunUh765du7YUTL7rrrv0FvVuW5RZV9FkXX2r2bx5s85izG37NQdRqFwCBrLNhOxYplC5eqxti5VbqlA5oF2s3JDcApaTXYeycI1ZppbAlCLSbVEXTNZVkFlN26LMzz77rEbRZE8j67zoKsY8btw4rX4tjSMXKl83c53lCpBXVvL3xkagjUwI2blGewqVq8fatli5OYXKAe26vIDtipULC7cNphSRbl3IG7hWMNlQUe+2N19f0eS2favRVYxZV7+WQBQq14GBbDMhO5rnMrdQOQCdxcrNKVQOwK6KlQuF2wZTiki3LuRdVVXVsl1qQWZAf9FkfX3rKsYsF6JQuQ4MZJsJ2blGewqVq49vW6zcnELlAOyrWLmxgrlmvixGRylAbi/Yqli5sYLPDlGAnIhIpSI6fpwoN1drl5Ady+IIxcpFAfJ2IlcR6c5OhylUbiDbTMiOZemIxcpFAXKBXeAQBcjVpKfz9759TR+UoEMhCpALBHJjRraZQAAIhSsQmI7apWCjMDeB4yIUrkBgKq6ugErFs80EAhMQClcgMBV1tpmoGiYwEaFwBQJTEQpXYCZC4QoEpiLWNhOYiVC4AoGpuLjweFyhcAUmIhSuRHQVkW6NLYosd5ii3o6GOvlBosIVsiNQIxSuRHQVkbYU5hZO7jBFvR0RV1fJClfIjkCNSO3VgdQi0rqKLBORVpHkkJAQrWLK8+bN0yicfN9992kdFxYWprUtPDxcYwy2Kurd6XF1BaqrtTYL2REYwrEU7pNPAifbt+IDEhIAA3nR6iLSu3fvhpeXF0aPHo1Ro0bpLEysLrL8xhtvYM2aNTh48KBGkeTjx49j2bJlGDduXEu7tWvXtlQsqqysxL59+8AYw7Zt27SOmzp1qta2Dz/8UGMMvXr1wuHDh9v3mXR05s+PQKrxFR9gpJ6shuyos82IuIsBQnYExnEshWsFTCkiravIsq4iybqKKQPXCicD0HmcvoLLrZGpFobAGK6uXNk2NraEiQnZERjDsRSuFapNmVJEWl1kefbs2S1FltVFkl966SUAvILUli1bWtqpCyAD1won6zvul19+0drWlg5R1Ftu1llwxQc1rWNxWylcITsCQ4hJszaYUkRaV5FlXUWSdRVTbouu44wVXAY6SFFvR0RH8oOQHYFRjBXMNfNlMTpKAXI5iikbK+rtSDhMAXI1jY1ER48SXb7csknITudDFCBvJ3IVkZajmHKHKertiDg7A05OGhaukB2BMYTCbcOuXbtk6Xf79u2y9CuwEYxpxeIK2REYQ/hwBQJzMSH5QSAAhMIVCMxHKFyBiTiEwiURL9ihkfH+qlQqVbvXRtOLOvnBzPRagWNjjtzavcJ1d3dHSUmJULodFCJCSUmJzuQAC5BaVFTUVTalK+ridlrMlVu7X7W3sbERFy9e1IpjFHQc3N3dER4e3jbGtN1K8vjx48EKhWItgAGQYFy4urpGmXSCujrgyhWgWzdAnh8MgR1jjtzavcIVdFrkcwXoxzS5zckBevQAPvkEePhheUYkcDTEMukCgSyEhwMKBaAjdVcg0IVQuAKBuTg7A5GRQHa2rUcicBCEwnUwvk75Gv/e9W8xiWgv9OghLFyBZESmmQORXpyOB7Y+gHplPWIDYnF/wv22HpKgZ09g2zZbj0LgIAgL14F4ee/L8HDxwPCw4Xjm92fQqNQuuSewMj178kiFq1dtPRKBAyAUroOgIhV2XdiFmX1m4oWxL6D4ajH2ZO+x9bAEPXvyd+HHFUhAKFwHIbUwFSW1JZjQYwKmRk+Fj6sPNp7eaOthCXr04O/CjyuQgFC4DsLe7L0AgPE9xsNd4Y6ZfWdi05lNwq1ga4SFKzABoXAdhP05+9HDtweifHky1J397kRZXRl2ZclTElAgkZAQnmUmLFyBBITCdRCSryRjSOi1Jbem9J6CLm5dhFvB1jAmQsMEkhEK1wGoaahBZmkmBgYPbNnmpnDDzD4zsSV9C1QkqlXZFKFwBRIRCtcBSCtKA4EwsNtAje3je4xHWV0Zzpect9HIBAC4H1coXIEEhMJ1AFIKUwBAw8IFgGHdhwEAjuYftfqYBK3o2RMoL+cvgcAAQuE6AClXUuDp4olefr00tscFxcHTxRNHLwmFa1N69+bvmZkmH9qgbEDy5WQRbdJJEArXAUgtSkW/oH5wdnLW2K5wUmBw6GBh4dqamBj+ft401w4RYdpX05CwOgETP5uIq40iW62jIxSuA3C2+Cz6BvbVuW9o6FCcuHwCSpXSyqMStKC2cE1UuPtz9mPnhZ2I7xaPA7kHsP7EehkGJ7AnJClcxtgPjLFpjDGhoK1MTUMN8irz0Cegj879A4IHoK6pDtnl2dYdmANgNbn19OS1cU1UuN+lfQcvFy8cfvAwhoQOwerjq0UVuA6OVEFcBeBuAOcZY8sYY7rNLYHFOV/Kv8T6FG5cUBwA4EzxGauNyYGwntxGR5uscPfl7MPoiNHwcPHA/fH3I6UwBVnlItqhIyNJ4RLRTiK6B8BgANkAfmeMHWKMPcAYczF8tKA9nCs5BwDoE6hb4apdDenF6VYbk6NgVbmNiTFJ4ZZcLUFqYSrGRY0DAEzuNRkAsDtrt0WHJbAvJD9qMcYCAMwD8BCAEwBWgAvy77KMTACA+28BINo/Wud+fw9/BHsF40yRsHB1YTW5jYkBSkqAsjJJzf+69BcAYEzkGAD8hzPEO0Qo3A6OpALkjLFNAPoC+BzAdCIqaN71LWPsmFyDEwBnS84ismskPF089baJC4xDeomwcNtiVblVRypkZADDhhltnnwlGQCQEJKgHivGRo5tUcSCjolUC3ctEfUjojfVQssYcwMAIhoq2+gEOF96HjH+MQbb9A3sizNFZ0BEqKkBVCLTV4315NbE0LCTl08iqmsUfN19W7YldEtAZlkmquqrLDo0gf0gVeG+pmPbYUsORKCb7PJsrYSHtsQFxqGsrgxrvylESAiQkABs3w6ICW8rym3v3ryQjUSFm3wlucW6VaP+/9SVUxYfnsA+MKhwGWMhjLEhADwYY4MYY4ObX+MB6H/GFViEmoYaFNYUoodvD4Pt1JEKC15KR9++fLWXadOAp5/unErXJnLr7g5EREhSuDUNNThbfBbx3eI1tseH8P9PXj4pyxAFtseYD/cm8AmHcADvttpeBeBFmcYkaCanIgcA0NO3p8F26kgFlf8ZfLdxHLp3B556CnjvPWDAAGD+fNmHam/YRm4lhoalFqaCQFoWbphPGAI8AoTC7cAYVLhE9CmATxljs4noByuNSdBMVhmPyezpZ1jhBruHAw1eiBmV3rLiywcfACkpwLPPArNmAf7+Mg/WjrCZ3MbEABuN1ydWT5ipLVo1jDEkhCTg5BWhcDsqBhUuY+xeIvoCQA/G2FNt9xPRuzoOE1gIdRC8MZfCvr1OQHEfePW9Fhrm5AR8+CH35775JvD223KO1L6wmdzGxPCwsNJSg79wJy+fRBe3Ljrva0JIAj7860M0qZqgcJIURCRwIIxNmnk1v3sD8NHxEshIdnk23BXu6ObVzWC7zZsB5/I4FJFmLO7AgcDddwMffQQUFso5UrvDNnIbG8vf0w2H6KUXp6NfUD846cg4TghJQL2yviX+WtCxMOZSWN38/op1hiNoTVZ5Fnr49gBjzGC77duBPtPikFb1JWoaauDl6tWy7z//Ab74Ali3Dnj+eblHbB/YTG779ePvZ84Ao0frbXah7ALGRo7V3UUQ7yO9OB39g/tbfIgC2yK1eM1bjLEujDEXxtguxlgxY+xeuQfX2ckuzzY6YVZYCOTlAYOjeBxoRmmGxv4+fYBRo4BvvpFtmHaL1eW2Rw8erXBGf9Zfg7IBeZV5ekP91BmFbe+joGMgNQ53ChFVArgVwEUAsQD+JduoBAD4pJkxhXviBH+/vj9XuOraC62ZOxdITjaoBzoq1pVbZ2f+C5eWprdJbkUuVKTSq3C7uHVBsFewULgdFKkKV13o4xYAXxNRqUzjETRTUVeBsroyoxNmaoV78whuGamri7Xm9tv5+5YtlhyhQ2B9ue3Xz+Av24WyCwBgMJkl2j8aGWVC4XZEpCrcnxhj6QCGAtjFGAsCUCffsATq+rbGQsKSkviSWhHBPgjxDtGpcLt3B4YMAbZtk2Okdo315TYuDsjOBmpqdO5WK1xDTy4x/jFiYdAOitTyjM8DGAVgKBE1AqgBMFPOgXV2pIaEnTgBDB7M/zb0Rb31VuDwYaC42JKjtG9sIrfqiTM9kQqXKi/BiTkh1CdUbxfR/tG4VHVJLLnTATGlEn4cgDmMsfsA3A5gijxDEgCtLFwDllBFBS9OpVa4sQGxOi1cgCtcIuCXXyw9UrvHunLbOlJBB/lV+ejm1c1gjK164kxtDQs6DlKjFD4HsBzAWADDml+iSpiMZJVlwcfVB/4eBgLomxOSBg3i7zH+MSisKURlfaVW28GDgZCQzuVWsIncRkcDCoXeibP86nx09+luuAsRqdBhkZrKMhRAPxILLlmN7IpsozG46gmzFpdCAI9UOF9yHkO6D9Fo6+TEC9p8/z3Q2Ai4dI51Oqwvty4uPOPMgIUb1TXKYBdqhSv8uB0PqS6FVAAhcg5EoElWWZakCbPQUKBbcyKaum6uIbdCRQVw4IBFh2rP2EZu4+L0W7hV+QjzCTN4uK+7LwI9A4WF2wGRauEGAkhjjP0FoF69kYhmyDKqTg4RIas8CxN6TDDYrvWEGQD09ufLdeuzjCZPBlxdgZ9+AiYY7rqjYBu57dcP2LoVqK8H3NxaNtc31aP4arFRlwLArVx9P5wCx0Wqwk2UcxACTUprS1HdUG3Qwr16lRtRt912bZuniyfCu4Tr/aJ6ewOTJgE//gi88w6vl93BSbTJWfv1A5RK4Nw5XtCimcvVlwFAssLdn7NftiEKbIPUsLB94KueujT/fRRAkozj6tRICQlLSeFL6agnzNTE+McYtIymTwcyM43WV+kQ2ExuBwzg76mpGpvzq/IBSFS4ftHIq8hDXZMId+9ISI1SeBjA9wBWN28KA9D58pashJSQsLYTZmqMBc1Pn87ff/qpPSN0DGwmt3368Mmz5GSNzSYpXP9oEKilJrKgYyB10mwhgDEAKgGAiM4DCJZrUJ0d9ZfMkIWblAT4+QGRkZrbYwJiUFJbgtJa3Vms4eHcKv7xR0uN1q6xjdy6unK3QjsUrjriREycdSykKtx6ImpQ/8MYUwAQIWIykV2eDT93P3R176q3jXrCrK0fVr3cTnqxfp/BjBk866yoyCLDtWdsJ7fx8ToVrouTCwI8A4weLmJxOyZSFe4+xtiL4Ivy3QjgOwCd4KHUNmSVGw4Ja2wETp3SdicAfAVfADhTpL+AyvTp3P+7Y0e7h2rv2E5u4+OBggKNX7X86nyE+oTqLDzeFn8Pf/i5+wmF28GQqnCfB1AEIAXAPwBsB/AfuQbV2VEXHtfHmTNAQ4P2hBnA3RDuCnekFekvEThoEI/d3b7dAoO1b2wnt9ddx99PXVvyPL/KeJZZa0RoWMdDUlgYEakYY1sAbCGijv8gakOICNnl2ZgWM01vm6TmeXZdFq6zkzP6BPTBmWL9Fq6TE3DzzTxUVKnkZVw7IjaV2/jmBSKTk3ksHrjCVbt8pBDtH40/L/4px+gENsKghcs4iYyxYgDpAM4yxooYY0usM7zOx5WaK6hrqjMYoZCUBHh58QxSXfQL6mfQwgW4wi0rA44cac9o7RO7kNugIJ4G2MqPm1+Vj+7eplm4ORU5aFA2GG8scAiMuRSeBJ/lHUZEAUTkD2AEgDGMscWyj64Tog4JM+RSOHGCr8brpOfuxQXGIaciBzUNumuyAsCNN3LLtoO6FexDbltNnF1tvIryunKTXAr+iIaKVJj9YBZuuw3YsIFXfBM4LsYU7n0A7iKilmBAIroA4N7mfQILow4J0zdpplLxKmG63AlqWi9EqA8/P77OYQdVuPYht9ddx9MBGxtRUFUAQFpIGADs2wc89zB/hDmSkYGUFOCBB4Bnn5VttAIrYEzhuhCRVsnqZn9Y56g3ZWXUWWb6KkplZADV1bonzNTEBTVHKhjw4wLALbdwazk/37yx2jH2IbcJCTykJC3NpBjc4mKesh3lzRXuC2+fw7lzwIIFwPLlvOKbwDExpnANOY+EY0kGssuzEewVrLHUeWsMTZipifaPhjNzNurHveUW/v7rr+aM1K6xD7kd0lwi89gxkxTu0qW8qtumLwLg5+6H86Xn4eQErFgBDB0K/POfQEmJnAMXyIUxhRvPGKvU8aoCMNDIsQIzyCo3vFJvUtK1RCZ9uDq7IiYgxqiFO3AgEBbWId0K9iG30dFA164mKdzSUmDtWmDePGDAAIaYgJiWlZhdXID/+z9uAb/9ttyDF8iBQYVLRM5E1EXHy4eIhEtBBrLLs41OmA0YYLyAeFxgnFELlzFu5f7+O3/y7SjYjdw6OXEr9+hR5Fflw13hDl93X4OHrFsH1NYCTz7J/48NiG1RuAB3C991F7ByJVBYKOfgBXJgyppmgvZQWQmsXs2/LVOmcBPmm294zdRmlColcspz9Fq4RHzCzJD/Vk2/oH7ILM1EfVO9wXZTpvChqYvhCCzMsGHAqVMoKruIUO9Qgyt4AMBXXwEjR16r6hjrH4u8yjzUNta2tFmyBKirA956S86BC+RAKFy5IeJmS48ewKOPAgcPcgfd9u1c+cbFtThRC6oL0Khq1Gvh5ufzx8mEBOOnjQuMg5KURjOVRo3i74cPm3BNAukMHQo0NsIj7bxRd0JODv/h+9vfrm3TVcSmTx/gnnuAVav4j6XAcRAKV06amoD584EHH+TPgkeO8G/VkSM8z377dsDDgz/Xv/kmskr5Kq36QsLUi0ZKUbjq0DBDNRUA7sONiBAKVzaGDQMAdEvPM6pw1SUzWxeVjw2IBaC9bNLChbwI/TffWG6oAvkRClculEpuhmzYALz0ErBrFzB8+LXyXs7OPN3r6FFg7lzgxRfRZdl7APTXwVUrXHWaviH6BPYBAzPqxwW4lfunyCCVh8hIIDAQvTPLEOodarDp/v1AVBSfa1OjXqeutR8X4KLUvz+fRBM4DkLhysXzzwMbN3JH29Kl+gsWeHoCX3wBzJuH+NVbMD8JiOwaqbPpyZNA795Aly7GT+/p4okevj2MRioA3GeYk8ONboGFYQxNQwYjPq/RoIVLxBf3HDNGc7uPmw9CvEO0FC5jwEMPAX/9pbWwhMCOEQpXDr78kkeoL1wI/Otfxts7OQFr1iA1vjs++hlwS0rW2ezkSWnuBDVxQcYjFQDhx5WbyoEx6F8ERCj89bbJzuY/eG0VLsDdCrp88ffey0ME162z4GAFsiIUrqXJyeEpQddfD7z3nvTjFAq8+HBPlPq6cVdEjWYdhKoqnmVmisLtF9gP50rOoUnVZLDdoEH8iyvcCvJwqX8knAnok1mht426iJD6x681Mf4xWhYuAAQGcq/Uxo085Vtg/wiFa0lUKp7wTgR8+qnxYNk2nGq8iHVPXM816wsvaO5rLqtqqoVbr6w3ui6WmxvPXBMWrjyci+GWbURqrt42yclcXPr3194XGxCLwppCVNRpK+w5c4BLl3jwi8D+EQrXkqxfD+zZwy3bnvqzxXRR31SPvMo8NF4/Bli0iEe2HzrUst+UCAU1A4L56rEphSlG244aBRw71rESIOyFXKcqpAYBfsf1u3eSk3mEoKur9j59kQoAX73DwwP49luLDVcgI0LhWorycm6Vjh7Nw8BMJLs8GypS8bWs3niDx2stWtTyrHjyJBAQwDdLZUDwADgxJyRf1u0Tbs2oUTyYPtl4U4GJ5Ffl488oZyj+Oqb32T85+VrN8rboi1QAAG9vYNo04LvveBSiwL4RCtdSvPIKz0pYuVJ7ZUcJqAPbe/v15tXF334bOH6cW83gVf4GDDCta08XT8QGxOLklZNG244cyd+FW8HyFFQX4EysH1hFBXD6tNb+4mKe1KJP4fb270loMRsAACAASURBVA0GhvMlupNY5szhab779lly1AI5EArXEpw/zxXtI48YLuNlALXCVa/WirlzgbFjgRdeAFVV48wZ/shpKvHd4iVZuBER3HoWE2eWJ78qH3kDm0P9dDhb1U8V+hSuu8IdUb5ROFeqbeECPG/Gy0u4FRwBoXAtQWIin3l65RWzu8gozUAXty4I9AzkGxjjoWVFRah+8wOUlZmvcHMqclBeV2607ciRwsKVg4LqAlCvnnzlTjMULqA/UgHgodwzZgA//CB88PaOULjtJTUV+Ppr7m/t1s3sbjLKMhDtH61Z3GTECGD6dLivfBtdUY6+0tcfbCEhhM+ynbxs3K0wahSQlQVcuWL6eQT64av1hvEnlgMHtPYnJ/Plz4KC9PcRGxCL8yXnQXrW2Jk7l5d23LXLUqMWyIFQuO0lMRHw8ZGW4GCAzNLMa+6E1rz6Klyqy/E03jHLwh0WxnP5/7r0l9G26hhQ4VawHDUNNaisr+RZZmPH8gyHvDyNNoYmzNTE+Megor4CRVd1Lz5800289K5wK9g3QuG2h6Qk/hz31FOAv/4sImM0qZqQVZ7FJ8zaEh+PEzF3YjHeQ7iH6WX+Az0D0duvt6TltgcP5rGgwq1gOQqqeb50qE8oMHEi37h7d8v+hgY+IWpM4aqXV9eXOejmBsyaBWzerFHxU2BnCIXbHl5/HfD1vVYt2kxyK3LRpGrSbeECWBW0BN6oAftwpVn9jwwfiT8v/qn3cVSNuzvPOhMK13JorPQwYABPD2v13J+ezv2uxhRufAhvYGgCdM4cXvnzt9/aP26BPAiFay5nz3JzYuFC/izXDrQiFNrwa15/HA+fwSMhqqtN7n9E2AgUVBfgYuVFo21HjeIFzISVZBnUq/WGeofymhkTJ3ILt/nHT8qEGQCEeIcg2CsYyVf0K9zJk/mDlijZaL8IhWsuy5fz57hFi9rdlSGFW13NXX6npr3AZ0XWrDG5/xHhIwAARy4dMdp20iS+xItIFbUMWmuZTZrEc3HP8YiD5GQuRrGxxvu6rtt1BhWuiwsvXv7jj/weCuwPoXDNoaAA+OwzXjchOLjd3WWUZsBD4aGzXurZs/y965SRwPjxwDvvmGx+JoQkwM3ZTZIfd8IEnl76yy8mnUKgh4LqArg5u11by0ztx212KyQnc0+DQmG8r/hu8ThdeNpgMaI5c/iPdAdcGLRDIBSuOaxYwfMon3nGIt1llGbwbCIdaWTp6fy9b1/wGruXLvHyjybg6uyKQaGDJFm43t680JlQuJaBh4R1v3Zve/fmRcl37waRtAgFNfHd4lGvrMfZ4rN624wfz20AEa1gnwiFayoVFXwxqTvuAHr1skiXmWV6QsIAnDnDa5dHR4Ov+DhoEE/7NbEe38iwkTiefxyNSuOR8TffzDNQ20QvCcygoLpAs/A4Y9ytsGcPLl9SoqjIBIWrnjgz4FZQKIDbbwe2bTPL3S+QGaFwTWX1ar5y33PPWaQ7pUrJY3D9dCvclBQgJqa5ihRjwNNPc7N3xw6TzjMifARqm2olVQ67+Wb+3ry2paAd5Ffl85Cw1tx0E1BaipyN/IlDqsLtG9gXLk4uOHXllMF2c+ZwH656jTSB/SAUrik0NHB3wuTJ0tYql0BmWSbqlfXoH6yjECr4Kq4ap7rjDqB7d+Ddd006z8hwXp3myEXjboW4OF5bQbgV2gcR4WLlRYT5tCnxdtNNgLMzVNt+BiBtjTqAu4b6BfUzaOECPL8iLIwnQArsC6FwTWHjRl7W6emnLdZlaiFfkKp/kLbCLSnhj/UaCtfVFXj8cWDnTm7+SiSqaxSCvYLx5yXjE2eMcSt3507+GyMwj7K6MlQ3VCOqa5TmDl9fYOxYhJ3YhshIwM9Pep/xIcaLETk5AXfdxX8wS0zPlRHIiFC4UiHihcXj4riFYiFOF/JyfXFB2nm7eouOP/IIr1hiwhI+jDGMCBshycIFuMKtqhJJEO0ht4Kv8BDlG6W9c9o0RJWfwqRY0xzl8d3iUVBdgKIa3Sm+au6+m8/r/vCDSd0LZEYoXKns389TeRcvNqverT5OF51GT9+e8Hb11tqnV+H6+wPz5vFoBRMqzYwMH4mzJWdRVltmtO3EiXwCRrgVzCenPAeA7lWY62+8FQAw08W0+K3runH/gzG3QkICj2z56iuTuhfIjFC4UnnvPZ6Wee+9Fu02tTBVr//2zz+5L1VnFaknnuDP+x99JPlcI8J4AoSUQjZdunBfoIjnNB+1hatL4aY29cUF9MSwwm0m9amu/pZUkGSwHWPcyt23j9fLEdgHQuFKISODp+8sWMAXkLIQ9U31OFdyTqf/lgjYu5fHVeokNpYvaLVqleS0omFhw8DAJMXjArz7lBR++QLTyanIgbvCHUGe2r+YyacYfsJ0hKTuNCl+K9AzED19e0r60Zw3j4cUmvCbLJAZoXClsGIFz5v85z8t2u3potNoVDVicKj2KhFpaXzplQkTDHSweDFQVCQ5EaKLWxf0C+onKeMM4PGcAF8vS2A6uRW5iOwaqTOhJTkZ2O4+G071dTxo1gRGhI+QpHAjIniq75o1QE2NSacQyIRQuMYoKwPWrePTviEhFu36eP5xAMCQ0CFa+/bs4e8GFe748dxZ9957LcVQjDEyfCSOXDpitHIYwBOiRo4UCtdccipytCMUmklOBqrjx/DK4yZ+wCPCRiCvMq+lToMhFi3i65uamJwokAmhcI2xZg1w9Sq3Ji1MUkESurp1RS8/7Yy1337jiWw9ehjogDE+rrQ0yTX5RoSNQGltaUvBHGPccQePBc7MlNRc0Aq1hdsWdUrvwARnYPZs7ig3wa1wQ9QNAIBdF4wv7zBmDA8r/OADyb/JAhkRCtcQjY1cUidOlJ4OZALHC45jcOhgrUfOujpe2+SWWyR0Mncut5IkJkKoEyCEW0Fe6prqcLn6sk6Fm5fHrc74ePBftDrT3AoJIQkI8gzCjkzj2YaMcSv39OlrT00C2yEUriG+/54Xi3nqKYt33aBswKkrp3T6b/fv50a1OsXWIK6uwGOPcQv3lOGUTwDoF9QPPq4+OHxRWoBtZCRfWm3jRknNBc2oaw/rcilo1MAd0+xWMOEDdmJOmNJ7CnZk7oBSpTTafu5cHmCzapXkUwhkQihcfRBxq7FPH4mazzSO5R9DvbIeYyLGaO3bvp2vvqA3QqEtCxbwdbKXLzfa1NnJGSPDR+JQ3iHJYxVuBdMxFBKmVrgDB4KHEdxxB/Dzz7zesURm9JmB4qvFkn443d35FMRPP/HaSwLbIRSuPg4cAI4d48vnOFn+Y9qXvQ8AMDZyrNa+X37hk2WenhI78/MDHnqIJ89LKPE1OmI0UgpTUFlfKal74VYwnezybAC6s8ySk3mVRh+f5g3z5vGYahOKH0yNngpXZ1dsSd8iqf3dd/Myyps3Sz6FQAaEwtXHO+/wjK777pOl+305+9AvqB+CvDRjNDMy+GIAkvy3rVm8mFvlK1YYbTo6YjRUpJIUWgQAUVHcrSAUrnTOl5yHwkmh18LVmBIYNIhv2LBBcv9d3LpgUs9J2Jy+WVLEyYgRfBJWRCvYFqFwdZGWBmzdyovESDYzpdOobMTBvIMYFzVOa586ldZkhRsVBdx5J/DJJ0afG0eEjQADM9mtkJQEXLhg4rg6KedLz6OXXy8onDSXcqip4T+qWnOw8+bxJ6rUVMnnmNV3Fi6UXWgpgGQIdebZ7t18wRKBbRAKVxdvvcUV7WOPydL93uy9qG6oxo29btTat307TyIzq7b5v/7FK86sXm2wWVf3rhgQPMAkhSvcCqZxvvQ8YgO0Fyo7dYo/iGgp3Hvu4cUr1q+XfI4ZfWaAgZnkVlCpxASoLREKty25ufy56+GH+dSuDHyf9j28Xb0xNXqqxvarV3nojsnWrZpBg/hqAitWGK2rODpiNP68+CdUJG3liKgoYPhwoXCloCIVzpecR4x/jNa+Eyf4++C2wSlBQcDMmcCnn0pO1Q7xDsGoiFHYclaawo2L44peLL9jO4TCbYs6nlWGUDAAaFI1YXP6Ztwaeys8XDTrMuzZwyc2zFa4ALdy8/ONlokaFT4KFfUVOFN0RnLXs2cDx4/z3ySBfi5VXkJtU61OhZuUBAQEAOHhOg587DFewNaEEl+z+sxCUkFSS1SEMebM4SU3RcSJbRAKtzXFxTyz7J57eACqDOzP2Y+iq0W4Pe52rX3bt3NPxg03tOMEU6bwJQSWLQOU+mM0R0eMBgCT3Aq33cbfxUy3YVqKyuuoAnfiBLdudVb4HDeO3zsT0sJm9JkBAPjprLT1dO67j0eiffyxpOYCCyMUbmtWruTP9c8+K9spvk/7Hp4unrg5RjO2l4gr3MmTATe3dpyAMeCll/j66gaeHaP9oxHoGYhDF6Ur3JgYHju6aVM7xtcJUCvcgcEDNbY3NPDqa3pXZ1KnhZ06xbNfJNAnsA9i/GPw0zlpCjcsjBe0WbuWlwkRWBehcNWUlgLvv8/NuH79ZDmFUqXEpjObMC1mGjxdNKMf0tN53VKL5Fj87W/AgAHA0qV6rVzGGEZHjDbJwlV3/ccfJtU973SkFKYgzCcMfh6aa+ecPs2zxbX8t625+24ejmjCah7TY6djT/YeVNVXSWr/73/zQJZXXpF8CoGFEApXzbvv8tV4ZZTCzembcaXmCu4acJfWvq1b+btFFK6TE7BkCbdyDUxJjw4fjXMl51B8tVhy17fdxq3xH3+0wDg7KCmFKRgQPEBru3rCzOD6ox4ePBxx61ZJqdoAML3PdDQoG/D7hd8ltY+PB/7xD/5Ad/CgpEMEFkIoXID7bles4HGsAwcab28GRIQ3D7yJGP+YFr/btX3AF18Ao0bxaACLMHs20L+/QSvXHD/uddfxkDWxVpZu6prqkFqYikEh2lo1KQnw9gaio4108sQTPA3t9dclnXNs5Fj4e/jj61TpmWpvvcWnKe6/36RCZYJ2IhQuALz9No9IT0yU7RQ7L+xEUkESnh3zLJydnDX2nTrFHzf//ncLntDJCXj5Ze6r0OPLHdp9KNyc3VrSjKXAGI/J3bXLpNT/TsOpK6fQpGrCsLBhWvuSkrh1azRT3M+PW7nffQecMR5FonBS4OHBD2PTmU24UCYtM8XHh0egXbjAA1sEVoKI5Hg5Dvn5RJ6eRPfcI+tpJmyYQN3f6U51jXVa+xYsIHJzIyoutvBJlUqi+HiiHj2I6rTPS0Q0+bPJ1P9//U3q9uhRIoBo3TpLDFIvcsmmrHL74ZEPCYmgnPIcje21tUSurkRPPy2xo6IiIi8vojvvlNT8UuUlcn3VlRb+vNCk8T7zDL+X27ebdJhAPwZlTFi4L73EZzJktG6PXDyCPdl78PSop+Gm0AxBKCvjlsbdd/P4TIvi5MSt9+xs4MMPdTa5qfdNOF10GpcqL0nudsgQXhhdJEFo8+elPxHsFYyILhEa248d41EK118vsaPAQB4LvnEjX03UCN19uuPegfdi3Yl1uFItfUbz1Ve55+nBB8UTizXo3Ao3OZkvn/P44xIca+bz5oE34efuh4cHP6y1Tx2JtmiRTCe/8UbgppuA117T+Y2a0nsKAEgqZq1G7VbYuVOEFrWGiLAvex/GRY3TKir/xx/8fYx2NU79PPssX9bpqackxeU+P/Z5NCgbsOzAMsmncHcHPv+cL423cKEJYxOYhzET2MyX/aNSEU2cSBQQQFRaKttpTheeJiSCluxeorWvqIjIx4do1izZTs85dYrIyYlo8WKtXSqViiLejaCZX880qcsjR/ij6Pr1FhqjNg7nUsgszSQkgv731/+09k2dShQXZ0ana9fyD/qbbyQ1n79lPrm96kZ5FXkmnea11/hpvv7ajDEKWmNQxuxScK3Cli388leulPU092++nzxe86CimiKtfY88QuTsTJSaKusQOA89RKRQEKWkaO1atH0Rub/mTlX1VZK7U6m4a3jKFEsOUgOHU7j/l/R/hETQ6cLTGtvr6/k0wULT3KucpiaiQYOIQkOJysqMNs8uyyaXpS70j5/+YdJpGhuJRoyQ3f7oDAiFq0VlJVF4ONGAAUQNDbKdJrc8lxRLFfT49se19h06xD/9p56S7fSaFBXxb9PYsXwyrRV7s/YSEkEbUzea1OVLLxExRpSba8mBtuBwCvfvm/5OwW8Hk0ql0ti+fz+/15s2mdnxsWP8CeXhhyU1X/jzQlIsVdD5kvMmnebkSX6axx4zZ5CCZoTC1eKJJ7imOHRI1tMs/nUxOb/iTFllWRrbGxt58EBYGNf9VmPdOn7L167V2NykbKLu73SnKZ+bZq5mZvLuXnvNkoNswaEUrto1c+d32lEFL7/MFVm7LMd//Yt/2Lt2GW2aX5lPPm/40I2f3ail/I3xz3/ysZ46Ze5AOz1C4Wpw9CiXqAULZD1NydUS8nrdi+7ddK/Wvnfe4Z/899/LOgRtlEqi668n8vcnKijQ2PXavtd0Pg4bY/x4ouho7mKwMA6lcM8WnyUkgj766yOtfWPHEg0d2p7eiaimhigmhv9KF2m7p9ry0V8fERJBG05sMOk0JSVcPMaNk+WedgaEwm2hpoaob1/J/rD28PSOp4klMkq5oukzzcnh4ZXTptlIoNPSiNzdiW65RWMAhdWF5P6aOz2w5QGTuvv0Uy5F+/dbeqCOpXDfOfQOIRGUXZatsb2qirvOn3uuPb03k5TEg3lvvdWo8ChVShq7bix5v+FNpy6bZq6uWsXv6QbTdLWAIxRuCwsW8Ev+7TdZT5NelE4uS13owa0Pau2bOZPIw4MoK0vWIRhm5Ur+OXykaY098csT5PyKs0m+v+pqHmkxb56lB+lYCnfChgk08KOBWtt//tnCIvfBB7zDN94w2vRixUXq/k53inovigqrCyWfoqmJ6IYbuGFw7lx7BtspEQqXiK5FJUhO9TGP2sZaSvg4gfz/60/5lfk6h7BsmaxDMI5KxeOUPDw0ohbyK/PJ/TV3um/zfSZ199BD/MtZJT3IQQoOo3DLastIsVRBL+x8QWvfo4/yz6a21tze26BSEc2dywXpu++MNj966Si5vupKd2y8w6TT5OVx18LgwXqTFAW6EQqX0tK4GWYF6Vn480JCImjb2W0a2ysqiCIiZA+MkE5+PlFICFHv3txx18xTvz5FTq840dnis5K7OniQdM3FtReHUbjfpHxDSAQdzD2osV2pJOrenWj2bHN71kNtLdGoUdw1JGHi9439bxASQd+mfmvSadQGgtUiaToGnVzhlpXxyYbgYO5AlZHPTn5GSAQ9vUPTilapiObM4XN1MgdGmMahQ0QuLjyYtrGRiIguV10mz9c9afa30rWESkU0cCBRbGxLN5bAYRTurG9mUdBbQdSkbNLYfuAA/4Z99pm5PRugsJD/WHbpwrNQDNCobKRhnwyjwLcC6Ur1FZNOs3AhiVoLptGJFe7Vq3wa3cWF6I8/ZD3VprRNpFiqoAkbJlB9U33LdpWK6IUXSKrbzfqsWcMHN39+y0TMq/teJSSCdmTskNzN5s28m48/ttjIHELh5pbnktMrTvTc79qzYvPmEXl7W9zV0urkuUS9enGle/iwwaanC0+T66uuNPvb2SaFil29yn9Mg4L4Q5HAKJ1U4dbX81AAxoi+/FLWU608spJYIqORa0dSeW15y/bCQl6EDOAx63YbZrNkCbU8O6pUVNtYS7ErYylkeQhdrLgoqQuVimjCBJ5RZaHMOYdQuAu2LSDFUoVWrHVREXeRP/KIOb2aQG4ut3Q9PIzGGS77YxkhEfTynpe1rHFDpKXx7ocOlT24pyPQCRVuTQ3RjBn88lavlu00hdWF9PdNfyckgmZ+PZNqGmqIiPto33qLGx7OzkSvvKKV3GVfqFREjz/OP68nnyRSKin1Sip5v+FNwz4ZRrWN0mZ88vOJunUjioqyiDVk9wo3oySDFEsVtGCbdkz3iy/y3/q0NFN7NYMrV7hPF+DC1qRbmTYpm1rkdcjqISaFi/30Ew9vi4oi+vVXC427Y9LJFG5JCdHo0Vza/6ddRMRSfH/6e/L/rz+5LHWhJbuXtFgMhw/ziTGAG9hW+cJZAqWSZ+ABRHffTVRfT5vPbCYkgu7ffL/kx9Bjx/isfEICnyhsB3avcOd8N4c8XvOgS5WXNLZnZ3NLf84cU3tsB7W1RPfey+/f+PE8zEAHKpWKvjr1FYUsDyGfN3zo1/PStefhw9xPDxDdfjvRRWkPP50Nx1O4hdWF9EXyFzT729nUZ2Uf6v+//hS/Kp4Grx5Mw9cMp9H/N5puWH8DTf9qOq06uoouVV7iCuGvv3hFFTc32dK4mpRN9NzvzxESQcPXDG/JzMrO5o+PjPEyDVu3ynJ6eVGpiN58k4vF8OFE2dn08p6XCYmg9w+/L7mbX37hlv2NN3LPjpnYtcL9PPlzQiIocU+i1r4ZM7jClXmOVhuVimcreHnxx6sVK/TOYuZV5NF1q64j51ec6cMjH0r+Qa2r46nc7u5Efn4OKufyYlDGGJHxOpvmVH2U2nBL+hZsSd+CnIoclNWWobS2FBcrL4JA6ObVDWMjxwIAmlRNaFI1QUnKlr/zKvKQWZYJhRJ48agb/vN7A4q6KnDfXDccj1Qg2j8aDAx5lXlQOCkQGxCLhG4JSAhJQC+/Xgj1CUWDsgFeLl4I6xIGJ6ZZHrhJ1YTCmkJcqb6CrPIsnCk6g29Pf4uUwhQ8OuRRvD/1fZSXuOH114GPP+bHLFjAS8/6+Fjuw7Q6mzYBDzwAODtDteJ9zHbZhJ/ObcOOe3dgUq9JkrrYsIF38fe/8wLrbcrDSsH0I9qPJLktqCpAv4/6oX9Qf+ybt09jyaTNm/nKxv/9Ly9naxMyMnhx299+44vQvfoqMH261k2orK/EXT/che3nt+O2vrdh3cx18HX3lXSK8+eBuXP5skGLFvE10tzcjB/XCTAotzZVuIt/XYz3j7yPYK9gxAbEwt/DH37ufoj2j8bU6KkYHDpYSwlqnIQIWVs3wOeZfyMoswB/DgvFykfi4RvaCwBwvvQ8GGOI6BKBJlUT0orSkFKYgrqmOq2+AjwCEBsQC08XT9Q01iCvIg8F1QVQkUqj3YDgAVhywxKMD74D774LfPABUF8PPPQQ8J//AOHhpnxMdkxGBnDvvcCRI2iaNAG3j8zFH13KcPTho+jl10tSF6+9xhfUeOEF4I03TB6B3Src2769Db9m/IrkR5MRGxDbsv3sWWDkSKBnT+DIEcDFRbZxGoeI/3A++yxfuGzQIP73bbdpaEYVqfDu4Xfxwq4X0MuvFzbduQn9g/tLOkV9Pe/ygw/40u/ffAPExMh1QQ6DfSrcvdl7MeHTCVgwdAE+uPkDKJwU0ntXqYDdu/m3eM8eruVWrgRmzTJ6aJOqCedKziG3IheXqy/DxckFVQ1V+PPin7hUdQlXG6/CQ+GBiK4RiOgSgVDvUIT6hCLMJxxNl/tg3+8++Pln4NAhLtN33cXXaoyNNXpqx0Op5Kb7iy+CqqqweaACn08NxYf/PoSwLmFGDycCHn0U+OQT4KOPuPVvAnapcNVy++akN/H82Odbtp85A0ycyD+yI0e40rULGhuBL7/kKwBnZPCle+6/n69QPXRoy4qWf+T8gdu/ux2ltaVYMHQBXh73MgI8pa35tHUrf5qpquL6/J57gAkTgC5d5Lwwu8Ww3BrzOZj5MohSpaSEjxMo6r0outpwVZpnRKUiOn2a6PXXeRgMwKfE33uPBwtamMpKot9/56e76y6eJcZVCK8H/dJLfDidgpISohdeoCZPDyKAksKdadvjU+nUX9uM+v4aG4mmT+dJH1u2mHRWu/PhqlQqGrtuLHV/p7uG3P76K1FgIBdHu5UJpZIP9G9/4+EGAK889uijfDWJS5eosLqQHv3pUXJ6xYl8l/nS2wfflvz9zMvjWfN+frxrZ2eiMWOIEhN5JqKewImOiP35cDed2YTZG2fj89s+x73X3au7UWMjXzv86FH++v13vhgiAIwbBzz8MHeWeXhYZMD5+cDBg8CBA/yVnMytFQCIiuILJ95yC3DzzUD37hY5peNRXo78//0X9R9/iJ4XqwEAp7srcDGhN2qHxEMxYiS8YgdgYNggBHoGthxWU8Otv1On+PLqo0dLOpvdWbg7MnZg6pdT8b9b/od/DvsnCgqApUv5Q8CAAcAPPzjIk05pKbBtG7BlC/fz1tTw7T17AvHxKOwVgjWNf+LbxpOoD+uGKYPvwIw+MzCuxzi4Orsa7LqhATh8mHf7++988UwiIDSUPw1OngyMGgX4SnMVOyL25VLILs/G3z6ZhD6XG/HFuA/gXF4BlJRwIbhyhfubLlwA8vJaNB75+kI56nrUjJ+G0pG3oNwnAtXV/BGmuhoafzc18bkB9Uul4ttav5RK/n71Knf+p6fzIQBcf48cCYwdyxf8GzGiQwuH2ZScOISsz1bA/fc96HW+GJ4N/JY3OAEX/IAr3X3gGt4D/j3i4B/VFwq/vnjqlQBcquqCResH4aYZbsZ8nHalcCvqKjBi7QjUNjRgeY8z2LrJDT/8wOVowQI+SWah337r0tQEnDjBV7k8fBhISeFfCtW1uYtydyC3C3C5qzOcAoOAgAAo/bpC5ecL5u8PZy8fNLm5QOHdBZ5dAuDZNRDeXYPg4xOA2queOPSXAlt/VuD3Pc5oaFKgCS4Ij3BFSHcXBAQo4NvFGb5dneDb1Qldu3JXRNeugJcX4OnJ3729+cvHh7/b1D9uGPtQuEqVEsHLg1FaW4qZuZ7Ysu6qZgMnJ9T6BOFUVU9koRcyqRdSqR+OYhgy0Rvt/f65uADOzoBCce3l5gb07g307Qv068d/eQcNsuubaZ80NeHqib9Qfmg36tNPoz4tBa4XcuBVVo3Aq4BzG2kIXeSGy94KuG39Himbp+qbaLELhXvwIDDu4xlQRv8EEAO++AXIvAl+fsAdpXpURgAABZtJREFUd/BJo969bTBSOamt5U7pjAwgNxeNWZkoSU9CfV4WnMsq4FPdiK51ltEbSgaoWt3ptr0SA6B0bd7HWrXT/bfChcFV/f1tHZVhRpiMXiIj+dO3bqyvcBljvwIINNIsEECxxU/uGIhrN04xEU2VezCtEXIrmc7+GRi6foNyK5eFaxTG2DEiGmqTk9sYce2Oe+2OPn5L0Nk/g/Zcv/4gV4FAIBBYFKFwBQKBwErYUuF+YsNz2xpx7Y6Lo4/fEnT2z8Ds67eZD1cgEAg6G8KlIBAIBFZCdoXLGJvKGDvLGMtgjD2vYz9jjH3QvP8UY2yw3GOyFhKufTxjrIIxdrL5tcQW47Q0jLF1jLFCxliqnv12f887s9wCnVd21cgmw8Zyf9vzAuAMIBNALwCuAJIB9GvT5hYAv4AHDI8EcETOMVnrJfHaxwPYZuuxynDtNwAYDCBVz367vuedWW5NuP4OKbutrk8WGZbbwh0OIIOILhBRA4BvAMxs02YmAPW6pn8C8GWMhco8Lmsg5do7JES0H0CpgSb2fs87s9wCnVh21cglw3Ir3DAAea3+v9i8zdQ2jojU6xrFGEtmjP3CGJNWiNTxsfd73pnlFhCyKwWz7r8JRWjNQldecduwCCltHBEp15UEIIqIqhljtwDYAqAzlHC293vemeUWELIrBbPuv9wW7kUAEa3+DweQb0YbR8TodRFRJRFVN/+9HYALY8xYLn9HwN7veWeWW0DIrhTMuv9yK9yjAGIYYz0ZY64A5gL4sU2bHwHc1zzrNxJABREVyDwua2D02hljIYzxMkaMseHg96PE6iO1PvZ+zzuz3AJCdqVg1v2X1aVARE2MsccA7ACf+VxHRKcZY4827/8YwHbwGb8MAFcBPCDnmKyFxGu/HcACxlgTgFoAc6l5CtSRYYx9DT6LHcgYuwjgZQAugGPc884st0Dnll01csmwyDQTCAQCKyEyzQQCgcBKCIUrEAgEVkIoXIFAILASQuEKBAKBlRAKVyAQCKyE3JlmHQbGmBJACnhoSBOATwG8T0QqgwcKBDZCyKz9IRSudGqJKAEAGGPBAL4C0BU8Pq9dMMaciUjZ3n4EgjYImbUzhEvBDIioEMAjAB5rzjRxZoy9zRg72lwb8x8AwBhzYox9xBg7zRjbxhjbzhi7vXlfNmNsCWPsAIA7GGO9GWO/MsaOM8b+YIz1bW4XxBj7obnvo4yxMc3bx7WqRXqCMeZjo49D4AAImbUTbF130lFeAKp1bCsD0A1ckP/TvM0NwDEAPcGzcbaD/7CFNLe/vbldNoBnW/W1C0BM898jAOxu/vsrAGOb/44EcKb5758AjGn+2xuAwtafkXjZ10vIrP29hEuhfagrBk0BcJ3aEgB/bIsBMBbAd8R9ZpcZY3vaHP8tADDGvAGMBvBdc3o6wL8EADAZQL9W27s0WwYHAbzLGPsSwCYiumjRKxN0VITM2hChcM2EMdYLgBJAIbgQP05EO9q0mWakm5rmdycA5dTsb2uDE4BRRFTbZvsyxtjP4PncfzLGJhNRuqnXIeg8CJm1PcKHawaMsSAAHwP4kPjz0Q7wQh4uzftjGWNeAA4AmN3sF+sGXgxDCyKqBJDFGLuj+XjGGItv3v0bgMdanVs9CdKbiFKI6L/gj4N9ZbhUQQdByKx9IBSudDyanf2nAewEF6pXmvetBZAGIInxRedWgz89/ABeN1O97QiACj393wPgQcZYMoDTuLakySIAQ5snNtIAPNq8/UnGWGpz+1rw9ZUEgtYImbUzRLUwmWGMeROvih8A4C/wSYPLth6XQKAPIbPyIXy48rONMeYLvvrpq0JwBQ6AkFmZEBauQCAQWAnhwxUIBAIrIRSuQCAQWAmhcAUCgcBKCIUrEAgEVkIoXIFAILASQuEKBAKBlfh/oHIvkxa5CdsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 410.976x246.586 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2)\n",
    "\n",
    "draw(results_relative_tall, \"tall\", 2, ax1, prior_property_tall)\n",
    "draw(results_absolute_dangerous, \"dangerous\", 2, ax2, prior_property_dangerous)\n",
    "\n",
    "ax1.set_ylim(0,7)\n",
    "ax2.set_ylim(0, 10)\n",
    "\n",
    "w = 5.708\n",
    "fig.set_size_inches(w, w*0.6)\n",
    "plt.savefig(\"lassiterGoodman2013.png\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [],
   "source": [
    "### CLOSED SCALE (FIG. 3)\n",
    "prior_property = stats.uniform()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### SOME PRESSURE (FIG. 4)\n",
    "myclip_a = 0\n",
    "myclip_b = 1\n",
    "my_mean = 0.25\n",
    "my_std = 0.15\n",
    "a, b = (myclip_a - my_mean) / my_std, (myclip_b - my_mean) / my_std\n",
    "prior_property = stats.truncnorm(a=a, b=b, loc=my_mean, scale=my_std)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:computational_models] *",
   "language": "python",
   "name": "conda-env-computational_models-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}