{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Bayesian estimation with multivariate normal disributions\n", "---------------------------------------------------------\n", "\n", "Copyright 2016 Allen Downey\n", "\n", "MIT License: http://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import print_function, division\n", "\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.stats import multivariate_normal, wishart\n", "\n", "from itertools import product, starmap\n", "\n", "import thinkbayes2\n", "import thinkplot\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook contains a solution to [a problem posted on Reddit](https://www.reddit.com/r/statistics/comments/4csjee/finding_pab_given_two_sets_of_data/); here's the original statement of the problem:\n", "\n", ">So, I have two sets of data where the elements correspond to each other. I'm trying to find out the probability that (91.9 <= A <= 158.3) and (56.4 <= B <= 100). I know that P(91.9 <= A <= 158.3) = 0.727098 and that P(56.4 <= B <= 100) = 0.840273, given that A is a normal distribution with mean 105.5 and standard deviation 21.7 and that B is a normal distribution with mean 76.4 and standard deviation 15.4. However, since they are dependent events, P(BA)=P(A)P(B|A)=P(B)P(A|B). Is there any way that I can find out P(A|B) and P(B|A) given the data that I have?\n", "\n", "The original poster added this clarification:\n", "\n", ">I'm going to give you some background on what I'm trying to do here first. I'm doing sports analysis trying to find the best quarterback of the 2015 NFL season using passer rating and quarterback rating, two different measures of how the quarterback performs during a game. The numbers in the sets above are the different ratings for each of the 16 games of the season (A being passer rating, B being quarterback rating, the first element being the first game, the second element being the second, etc.) The better game the quarterback has, the higher each of the two measures will be; I'm expecting that they're correlated and dependent on each other to some degree. I'm assuming that they're normally distributed because most things done by humans tend to be normally distributed.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a first step, let's look at the data. I'll put the two datasets into NumPy arrays." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "16" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.array([122.8, 115.5, 102.5, 84.7, 154.2, 83.7,\n", " 122.1, 117.6, 98.1, 111.2, 80.3, 110.0,\n", " 117.6, 100.3, 107.8, 60.2])\n", "b = np.array([82.6, 99.1, 74.6, 51.9, 62.3, 67.2,\n", " 82.4, 97.2, 68.9, 77.9, 81.5, 87.4,\n", " 92.4, 80.8, 74.7, 42.1])\n", "n = len(a)\n", "n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And make a scatter plot:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF1RJREFUeJzt3X+QXXV9//HniyxRkkJ+ULLBECigQIxiRb6Qyhe4EgHB\nkminhlgGCRSdEb/U0Y4l8O232bbjGNIyVDvCSGtjdERELJAiNoEhV0D5ESCRYEKMYCBEshkSCJUY\nwibv7x/nhL273DW79557z+45r8fMTu757D33vj+EvO7nfs7nnKOIwMzMyuGAvAswM7P2ceibmZWI\nQ9/MrEQc+mZmJeLQNzMrEYe+mVmJ7Df0JX1TUrekJ2va/lzSU5L2SDqp3/OvlrRB0jpJ57SiaDMz\na8xgRvqLgXP7ta0BPg78pLZR0jRgDjANOA+4QZIyqNPMzDKw39CPiAeBl/u1rY+IDUD/QJ8N3BIR\nPRGxEdgAnJJRrWZm1qSs5/SnAJtqtjenbWZmNgz4QK6ZWYl0ZPx6m4GpNdtHpG1vIckX/TEza0BE\nNHysdLAjffHW+fva3+2zFJgrabSko4F3Ao8O9KIRUdifBQsW5F6D++f+lbF/Re5bRPNj5f2O9CXd\nDFSAQyU9DywgObD7r8AfAndJWh0R50XEWkm3AmuBN4ArIosqzcwsE/sN/Yj4iwF+dccAz/8K8JVm\nijIzs9bwgdwWqVQqeZfQUu7fyFbk/hW5b1lQXrMvkjzzY2Y2RJKINhzINTOzAnDom5mViEPfzKxE\nHPpmZiXi0DczKxGHvplZiTj0zcxKxKFvZlYiDn0zsxJx6JuZlYhD38ysRBz6ZmYl4tA3MysRh76Z\nWYk49M3MSsShb2ZWIg59M7MS2W/oS/qmpG5JT9a0TZC0XNJ6Scskjav53dWSNkhaJ+mcVhVuZmZD\nN5iR/mLg3H5t84F7I+J44D7gagBJ7wbmANOA84AbJDV8Wy8zM8vWfkM/Ih4EXu7XPBtYkj5eAnws\nfTwLuCUieiJiI7ABOCWbUs3MrFmNzulPiohugIjYAkxK26cAm2qetzltMzOzYaAjo9eJjF7HzAZp\n9WpYtgwOOADOPx+mT8+7IhsJGg39bkmdEdEtaTKwNW3fDEyted4RaVtdXV1dbz6uVCpUKpUGyzEr\nl8cfh7//e4h0uPXgg7BwIUyblm9dlr1qtUq1Ws3s9RSx/0G6pD8C/isi3ptuXwtsj4hrJV0FTIiI\n+emB3O8Cp5JM69wDvCvqvImkes1mNgj/+I/w6KN92z70IfjiF/Opx9pHEhHR8AKZ/Y70Jd0MVIBD\nJT0PLAAWAj+QdBnwHMmKHSJiraRbgbXAG8AVTnaz7PlflTVqUCP9lryxR/pmDVu5Mhnt7/snJMFX\nvuJ5/TJodqTv0DcboR5/HP77v5MDuR/9KJx4YtK+Y0cy9XPwwXDyydCR1XINGxYc+mb2pnXr4O/+\nDnbtSrbf+c7kG8Db355vXZadZkPf194xK5Bvf7s38AF+9StYsSK/emz4ceibFciWLW9te/HF9tdh\nw5dD36xATj55cG1WXj7EY1Yg8+bBb38LP/sZjB0LF17Ye4DXDHwg16yQenpg1KhkKacVS8tPzjKz\nkcfLNG0gntM3MysRh76ZWYk49M3MSsShb2ZWIg59M7MSceibFdDevb78stXnhV1mBbJzJ3z96/DT\nnyYnZ82dCxdckHdVNpx4pG9WIN/6Ftx/P+zZA6++CjfdBE89lXdVNpw49M0KZOXKwbVZeTn0zQqk\ns3NwbVZeDn2zArn4Ynjb23q3jzkGzjorv3ps+GnqgmuSPg9cnm7+W0R8TdIE4PvAUcBGYE5E7Kiz\nry+4ZtYC27f33i7x1FN9HZ6iye12iZKmA98D/hfQA/wY+CzwGWBbRCySdBUwISLm19nfoW9mNkR5\n3i5xGvBIRLweEXuA+4E/A2YBS9LnLAE+1sR7mJlZhpoJ/aeA0yVNkDQGOB+YCnRGRDdARGwBJjVf\nppmZZaHh2b6IeFrStcA9wG+BVcCeek8d6DW6urrefFypVKhUKo2WY2ZWSNVqlWq1mtnrZXbnLElf\nBjYBnwcqEdEtaTKwIiKm1Xm+5/TNzIYozzl9JB2W/nkk8HHgZmApMC99yiXAnc28h5mZZafZJZv3\nAxOBN4AvRERV0kTgVpL5/edIlmy+Umdfj/TNzIYotyWbzXLomyVeew3GjPFNzG1wfGN0sxFq0ya4\n7jp45hk49FD4zGfggx/MuyorOo/0zXJy5ZWwcWPvdkcHLF4M48fnVpKNALkeyDWzxrz8ct/AB+jp\ngTVrcinHSsShb5aDgw+GQw55a/uUKe2vxcrFoW+Wg44OuOwyOKDmX+A55yRXxTRrJc/pm+Vo69Zk\nSmfqVDjuuLyrsZHASzbNzErEB3LNzGzQHPpmZiXi0DczKxGHvplZiTj0zcxKxKFvZlYivuCamfHG\nG8kF4B54AF58Ed7//uRkMV/5s3gc+mYlt2IFfOMb8NBDyfaxx8JPf5p8CFx+eb61WfY8vWNWYtu2\nwVe/Ci+8ALt2JT/PPAMR8OMfJ98ArFg80rchW7UKfvaz5Brw550H48blXZE1at062LMnCfl9du1K\nwv7AA/u2WzE49G1I7r4bbryxd/uee+BrX4OxY/OryRp35JHJn+PHw+jRsHt3EvYdHVCpJG1WLJ7e\nsSG57ba+21u3Jgf/iqKnB370I/inf0r6+rvf5V1Rax15JHz848nVPk84ATo74Ywz4FOfgiuuyLs6\na4WmRvqSvgD8JbAXWANcCowFvg8cBWwkuTH6jubKtOGiXggWKRj/5V/gJz9JHt9/Pzz6KCxalG9N\nrXbZZXD22cmB23e/23fuKrqGR/qS3gFcCZwUESeSfIB8EpgP3BsRxwP3AVdnUagND2ed1Xd79Gg4\n7bR8asna9u1J0Ndatw5++ct86mmnqVOT+/M68Iuv2Tn9UcBYSXuBg4DNJCF/Zvr7JUCV5IPACuDS\nS2HMGHj4YZg4EebOhUmT8q4qGz099Q9c7t7d/lrMWqWp6+lL+ivgy8BOYHlEXCzp5YiYUPOc7REx\nsc6+vp6+DTt/+7fw85/3bk+ZAjfc0PcOV2Z5avZ6+g2P9CWNB2aTzN3vAH4g6SKgf5IPmOxdXV1v\nPq5UKlQqlUbLMcvE1VfDLbfA2rVw9NHJNxkHvuWpWq1SrVYze72GR/qS/hw4NyI+nW5fDMwAzgIq\nEdEtaTKwIiKm1dnfI30zsyHK885ZzwMzJL1dkoCZwFpgKTAvfc4lwJ1NvIeZmWWo2Tn9BcBc4A1g\nFXA5cDBwKzAVeI5kyeYrdfb1SN+sxm9+A/fdl0wnzZyZrJk36883RjcrgA0bYP783pVCBx0E112X\nLKU0q+Ubo5sVwO23910a+rvfwdKl+dVjxeXQNxsG/ud/Btdm1iyHvtkwUG+18hlntL0MKwFfZdNs\nGJg5E3buTK5hP2oUXHBBclkEs6z5QK6Z2QjiA7lmZjZoDn0zsxLxnL7ZIO3YkVyMrbMTjj9+8Pvt\n3g3LliVr8Y8/Hs49N7kzlVkePKdvNggrV8LChb1r6U8/Hb70JdAgZlb/4R+S/ff5kz+Ba65pTZ1W\nfJ7TN2uDm27qe/LUAw/AmjX732/Tpr6BD/DQQ7BlS7b1mQ2WQ99sP3p66of0Cy/sf9+BbsDy+uvN\n1WTWKIe+2X50dMB73tO3TYITT9z/vscck1yXv9a73gVHHZVdfWZD4Tl9s0HYsiW5QfqGDXDwwXDJ\nJckB2cHYvh1uvrn3QO5FF8G4ca2t14rLV9k0a6MdO2DsWK++sfw49M3MSsSrd8zMbNAc+mZmJeLQ\nNzMrEYe+mVmJNBz6ko6TtErSE+mfOyT9laQJkpZLWi9pmSQvTjMzGyYyWb0j6QDgBeBU4P8A2yJi\nkaSrgAkRMb/OPl69Y2Y2RMNl9c6HgWciYhMwG1iSti8BPpbRe5iZWZOyCv0LgZvTx50R0Q0QEVuA\nSRm9h5mZNanp8wolHQjMAq5Km/rP2Qw4h9PV1fXm40qlQqXe3aHNWiACVq+GjRth+nQ47ri8KzKr\nr1qtUq1WM3u9puf0Jc0CroiIj6Tb64BKRHRLmgysiIhpdfbznL7l5vrr4b77erc/9Sn4xCfyq8ds\nsIbDnP4nge/VbC8F5qWPLwHuzOA9zDLz3HN9Ax/glltg58586jFrp6ZCX9IYkoO4/1nTfC1wtqT1\nwExgYTPvYZa1l156a9vu3fDqq+2vxazdmprTj4idwGH92raTfBCYDUvTpydXynzttd62o46CyZPz\nq8msXXyVTWuJF1+E22+Hbdtgxgz48IcHdz/Zdlm7Fr7xDfj1r5MbpHzuczBlSt5Vme2fL61sw86r\nr8JnP9t3uuSii2Du3PxqMiuK4XAg16yPBx546/z4XXflU4uZ9eXQt8zV+wLnL3Vmw4ND3zJ3+unJ\nfWRrnX9+PrWYWV+e07eW2LwZfvjD3gO5H/nI8DqQazZS+UCumVmJ+ECumZkNmkPfzKxEHPpmZiXi\n0DczKxGHvplZiTj0zcxKxKFvZtYmPT2wZ0++NTR9u0QzM/v9enrgxhuTm/d0dMCf/mlyt7Y8Tlj0\nSN/MrMVuuw2WL0/Cf9euZLv/3dvaxaFvZtZijz8+uLZ2cOibmbXY4YcPrq0dHPpmZi124YUwblzv\n9uGHwwUX5FNLUxdckzQO+HfgPcBe4DLgl8D3gaOAjcCciNhRZ19fcM3MSmPnTnjkETjwQDjlFBg9\nurHXyfUqm5K+BfwkIhZL6gDGAtcA2yJikaSrgAkRMb/Ovg59M7Mhyi30JR0CrIqIY/u1Pw2cGRHd\nkiYD1Yg4oc7+Dn0zsyHK89LKRwMvSVos6QlJN0kaA3RGRDdARGwBJjXxHmZmlqFmTs7qAE4CPhcR\nj0m6HpgP9B++Dzic7+rqevNxpVKhUqk0UY6ZWfFUq1Wq1Wpmr9fM9E4n8FBEHJNu/2+S0D8WqNRM\n76yIiGl19vf0jpnZEOU2vZNO4WySdFzaNBP4BbAUmJe2XQLc2eh7mJlZtppdvfM+kiWbBwLPApcC\no4BbganAcyRLNl+ps69H+mZmQ+Qbo5uZlYhvjG5mZoPm0DczKxGHvplZiTj0zcxKxKFvZlYiDn0z\nsxJx6JuZlYhD38ysRBz6ZmYl4tA3MysRh76ZWYk49M3MSsShb2ZWIg59M7MScehbW0TA3r15V2Fm\nzdwj12y/IuDb34a77kpC/+yz4dOfhlGj8q7MrJw80reWuvdeuO022LULdu+GH/0I7rgj76rMysuh\nby312GODazOz9nDoW0t1dg6uzczao6nQl7RR0s8lrZL0aNo2QdJySeslLZM0LptSbSSaPRsOO6x3\ne9w4+MQn8qvHrOyaujG6pGeBD0TEyzVt1wLbImKRpKuACRExv86+vjF6SezaBQ8/DHv2wIwZMHZs\n3hWZjVzN3hi92dD/NXByRGyraXsaODMiuiVNBqoRcUKdfR36ZmZD1GzoNzunH8A9klZKujxt64yI\nboCI2AJMavI9zMwsI82u0z8tIl6UdBiwXNJ6kg+CWgMO57u6ut58XKlUqFQqTZZjZlYs1WqVarWa\n2es1Nb3T54WkBcBvgcuBSs30zoqImFbn+Z7eMTMbotymdySNkfQH6eOxwDnAGmApMC992iXAnY2+\nh5mZZavhkb6ko4HbSaZvOoDvRsRCSROBW4GpwHPAnIh4pc7+HumbmQ1Rrqt3muHQNzMburxX75iZ\n2Qji0DczKxGHvplZiTj0zcxKxKFvZlYiDn0zsxJx6JuZlYhD38ysRBz6ZmYl4tA3MysRh76ZWYk4\n9M3MSsShb2ZWIg59M7MSafZ2idZPBNxzDzz+OEyeDLNnw8SJeVdlZpbw9fQz9h//Abff3rs9aRLc\neCOMHp1fTWZWHL6e/jCyZw/cfXfftq1b4ZFH8qnHzKw/h36GImDv3re29/S0vxYzs3qaDn1JB0h6\nQtLSdHuCpOWS1ktaJmlc82WODB0dMHNm37bx42HGjHzqMTPrr+k5fUlfAD4AHBIRsyRdC2yLiEWS\nrgImRMT8OvsVck6/pwfuuKP3QO6cOXD44XlXZWZFkeuN0SUdASwGvgx8MQ39p4EzI6Jb0mSgGhEn\n1Nm3kKFvZtZKeR/IvR74ElCb3p0R0Q0QEVuASU2+h5mZZaTh0Jf0UaA7IlYDv+9Tx8N5M7NhopmT\ns04DZkk6HzgIOFjSd4Atkjprpne2DvQCXV1dbz6uVCpUKpUmyjEzK55qtUq1Ws3s9TI5OUvSmcBf\np3P6i0gO5F5bxgO5ZmatlPecfj0LgbMlrQdmpttmZjYM+DIMZmYjyHAc6ZuZ2TDl0DczKxGHvplZ\niTj0zcxKxKFvZlYiDn0zsxJx6JuZlYhD38ysRBz6ZmYl4tA3MysRh76ZWYk49M3MSsShb2ZWIg59\nM7MSceibmZWIQ9/MrEQc+mZmJeLQNzMrEYe+mVmJNBz6kt4m6RFJqyStkbQgbZ8gabmk9ZKWSRqX\nXblmZtaMhkM/Il4HPhQR7wf+GDhP0inAfODeiDgeuA+4OpNKR5hqtZp3CS3l/o1sRe5fkfuWhaam\ndyJiZ/rwbUAHEMBsYEnavgT4WDPvMVIV/X88929kK3L/ity3LDQV+pIOkLQK2ALcExErgc6I6AaI\niC3ApObLNDOzLDQ70t+bTu8cAZwiaTrJaL/P05p5DzMzy44isslkSf8P2AlcDlQiolvSZGBFREyr\n83x/GJiZNSAi1Oi+DYe+pD8E3oiIHZIOApYBC4Ezge0Rca2kq4AJETG/0QLNzCw7zYT+e0kO1B6Q\n/nw/Ir4saSJwKzAVeA6YExGvZFSvmZk1IbPpHTMzG/7adkZuutLnCUlL0+1CncQlaZykH0haJ+kX\nkk4tSh8lfUHSU5KelPRdSaNHct8kfVNSt6Qna9oG7I+kqyVtSP9uz8mn6sEboH+L0vpXS/qhpENq\nfjfi+1fzu7+WtDedcdjXVoj+Sboy7cMaSQtr2ofUv3ZehuHzwNqa7aKdxPVV4O70oPX7gKcpQB8l\nvQO4EjgpIk4kOR/jk4zsvi0Gzu3XVrc/kt4NzAGmAecBN0hq+CBam9Tr33JgekT8MbCB4vUPSUcA\nZ5NMK+9rm0YB+iepAlwAvDci3gv8c9o+5P61JfTTv4zzgX+vaS7MSVzpqOn0iFgMEBE9EbGD4vRx\nFDBWUgdwELCZEdy3iHgQeLlf80D9mQXckv6dbiQJzFPaUWej6vUvIu6NiL3p5sMky6yhIP1LXQ98\nqV/bbIrRv88CCyOiJ33OS2n7kPvXrpH+vr+M2gMIRTqJ62jgJUmL0ymsmySNoQB9jIjfANcBz5OE\n/Y6IuJcC9K2fSQP0ZwqwqeZ5m9O2kewy4O70cSH6J2kWsCki1vT7VSH6BxwHnCHpYUkrJH0gbR9y\n/1oe+pI+CnRHxGrg933tGMlHlDuAk4CvR8RJwGsk0wUj/kQ1SeNJRhNHAe8gGfFfRAH6th9F6w8A\nkv4vyVLr7+VdS1bSJePXAAvyrqWFOkiWv88A/gb4QaMv1I6R/mnALEnPAt8DzpL0HWCLpE6A9CSu\nrW2opVVeIBllPJZu/5DkQ6C7AH38MPBsRGyPiD3A7cAHKUbfag3Un80ky4/3OSJtG3EkzSOZZv2L\nmuYi9O9Y4I+An0v6NUkfnpA0iaQvR9Y8dyT2D5LR/H8CpJe72SPpUBroX8tDPyKuiYgjI+IYYC5w\nX0RcDPwXMC992iXAna2upVXSaYFNko5Lm2YCvwCWMvL7+DwwQ9Lb0wNEM0kOyI/0vom+3zwH6s9S\nYG66Yulo4J3Ao+0qsgl9+ifpIyRTrLPSK+TuM+L7FxFPRcTkiDgmIo4mGYS9PyK2kvTvwpHcv9Qd\nwFkAac6MjohtNNK/iGjbD8nZukvTxxOBe4H1JCsLxrezlhb07X3ASmA1ySfyuKL0keRr8zrgSZKD\nnAeO5L4BNwO/AV4n+VC7FJgwUH9IVrr8Kv1vcE7e9TfYvw0kq1qeSH9uKFL/+v3+WWBikfpHMr3z\nHWAN8BhwZqP988lZZmYl4tslmpmViEPfzKxEHPpmZiXi0DczKxGHvplZiTj0zcxKxKFvZlYiDn0z\nsxL5/3C5wq/eUFgpAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Scatter(a, b, alpha=0.7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like modeling this data with a bi-variate normal distribution is a reasonable choice.\n", "\n", "Let's make an single array out of it:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X = np.array([a, b])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And compute the sample mean" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 105.5375 76.4375]\n" ] } ], "source": [ "x̄ = X.mean(axis=1)\n", "print(x̄)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sample standard deviation" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 21.04040384 14.93640163]\n" ] } ], "source": [ "std = X.std(axis=1)\n", "print(std)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Covariance matrix" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 472.21183333 161.33583333]\n", " [ 161.33583333 237.96916667]]\n" ] } ], "source": [ "S = np.cov(X)\n", "print(S)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And correlation coefficient" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 0.4812847]\n", " [ 0.4812847 1. ]]\n" ] } ], "source": [ "corrcoef = np.corrcoef(a, b)\n", "print(corrcoef)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's start thinking about this as a Bayesian estimation problem.\n", "\n", "There are 5 parameters we would like to estimate: \n", "\n", "* The means of the two variables, μ_a, μ_b\n", "\n", "* The standard deviations, σ_a, σ_b\n", "\n", "* The coefficient of correlation, ρ.\n", "\n", "As a simple starting place, I'll assume that the prior distributions for these variables are uniform over all possible values.\n", "\n", "I'm going to use a mesh algorithm to compute the joint posterior distribution, so I'll \"cheat\" and construct the mesh using conventional estimates for the parameters.\n", "\n", "For each parameter, I'll compute a range of possible values where\n", "\n", "* The center of the range is the value estimated from the data.\n", "\n", "* The width of the range is 6 standard errors of the estimate.\n", "\n", "The likelihood of any point outside this mesh is so low, it's safe to ignore it.\n", "\n", "Here's how I construct the ranges:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(89.757197120005102, 121.31780287999489)\n", "(65.235198775056304, 87.639801224943696)\n", "(9.5161000378105989, 32.564707642175762)\n", "(6.7553975307657055, 23.117405735750815)\n", "(0.031284703359568844, 0.93128470335956881)\n" ] } ], "source": [ "def make_array(center, stderr, m=11, factor=3):\n", " return np.linspace(center-factor*stderr, \n", " center+factor*stderr, m)\n", "\n", "μ_a = x̄[0]\n", "μ_b = x̄[1]\n", "σ_a = std[0]\n", "σ_b = std[1]\n", "ρ = corrcoef[0][1]\n", "\n", "μ_a_array = make_array(μ_a, σ_a / np.sqrt(n))\n", "μ_b_array = make_array(μ_b, σ_b / np.sqrt(n))\n", "σ_a_array = make_array(σ_a, σ_a / np.sqrt(2 * (n-1)))\n", "σ_b_array = make_array(σ_b, σ_b / np.sqrt(2 * (n-1)))\n", "#ρ_array = make_array(ρ, np.sqrt((1 - ρ**2) / (n-2)))\n", "ρ_array = make_array(ρ, 0.15)\n", "\n", "def min_max(array):\n", " return min(array), max(array)\n", "\n", "print(min_max(μ_a_array))\n", "print(min_max(μ_b_array))\n", "print(min_max(σ_a_array))\n", "print(min_max(σ_b_array))\n", "print(min_max(ρ_array))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although the mesh is constructed in 5 dimensions, for doing the Bayesian update, I want to express the parameters in terms of a vector of means, μ, and a covariance matrix, Σ.\n", "\n", "`Params` is an object that encapsulates these values. `pack` is a function that takes 5 parameters and returns a `Param` object." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class Params:\n", " def __init__(self, μ, Σ):\n", " self.μ = μ\n", " self.Σ = Σ\n", " \n", " def __lt__(self, other):\n", " return (self.μ, self.Σ) < (self.μ, self.Σ)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def pack(μ_a, μ_b, σ_a, σ_b, ρ):\n", " μ = np.array([μ_a, μ_b])\n", " cross = ρ * σ_a * σ_b\n", " \n", " Σ = np.array([[σ_a**2, cross], \n", " [cross, σ_b**2]])\n", " return Params(μ, Σ)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can make a prior distribution. First, `mesh` is the Cartesian product of the parameter arrays. Since there are 5 dimensions with 11 points each, the total number of points is `11**5 = 161,051`." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mesh = product(μ_a_array, μ_b_array, \n", " σ_a_array, σ_b_array, ρ_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is an iterator. We can use `itertools.starmap` to apply `pack` to each of the points in the mesh: " ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mesh = starmap(pack, mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need an object to encapsulate the mesh and perform the Bayesian update. `MultiNorm` represents a map from each `Param` object to its probability.\n", "\n", "It inherits `Update` from `thinkbayes2.Suite` and provides `Likelihood`, which computes the probability of the data given a hypothetical set of parameters.\n", "\n", "If we know the mean is `μ` and the covariance matrix is `Σ`:\n", "\n", "* The sampling distribution of the mean, `x̄`, is multivariable normal with parameters μ and `Σ/n`.\n", "\n", "* The sampling distribution of `(n-1) S` is Wishart with parameters `n-1` and `Σ`.\n", "\n", "So the likelihood of the observed summary statistics, `x̄` and `S`, is the product of two probability densities:\n", "\n", "* The pdf of the multivariate normal distrbution evaluated at `x̄`.\n", "\n", "* The pdf of the Wishart distribution evaluated at `(n-1) S`." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class MultiNorm(thinkbayes2.Suite):\n", " \n", " def Likelihood(self, data, hypo):\n", " x̄, S, n = data\n", "\n", " dist_x̄ = multivariate_normal(hypo.μ, hypo.Σ/n)\n", " dist_S = wishart(n-1, hypo.Σ)\n", " \n", " return dist_x̄.pdf(x̄) * dist_S.pdf((n-1) * S)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can initialize the suite with the mesh." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "suite = MultiNorm(mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And update it using the data (the return value is the total probability of the data, aka the normalizing constant). This takes a minute or two on my machine (with m=11)." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1min 21s, sys: 20 ms, total: 1min 21s\n", "Wall time: 1min 21s\n" ] }, { "data": { "text/plain": [ "2.9895441818707563e-14" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time suite.Update((x̄, S, n))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now to answer the original question, about the conditional probabilities of A and B, we can either enumerate the parameters in the posterior or draw a sample from the posterior.\n", "\n", "Since we don't need a lot of precision, I'll draw a sample." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sample = suite.MakeCdf().Sample(300)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a given pair of values, μ and Σ, in the sample, we can generate a simulated dataset.\n", "\n", "The size of the simulated dataset is arbitrary, but should be large enough to generate a smooth distribution of P(A|B) and P(B|A)." ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def generate(μ, Σ, sample_size):\n", " return np.random.multivariate_normal(μ, Σ, sample_size)\n", "\n", "# run an example using sample stats\n", "fake_X = generate(x̄, S, 300)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function takes a sample of $a$ and $b$ and computes the conditional probabilites P(A|B) and P(B|A)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.7411764705882353, 0.883177570093458)" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def conditional_probs(sample):\n", " df = pd.DataFrame(sample, columns=['a', 'b'])\n", " pA = df[(91.9 <= df.a) & (df.a <= 158.3)]\n", " pB = df[(56.4 <= df.b) & (df.b <= 100)]\n", " pBoth = pA.index.intersection(pB.index)\n", " pAgivenB = len(pBoth) / len(pB)\n", " pBgivenA = len(pBoth) / len(pA)\n", " return pAgivenB, pBgivenA\n", "\n", "conditional_probs(fake_X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can loop through the sample of parameters, generate simulated data for each, and compute the conditional probabilities:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def make_predictive_distributions(sample):\n", " pmf = thinkbayes2.Joint()\n", "\n", " for params in sample:\n", " fake_X = generate(params.μ, params.Σ, 300)\n", " probs = conditional_probs(fake_X)\n", " pmf[probs] += 1\n", "\n", " pmf.Normalize()\n", " return pmf\n", "\n", "predictive = make_predictive_distributions(sample)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then pull out the posterior predictive marginal distribution of P(A|B), and print the posterior predictive mean:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.7403697580000659" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFBhJREFUeJzt3X+sZHV5x/H3g+4SfyCW0hJdCm3wB8VUGoMrpk2dQiO7\npnUNbQJLYhXRYiLGpLGFJhrvHyZq+UcsVkpDNaahkIjtbq0UqMuNwYIuFRbUXX74g8JCtFJJqsYW\nydM/Zu4ye+7M3HNnzsw5c877lWyYmXtm5rnnzv3cL8/5nu+JzESS1H7H1F2AJGkxDHxJ6ggDX5I6\nwsCXpI4w8CWpIwx8SeqIDQM/Iq6LiO9HxH0TtvlERDwUEfdGxG9WW6IkqQplRvifBs4b98WI2Amc\nlpkvBy4FrqmoNklShTYM/My8A/jRhE12AZ8dbPtV4PiIOKma8iRJVamih78NeHTo/uHBY5KkBvGg\nrSR1xHMreI3DwK8M3T958Ng6EeHCPZI0hcyMWV+jbODH4N8oe4H3ADdGxNnAU5n5/XEv1LTF2lZW\nVlhZWam7jHWaWJc1lWNN5TWxrmlq2rPvAJ/dc+fM733TVe8e+XjEzFkPlAj8iLge6AG/GBH/CXwI\n2ApkZl6bmV+MiDdFxMPAT4CLK6lMkhpsz74D3Hjz3fzv/z294bbHbt3CBTvPYtc5Zy6gsvE2DPzM\nvKjENpdVU44kLYdxYf/Hu15fe7CPU0UPf6n1er26SxipiXVZUznWVF4T65pU06RRfVNG8ZPEInvq\nEZFN6+FLUhnj+vTHbt3C9VdeMtf3joiFHrSVpE4qM6pfFga+JA0pczC2yX36SQx8SZ22jLNtpmXg\nS+qczYQ8LH/QrzHwJXVG2aBvS8AXGfiSOmHS2bBtDfgiA19Say37vPmqGfiSWmnSiH5ZZ9nMysCX\n1Eo33nz3use6OKofZuBLaqXhNk5XR/RFXgBFUuvs2XfgqPuGfZ+BL6l1hts5x27dUmMlzWLgS2qV\nPfsOHNXOWaa1bubNwJfUKsXRve2cZxn4klrD0f1kztKRtPRGnWDl6H49A1/SUht3gpWj+/UMfElL\nrXiCVddPrprEwJe0tIo9e0+wmszAl7R07NlPx8CXtFTs2U/PwJe0NEaFvT378gx8SUujeIDWnv3m\neOKVpKXgAdrZGfiSGq/YyvEA7XQMfEmNNqpv7wHa6Rj4khprVNjbypmegS+pkQz76hn4khrJGTnV\nM/AlNZIzcqpn4EtqHK9JOx8GvqRGGTUFU9Uw8CU1hlMw58vAl9QYHqidr1KBHxE7IuJQRDwYEZeP\n+PqLImJvRNwbEfdHxNsrr1RSq7l0wvxtuHhaRBwDXA2cCzwO7I+IPZl5aGiz9wDfzMw3R8SJwAMR\n8feZ+fO5VC2pNVzbfnHKrJa5HXgoMx8BiIgbgF3AcOAncNzg9nHAk4a9pElGBf0a+/bzUSbwtwGP\nDt1/jP4fgWFXA3sj4nHghcAF1ZQnqW0mBb1r289XVevhnwfck5nnRMRpwG0R8erM/HFxw5WVlSO3\ne70evV6vohIkLYNRYW/QH211dZXV1dXKXzcyc/IGEWcDK5m5Y3D/CiAz82ND23wB+EhmfmVw/0vA\n5Zl5d+G1cqP3k9Ruf/i+a47cNujLiQgyM2Z9nTIj/P3AyyLiVOAJ4EJgd2GbR4DfA74SEScBrwC+\nM2txktqleAbt9VdeUlMl3bRh4GfmMxFxGXAr/Wmc12XmwYi4tP/lvBb4MPCZiLhv8LQ/z8z/nlvV\nkpaOZ9DWb8OWTqVvZktH6iSXOp5NVS0dz7SVNFeGfXMY+JLmyuUSmsPAlzQ3LpfQLAa+pLkYdZDW\nsK+XgS9pLoqtHJdLqF9VZ9pKEjB66QRbOc3gCF9SpVz5srkMfEmVKR6kXVs6Qc1gS0dSJUYdpHXp\nhGZxhC+pEh6kbT4DX1IlPEjbfAa+pJkVV8E07JvJwJc0E1fBXB4GvqSZ2LtfHga+pKm5Vs5yMfAl\nTW14dO8JVs1n4Eua2vDo3lZO8xn4kqbizJzlY+BLmkqxnaPmM/AlTcV2zvIx8CVtmu2c5WTgS9o0\n2znLycCXtGm2c5aTgS9pU2znLC8DX9Km2M5ZXga+pNKKSynYzlkuBr6kUkatimk7Z7kY+JI2VAx7\ncHS/jAx8SRONCntXxVxOXsRc0kh79h3gxpvvPqpnD4b9MnOEL2kkw759HOFLWqc4G+fYrVu4YOdZ\nhv2SM/AlrVOca3/9lZfUWI2qYktH0lGca99eBr6ko3jZwvYqFfgRsSMiDkXEgxFx+ZhtehFxT0R8\nIyJur7ZMSYvg6L7dNuzhR8QxwNXAucDjwP6I2JOZh4a2OR74JPDGzDwcESfOq2BJ8+Povt3KHLTd\nDjyUmY8ARMQNwC7g0NA2FwE3ZeZhgMz8YdWFSpqfUXPuHd23T5mWzjbg0aH7jw0eG/YK4ISIuD0i\n9kfEW6sqUNJ8rZ1JW5yG6ei+faqalvlc4DXAOcALgDsj4s7MfLii15c0J8NtHHh2zr3ap0zgHwZO\nGbp/8uCxYY8BP8zMnwE/i4gvA2cC6wJ/ZWXlyO1er0ev19tcxZIqUzxI65m0zbC6usrq6mrlrxuZ\nOXmDiOcAD9A/aPsE8DVgd2YeHNrmdOCvgB3AscBXgQsy81uF18qN3k/S4lz0Z9cdCXxPsGquiCAz\nY9bX2XCEn5nPRMRlwK30e/7XZebBiLi0/+W8NjMPRcQtwH3AM8C1xbCX1DwepO2WDUf4lb6ZI3yp\nMYrLHt901btrrEaTVDXC90xbqYNGXb1K7WfgSx3j1au6y8CXOqY4DdOZOd1h4Esd4zTM7jLwpQ7Z\ns+/AUfcN+27xAihSB4xaK8cDtd3jCF/qgFHXp/VAbfc4wpdazuvTao2BL7Wc16fVGls6Usu5fILW\nGPhSh9jG6TYDX2qx4jRMdZuBL7VYsX+vbjPwpRazf69hBr7UUp5VqyIDX2op2zkqMvClFiqebGU7\nR2DgS61UHN3bzhEY+FLrOLrXOAa+1CKjLl3o6F5rDHypRYpXs3J0r2EGvtQSxVaOV7NSkYEvtYQH\narURA19qCQ/UaiOuhy8tubXLFw5zdK9RHOFLS85r1aosA19aYuMuXyiNYktHWlKj5tx7+UJN4ghf\nWlLOuddmGfjSEnLOvaZh4EtLyDn3moaBLy0h59xrGga+tGS8kpWmZeBLS2TUzBypLANfWiLOzNEs\nDHxpSTgzR7MqFfgRsSMiDkXEgxFx+YTtXhsRT0fE+dWVKAmcmaPZbRj4EXEMcDVwHvAqYHdEnD5m\nu48Ct1RdpCRn5mh2ZUb424GHMvORzHwauAHYNWK79wKfA35QYX2ScGaOqlEm8LcBjw7df2zw2BER\n8VLgLZn5KSCqK0+SM3NUlaoO2n4cGO7tG/pSBYphD7ZzNL0yq2UeBk4Zun/y4LFhZwE3REQAJwI7\nI+LpzNxbfLGVlZUjt3u9Hr1eb5MlS91RnIbpzJxuWF1dZXV1tfLXjcycvEHEc4AHgHOBJ4CvAbsz\n8+CY7T8N/HNmfn7E13Kj95PUVxzdG/bdFRFk5sydkw1H+Jn5TERcBtxKvwV0XWYejIhL+1/Oa4tP\nmbUoqetG9e0Ne82q1AVQMvNfgVcWHvubMdu+o4K6pE7zjFrNg2faSg3kGbWaBwNfahjn3GtevKat\n1BB79h3gxpvvXndRcqkqjvClhiiGPdi7V7Uc4UsNUFwJ89itW7hg51m2c1QpA19qgOJKmNdfeUmN\n1aitDHypRqP69rZxNC8GvlSTUevkeIKV5smDtlINxoW9o3vNkyN8aYFGtXDAk6u0GI7wpQUy7FUn\nR/jSgjj1UnUz8KUFceql6mZLR1qA4ujeg7Oqg4EvLUBxdG8bR3Uw8KU5c3SvpjDwpTnyylVqEgNf\nmpNRJ1c5uledDHxpToqXKXS+vepm4EtzUOzbG/ZqAgNfmgNn5aiJDHypYs7KUVMZ+FKFnJWjJjPw\npQoVD9Q6uleTGPhSRTxQq6Zz8TRpBuPWt7eVoyZyhC/NYFTYg60cNZMjfGlKxRYOuMa9ms3Al6bk\n+vZaNrZ0pCk4117LyMCXpuCZtFpGBr60SY7utazs4UslOP1SbeAIXyrB6ZdqA0f4UglOv1QbGPjS\nBvbsO3DU/ZuuendNlUizKdXSiYgdEXEoIh6MiMtHfP2iiDgw+HdHRPxG9aVK9SjOyJGW1YaBHxHH\nAFcD5wGvAnZHxOmFzb4D/E5mngl8GPjbqguV6uKMHLVFmZbOduChzHwEICJuAHYBh9Y2yMy7hra/\nC9hWZZFSHdZm5gyzZ69lVqalsw14dOj+Y0wO9HcCN89SlNQExZk5tnO07Co9aBsRvwtcDPz2uG1W\nVlaO3O71evR6vSpLkCpRPLlqbVaOtAirq6usrq5W/rqRmZM3iDgbWMnMHYP7VwCZmR8rbPdq4CZg\nR2Z+e8xr5UbvJ9Vt1GUKXRhNdYoIMjNmfZ0yLZ39wMsi4tSI2ApcCOwtFHMK/bB/67iwl5ZBMezB\nA7Vqjw1bOpn5TERcBtxK/w/EdZl5MCIu7X85rwU+CJwA/HVEBPB0Zm6fZ+FS1UaFvZcpVJts2NKp\n9M1s6aiBxq2TY9irKRbZ0pFazbBXV7i0gjpv1Gwcw15tZOCr04rr5DgbR21m4KuTRvXtPbFKbWcP\nX500qm/v9Eu1nSN8dcq4kb19e3WBga9OmHSJQvv26gpbOuqEcWFvG0dd4ghfrTduITRbOOoaA1+t\nNa5fbwtHXWVLR63lTBzpaI7w1TrOxJFGM/DVKqNWvLSNI/XZ0lFrjAt72zhSnyN8tUbxguOueCkd\nzRG+WqE49dKwl9Yz8LX0Rl2D1rCX1jPwtfSKrRx79tJoBr6Wmq0cqTwP2mqpjFsEDWzlSBsx8LUU\nJgX9Gls50mQGvhpv1Pz6YZ5FK5Vj4KtxNhrNG/DSdAx8Nc6ksPegrDQ9A1+NUpx1s8ZRvTQ7A1+1\n2mjWjYueSdUx8FWbjQ7GOutGqpaBr1pMCnvbN9J8GPiauzJz6D0YK82fSyto7gx7qRkc4Wuuxs26\nAVs30qIZ+KpMmROmnHUj1cfA19TK9OaHOetGqpeBrw1tNtiLbN1IzWDga+ZAH2a4S81VKvAjYgfw\ncfqzeq7LzI+N2OYTwE7gJ8DbM/PeKgtVdQx4qZs2DPyIOAa4GjgXeBzYHxF7MvPQ0DY7gdMy8+UR\n8TrgGuDsOdVcqdXVVXq9Xt1lrDNNXVUGedGxW7dwxkuSD/zpOyt/7Vk08ednTeU1sa4m1lSVMiP8\n7cBDmfkIQETcAOwCDg1tswv4LEBmfjUijo+IkzLz+1UXXLWm/nDH1TWvUC8zUl9ZWan0PavQxJ+f\nNZXXxLqaWFNVygT+NuDRofuP0f8jMGmbw4PHGhv4a8F5z5fv5v4fXVN3Oet866751mUrRuqezh60\nnVfro24GuaRxIjMnbxBxNrCSmTsG968AcvjAbURcA9yemTcO7h8C3lBs6UTE5DeTJI2UmTHra5QZ\n4e8HXhYRpwJPABcCuwvb7AXeA9w4+APx1Kj+fRUFS5Kms2HgZ+YzEXEZcCvPTss8GBGX9r+c12bm\nFyPiTRHxMP1pmRfPt2xJ0mZt2NKRJLVDZcsjR8SOiDgUEQ9GxOUTtnttRDwdEecPPfa9iDgQEfdE\nxNcWVVNEvCEinoqIrw/+fWCz38+Ca6plPw226Q3e9xsRcftmnltDTXPZT2Xqioj3D9736xFxf0T8\nPCJeXPZ7qqGmun73XhQReyPi3kFNby/73Jpqqms/vTgiPj9477si4oyyzx0pM2f+R/8Px8PAqcAW\n4F7g9DHbfQn4AnD+0OPfAX6hilo2UxPwBmDvtN/PImuqeT8dD3wT2Da4f2ID9tPImua1n6b5foHf\nB/6t7n01rqaaP1N/AXxk7WcHPEm/xVznZ2pkTTXvp78EPji4/cpZP09VjfCPnJyVmU8DaydnFb0X\n+Bzwg8LjQfUXYylb06gDyWWfu8ia1h6vYz9dBNyUmYcBMvOHm3juomuC+eynsnUN2w38w5TPXURN\nUN9nKoHjBrePA57MzJ+XfO6ia4L69tMZwD6AzHwA+NWI+KWSz12nqm9g1MlZ24Y3iIiXAm/JzE+x\nPtASuC0i9kfEuxZV08DrB/8L9y9D/7tU9rmLrAnq20+vAE6IiNsH7/3WTTx30TXBfPZT2boAiIjn\nATuAmzb73AXWBPV9pq4GzoiIx4EDwPs28dxF1wT17acDwPkAEbEdOAU4ueRz11nkiVcfB4b7TMOh\n/1uZ+cTgL9dtEXEwM+9YQE3/AZySmT+N/npA/0Q/SOo0qaa69tNzgdcA5wAvAO6MiNFXIF+ckTVl\n5sPUt5+G/QFwR2Y+teD3nWRUTXXtq/OAezLznIg4bfDer17A+266psz8MfXtp48CV0XE14H7gXuA\nZ6Z9sapG+Ifp/+VZc/LgsWFnATdExHeBPwI+GRFvBsjMJwb//S/gH1m/dMNcasrMH2fmTwe3bwa2\nRMQJJb+fRddU236iP3q4JTN/lplPAl8Gziz53EXXNK/9VLauNRdydOukzn01rqY6P1MXA58fvPe3\nge8Cp5d87qJrqjOj/icz35GZr8nMtwG/TP94wnT7qaKDD8/h2QMIW+kfQPj1Cdt/msFBW+D5wAsH\nt18AfAV44yJqAk4aur0d+N4038+CaqpzP50O3DbY9vn0Rxpn1LyfxtU0l/20mc8F/QPKTwLPm/Z3\nZEE11fmZ+iTwobXPPP32xAk1f6bG1VTnfjoe2DK4/S7gM7N8nmb+JRgqbAfwAPAQcMXgsUuBPxmx\n7d/xbOD/2qDYewa/tFcsqib6Zwd/Y/De/w68btJz66ypzv00uP9++rNi7gPeW/d+GlfTPPfTJup6\nG3B9mefWWVPNv3svAW4Z/OzuA3bXvZ/G1VTzfjp78PWD9Ce8HD/LfvLEK0nqiHlMXZMkNZCBL0kd\nYeBLUkcY+JLUEQa+JHWEgS9JHWHgS1JHGPiS1BH/D+4svyXt8v7iAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Cdf(predictive.Marginal(0).MakeCdf())\n", "predictive.Marginal(0).Mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then pull out the posterior predictive marginal distribution of P(B|A), with the posterior predictive mean" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.8484086164036705" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE8pJREFUeJzt3W+sZHV9x/H3F90lanEt1ZB2EWxQ3GIqxFhcU9OO2Mgu\nabONfSBLohWphUSsT7TQB4b7wEQtaSJ2rUi70fBgC4nY7rYVhbrcGAzoEmFB3YVFWmQXSsVqUk2w\nSL59MHOX2bMzd87MnJk5Z877lWwyf35n7pfD8rk/vud3fhOZiSRp+Z2y6AIkSfNh4EtSSxj4ktQS\nBr4ktYSBL0ktYeBLUkuMDPyI2B0RT0fEg+uM+UxEHImIByLigmpLlCRVocwM/wvAxcPejIjtwDmZ\n+TrgSuDGimqTJFVoZOBn5t3AT9YZsgO4uTf2W8CmiDijmvIkSVWpooe/GXii7/mx3muSpBrxoq0k\ntcSLK/iMY8Cr+56f2XvtJBHhxj2SNIHMjGk/o+wMP3p/BtkHvBcgIrYCP83Mp4d9UGbW/s911123\n8Bqs0zqbWuOy1PnPX3+AnR/5B971F5876c+866zKyBl+ROwBOsCvRcQPgeuAjd3szpsy8ysRcUlE\nPAr8HLi8suokaY727j/Irbffxy/+77mB7793x1vZcdH5c66qOiMDPzMvKzHm6mrKkaT5GxX0p27c\nwLu3v7nRYQ/V9PCXTqfTWXQJpVhntZpQZxNqhGbVuXf/QW7ee8/A95cl6NdElf2hkT8sIuf58yQJ\nRs/g+9Ux5COCrOCirTN8SUtrnKCH5vfoRzHwJS2l9Vo1RXWc1c+CgS+p0crO4tsS6usx8CU10jjt\nmmVv1ZRl4EtqJGf14zPwJdVamZm8wV6OgS+plsbpze+5/oo5VdVsBr6kWplkzbzKMfAl1cqgsLdl\nUw0DX1Kt9Ie9QV8tA1/Swg1r49ibr5aBL2mhht0Re+rGDQuoZrn5FYeSFurW2+876TUvxs6GM3xJ\nCzGojeMdsbPlDF/SQhTD/tSNGwz7GTPwJc3d3v0HB67G0WzZ0pE0d/19e++UnR9n+JLmrn9278x+\nfgx8SQtl335+bOlImpu1lTlaDANf0lwMusHKm6vmy5aOpLkozuxdmTN/zvAlzVxxGaY3WC2GM3xJ\nM1dchmnYL4aBL2nmXIZZD7Z0JM3MoFU5zu4Xxxm+pJkZtF+OFsfAlzQT7pdTP7Z0JM2E++XUjzN8\nSTPhhdr6MfAlVW7v/oMnPPdCbT0Y+JIqV2znqB4MfEmVs51TTwa+pJmynVMfpQI/IrZFxOGIeCQi\nrhnw/ssjYl9EPBARD0XE+yqvVFIjFPv3qo+RgR8RpwC7gIuBNwA7I2JLYdgHge9l5gXA24G/iQiX\nfEotU9wC2f59vZQJ5QuBI5n5OEBE3ALsAA73jUngtN7j04AfZ+YvqyxUUn2tbaHQ37sH+/d1U6al\nsxl4ou/50d5r/XYB50XEk8BB4MPVlCepCQaFvVsg109VbZeLgfsz86KIOAe4MyLemJk/Kw5cWVk5\n/rjT6dDpdCoqQdKiDNpCwbCf3OrqKqurq5V/bmTm+gMitgIrmbmt9/xaIDPzU31j/hX4RGZ+s/f8\n68A1mXlf4bNy1M+T1CzFvv1tN1y1wGqWU0SQmTHt55Rp6RwAXhsRZ0fERuBSYF9hzOPAH/QKOwM4\nF3hs2uIk1Z83WTXHyJZOZj4fEVcDd9D9BbE7Mw9FxJXdt/Mm4OPAFyPiwd5hf5mZ/zOzqiXVhjdZ\nNUepHn5mfhV4feG1z/c9fopuH19Si7hnTrN4p62kibjmvnkMfEkTKX51oe2c+jPwJU2kv3fvmvtm\nMPAlTc2wbwYDX5Jawg3OJJU2bM8cNYMzfEmlDQp7V+c0h4EvqZS9+w8ODHtX5zSHLR1JpRS3UNhz\n/RULrEaTcIYvqRS3UGg+A1/SSG6hsBwMfEkjuSPmcjDwJa2reLHWdk5zGfiS1lWc3dvOaS4DX9JQ\nzu6Xi4EvaShn98vFdfiSTjJoCwVn983nDF/SSYph7+x+ORj4kk5Q7Nu7fcLysKUj6QRuobC8nOFL\nOs5VOcvNwJd0nKtylpuBL+k4Z/fLzcCXBLhBWhsY+JIAN0hrAwNfEmA7pw0MfEknsZ2znAx8SSf1\n77WcDHxJ9u9bwsCXZP++JQx8SSewf7+8DHyp5ezft4eBL7Wc/fv2MPCllrN/3x4GvtRibqfQLqUC\nPyK2RcThiHgkIq4ZMqYTEfdHxHcj4q5qy5RUtb37D3Lz3nuOP7eds/xGfgFKRJwC7ALeATwJHIiI\nvZl5uG/MJuCzwDsz81hEvHJWBUuqRn/vHmzntEGZGf6FwJHMfDwznwNuAXYUxlwG3JaZxwAy85lq\ny5RUpeIXnbx3x1tt57RAmcDfDDzR9/xo77V+5wKnR8RdEXEgIt5TVYGSqucXnbRTVd9p+2LgTcBF\nwMuAeyLinsx8tKLPl1QhV+a0U5nAPwac1ff8zN5r/Y4Cz2Tms8CzEfEN4HzgpMBfWVk5/rjT6dDp\ndMarWNJUXJlTf6urq6yurlb+uZGZ6w+IeBHwMN2Ltk8B3wZ2ZuahvjFbgL8FtgGnAt8C3p2Z3y98\nVo76eZJm67KP7j4+wz914wb2XH/FgivSKBFBZsa0nzNyhp+Zz0fE1cAddHv+uzPzUERc2X07b8rM\nwxHxNeBB4HngpmLYS6oH2zntNXKGX+kPc4YvLdyffPjG449vu+GqBVaisqqa4XunrdQibpTWbga+\n1CJulNZuVS3LlFRTe/cf5Nbb7zuhdw/279vIGb605AaFvTdbtZOBLy2x4hYK0A17Z/ftZEtHWlKD\ndsN0zX27OcOXllAx7MGevQx8aSkVtz52N0yBgS8tHbc+1jAGvrRk3PpYwxj40pJxrxwNY+BLS8zZ\nvfoZ+NISca8crcfAl5bEoHX3Uj8DX1oSxaWY9u9VZOBLS8ClmCrDrRWkBhu0E6ZLMTWMM3ypwdz2\nWONwhi81VLGNs7YLprN7DWPgSw3kTpiahC0dqWHcCVOTMvClhnEnTE3KwJcaxOWXmoaBLzWIO2Fq\nGga+1BDF2b19e43LwJcawtm9pmXgSw3h7F7TMvClBihue+zsXpPwxiupxobtlSNNwhm+VGPulaMq\nOcOXasq9clQ1A1+qqeKqHPfK0bRs6Ug15aocVc3AlxrANo6qYOBLNVRchilVwcCXaqjYv5eqUCrw\nI2JbRByOiEci4pp1xv1ORDwXEe+qrkSpXdwzR7MyMvAj4hRgF3Ax8AZgZ0RsGTLuk8DXqi5SaotB\n32Rl/15VKTPDvxA4kpmPZ+ZzwC3AjgHjPgR8CfjvCuuTWqX45SbO7lWlMoG/GXii7/nR3mvHRcRv\nAH+cmZ8DorrypHbxy000S1VdtP000N/bN/SlKRn2qlqZO22PAWf1PT+z91q/NwO3REQArwS2R8Rz\nmbmv+GErKyvHH3c6HTqdzpglS8vJpZhas7q6yurqauWfG5m5/oCIFwEPA+8AngK+DezMzENDxn8B\n+JfM/PKA93LUz5Pa6rKP7j7e0nErBfWLCDJz6s7JyBl+Zj4fEVcDd9BtAe3OzEMRcWX37bypeMi0\nRUlt5FJMzVqpzdMy86vA6wuvfX7I2PdXUJfUKn7BiebBO22lBRu09l6aBQNfWqBi2IPtHM2OgS8t\nyKCwd+29ZsnAlxakeFetYa9ZM/ClBfGuWs2bgS8tgKtytAgGvrQA7nevRTDwpTlzv3stSqkbryRN\nb+/+g9x6+30nhL373WuenOFLc1IMe3B2r/lyhi/NSXFm/+7tb3Z2r7ky8KUFcCdMLYKBL83YWu9e\nWjR7+NKMDbpQKy2CgS/N2KDevbQItnSkGSreUWvvXovkDF+aEfe5V90Y+NKMFC/U2srRotnSkSo2\n6I5ad8NUHRj4UoUGfamJ2yeoLmzpSBUqtnFclaM6cYYvVcA2jprAGb5UAXfBVBMY+NKUivvb28ZR\nXdnSkaYwaK29N1eprpzhS1Nwrb2axMCXJlRs5XiRVnVn4EsTKn4RuWGvujPwpQn5ReRqGi/aSmMa\n9IUmzu7VBAa+VNKgm6vAXTDVHAa+VMKgPXLANfdqFgNfKmHYHjm2ctQkBr40gssvtSxcpSON4PJL\nLQsDX1pHcXZvv15NVirwI2JbRByOiEci4poB718WEQd7f+6OiN+uvlRp/pzda5mMDPyIOAXYBVwM\nvAHYGRFbCsMeA34vM88HPg78fdWFSovg7F7LpMxF2wuBI5n5OEBE3ALsAA6vDcjMe/vG3wtsrrJI\naZ6Grbd3dq+mK9PS2Qw80ff8KOsH+p8Bt09TlLRI3lylZVXpssyIeDtwOfC2YWNWVlaOP+50OnQ6\nnSpLkKY2KOxt52ieVldXWV1drfxzIzPXHxCxFVjJzG2959cCmZmfKox7I3AbsC0zfzDks3LUz5MW\nqXhH7W03XLXAaqSuiCAzY9rPKdPSOQC8NiLOjoiNwKXAvkIxZ9EN+/cMC3up7gZ9e5W0TEa2dDLz\n+Yi4GriD7i+I3Zl5KCKu7L6dNwEfA04H/i4iAnguMy+cZeFSlQbtlWMbR8tmZEun0h9mS0c1ddlH\nd7t9gmqrqpaOe+modYYtu1xj2GtZubWCWme9sPduWi0zZ/hqjVEze5dfatkZ+GqNYtifunEDe66/\nYoEVSfNlS0etUNz10tm82sgZvpbaoDaOM3u1lYGvpTXse2id2autDHwtnWEXZ/0eWrWdga+l4Np6\naTQDX403rHUDzuqlfga+Gm1Y2Bv00skMfDXWoLC3dSMN5zp8NZJhL43PwFfjGPbSZAx8NYphL03O\nwFdjGPbSdLxoq9pybb1ULQNftTMq6MGwlyZh4Kt2yuxZb9hL4zPwVSvDtjE24KXpGfiqjeJFWbcx\nlqpl4GshyvTp3cZYqpaBr7koE/D9vCgrVc/A10yNG/T27KXZMfA1E2WC3nCX5svAV+XcsliqJwNf\nlRoU9ga9VA8GvirjXjdSvRn4mtqwfr1hL9WLu2Vqaoa91AzO8DWWUatv7NdL9WXgt9i4a+RHcSsE\nqd4M/JYatnRyUmsze0n1ZeC3SFUzets2UjMZ+A1XRYh7gVVqh1KBHxHbgE/TXdWzOzM/NWDMZ4Dt\nwM+B92XmA1UW2nZV99vBmbrUNiMDPyJOAXYB7wCeBA5ExN7MPNw3ZjtwTma+LiLeAtwIbJ1RzTO3\nurpKp9OZ+PhZhPMgPzr6CK8689yxjllEyE97PuelCXU2oUawzroqM8O/EDiSmY8DRMQtwA7gcN+Y\nHcDNAJn5rYjYFBFnZObTVRc8D2t/CeYV3JNaC/y6z9Sb8h9VE+psQo1gnXVVJvA3A0/0PT9K95fA\nemOO9V6rXeCXCfHv33sfD/3kxjlWVU4x2FdW/ouVlasWXJWkpmjdRdt5zdjrPuuW1D6RmesPiNgK\nrGTmtt7za4Hsv3AbETcCd2Xmrb3nh4HfL7Z0ImL9HyZJGigzY9rPKDPDPwC8NiLOBp4CLgV2Fsbs\nAz4I3Nr7BfHTQf37KgqWJE1mZOBn5vMRcTVwBy8syzwUEVd2386bMvMrEXFJRDxKd1nm5bMtW5I0\nrpEtHUnScqhse+SI2BYRhyPikYi4ZsiYTkTcHxHfjYi7xjm2BjX+Z0Qc7L337VnVWKbOiPhIr47v\nRMRDEfHLiHhFmWNrVGedzufLI2JfRDzQq/N9ZY+tUZ11Op+viIgv9+q5NyLOK3tsjeqcy/mMiN0R\n8XREPLjOmM9ExJHev/cL+l4f/1xm5tR/6P7ieBQ4G9gAPABsKYzZBHwP2Nx7/sqyxy66xt7jx4Bf\nrbquSeosjP9D4N/neS6nrbNu5xP4K+ATa//OgR/TbXfW6nwOq7OG5/OvgY/1Hr++rn8/h9U55/P5\nNuAC4MEh728H/q33+C3AvdOcy6pm+MdvzsrM54C1m7P6XQbclpnHADLzmTGOXXSNAMF8vjBm3POx\nE/jHCY9dVJ1Qr/OZwGm9x6cBP87MX5Y8tg51Qr3O53nAfoDMfBh4TUS8quSxdagT5nQ+M/Nu4Cfr\nDDnhplZgU0ScwYTnsqp/oEE3Z20ujDkXOD0i7oqIAxHxnjGOXXSN0P2P7c7e6x+YQX3j1AlARLwE\n2AbcNu6xFZimTqjX+dwFnBcRTwIHgQ+PcWwd6oR6nc+DwLsAIuJC4CzgzJLH1qFOmN/5HGXYP8dE\n53KeN169GHgTcBHwMuCeiKhuQ/ZqDKwxMx8Ffjczn+rNAO6MiEO9386L9EfA3Zn50wXXMcqgOut0\nPi8G7s/MiyLinF49b1xQLesZWGdm/ox6nc9PAjdExHeAh4D7gecXVMt61quzTuez31RL26sK/GN0\nfzuuObP3Wr+jwDOZ+SzwbER8Azi/5LGLrvHRzHwKIDN/FBH/RPd/qWbxF2Cc83EpJ7ZJ5nUux/1Z\nxTqp2fm8HPhEr54fRMR/AFtKHluHOu+r0/nMzP8F3r/2vFfnY8BLRx1bkzrn+fdzlGPAq/uer/1z\nbGSSc1nRhYcX8cIFhI10LyD8VmHMFuDO3tiX0v2Nel6ZY2tQ40uBX+mNeRnwTeCdVddYts7euE10\nL9q9ZNxja1Bnrc4n8Fngut7jM+j+r/LpdTuf69RZt/O5CdjQe/wB4It1/Pu5Tp1zO5+9n/Ea4KEh\n713CCxdtt/LCRduJzmWVRW8DHgaOANf2XrsS+PO+MR+huwrmQeBD6x07oxM7UY3Ab/ZO6P10fwnM\nrMYx6vxTYE+ZY+tWZ93OJ/DrwNd6/84fBHbW8XwOq7OG53Nr7/1DwJeATTU9nwPrnOf5BPbQ3Xb+\nF8AP6f5fXPG/oV10w/0g8KZpzqU3XklSS8xjGZckqQYMfElqCQNfklrCwJekljDwJaklDHxJagkD\nX5JawsCXpJb4f1msJmkLLtCiAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Cdf(predictive.Marginal(1).MakeCdf())\n", "predictive.Marginal(1).Mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We don't really care about the posterior distributions of the parameters, but it's good to take a look and make sure they are not crazy." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function takes μ and Σ and unpacks them into a tuple of 5 parameters:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def unpack(μ, Σ):\n", " μ_a = μ[0]\n", " μ_b = μ[1]\n", " σ_a = np.sqrt(Σ[0, 0])\n", " σ_b = np.sqrt(Σ[1, 1])\n", " ρ = Σ[0, 1] / σ_a / σ_b\n", " return μ_a, μ_b, σ_a, σ_b, ρ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we can iterate through the posterior distribution and make a joint posterior distribution of the parameters:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def make_marginals(suite):\n", " joint = thinkbayes2.Joint()\n", " for params, prob in suite.Items():\n", " t = unpack(params.μ, params.Σ)\n", " joint[t] = prob\n", " return joint\n", "\n", "marginals = make_marginals(suite)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here are the posterior marginal distributions for `μ_a` and `μ_b`" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEPBJREFUeJzt3X+s3XV9x/Hnu/THJlPUOdlohRkQUDIFwpAsuh2VQGuW\nlukSKItRXE1JQMn+UDT74e1iok5dhkOFrg1iImujdWm3Ca0GT4yJKNXSorRQZCv9ATgmkmhmuZb3\n/jjf1uNpb++5535Pz/l++nwkN/3++PTbV097X+d7P+d8zzcyE0lSueaMOoAkabgsekkqnEUvSYWz\n6CWpcBa9JBXOopekwk1b9BGxNiKeiogdxxnz6YjYHREPRMSF9UaUJM1GP2f0dwBXTrUzIpYAZ2fm\nq4CVwG01ZZMk1WDaos/MbwHPHGfIMuAL1djvAKdFxOn1xJMkzVYdc/QLgb1d6/urbZKkMeCLsZJU\nuLk1HGM/8Iqu9UXVtqNEhB+sI0kDyMwY9Pf2W/RRfR3LJuAGYH1EXAb8NDOfmupATf4QtYmJCSYm\nJkYdY2CD5l+xZmv9YbqsWXFJX+NO1sd/HAwz+9tvGv77Nx667z94zWV/OuX+DbdcP/QMsxExcMcD\nfRR9RNwFtIDfjojHgQ8D84HMzNWZ+dWIeGtEPAr8HLhuVok01votZWkQwyrciYknmZgY7zIfpmmL\nPjOv7WPMjfXEkSTVrY45+pNGq9UadYRZ6c6/eceTbNp2gIOTz48u0AyV9Pg3zXTZN967nfV3b+Xg\nc5MnJtAMNfmxr0OcyDnziMgmz9GX5IY7vz/jkl8wbw6feefFQ0qkJrv2/WtnXfIL5s/jrk/8ZU2J\nyhIRs3ox1rdXnqQGKfmlF50xpDRqujpK/uolvv4zLE7dyBdYVatxfwfLycgzekkqnEUvSYWz6CWp\ncBa9JBXOopekwln0klQ4i16SCmfRS1LhLHpJKpxFL0mFs+glqXAWvSQVzqKXpMJZ9JJUOItekgrn\n59EXpom3CNR4GPfbAWpwntEXZqYlv2Ce/wXUUUfJL5g/r6Y0qpPf5YWZacl7e0Ad5u0Ay+XUTcG8\nRaAG5e0Ay+IZvSQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6SSqcRS9JhbPoJalwFr0kFc6il6TCWfSS\nVDiLXpIKZ9FLUuH6KvqIWBwRuyLikYi4+Rj7XxQRmyLigYh4MCLeVXtSSdJApi36iJgD3ApcCVwA\nLI+I83uG3QD8MDMvBN4EfCoi/AhkSRoD/ZzRXwrszsw9mTkJrAOW9YxJ4IXV8guB/83MX9YXU5I0\nqH6KfiGwt2t9X7Wt263AayLiALAduKmeeJKk2apreuVKYFtmvjkizga+FhGvzcyf9Q6cmJg4stxq\ntWi1WjVFkKQytNtt2u12bcfrp+j3A2d2rS+qtnW7DvgoQGb+KCL+Czgf2Np7sO6ilyQdrfckeNWq\nVbM6Xj9TN/cD50TEWRExH7gG2NQzZg9wOUBEnA6cCzw2q2SSpFpMe0afmYci4kZgC50nhrWZuTMi\nVnZ252rgI8DnI2JH9ds+kJk/GVpqSVLf+pqjz8x7gPN6tt3etfwEnXl6SdKY8cpYSSqcRS9JhbPo\nJalwFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4i16S\nCmfRS1LhLHpJKpxFL0mFs+glqXAWvSQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6SSrc3FEH0PQ273iS\nTdsOcHDy+VFH0ZjbeO921t+9lYPPTY46isaIZ/QNMEjJL5jnP+3JqI6SXzB/Xk1pNC5sgwYYpOSX\nXnTGkNJonNVR8lcvuaSmNBoXTt00zJoVfhOqPxtuuX7UETQmPKOXpMJZ9JJUOItekgpn0UtS4Sx6\nSSqcRS9Jheur6CNicUTsiohHIuLmKca0ImJbRPwgIr5Rb0xJ0qCmfR99RMwBbgXeAhwA7o+IjZm5\nq2vMacBngCsyc39EvGxYgSVJM9PPGf2lwO7M3JOZk8A6YFnPmGuBDZm5HyAzn643piRpUP0U/UJg\nb9f6vmpbt3OBl0bENyLi/oh4R10BJUmzU9dHIMwFLgbeDJwKfDsivp2Zj9Z0fEnSgPop+v3AmV3r\ni6pt3fYBT2fmL4BfRMQ3gdcBRxX9xMTEkeVWq0Wr1ZpZYkkqXLvdpt1u13a8for+fuCciDgLeAK4\nBljeM2Yj8M8RcQqwAHg98I/HOlh30UuSjtZ7Erxq1apZHW/aos/MQxFxI7CFzpz+2szcGRErO7tz\ndWbuiojNwA7gELA6Mx+aVTJJUi36mqPPzHuA83q23d6z/kngk/VFkyTVwStjJalwFr0kFc6il6TC\nWfSSVDiLXpIKZ9FLUuEsekkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4i16SCmfRS1LhLHpJKpxF\nL0mFs+glqXAWvSQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6SSqcRS9JhbPoJalwFr0kFc6il6TCWfSS\nVDiLXpIKZ9FLUuEsekkqnEUvSYWz6CWpcBa9JBXOopekwln0klS4voo+IhZHxK6IeCQibj7OuD+M\niMmIeFt9ESVJszFt0UfEHOBW4ErgAmB5RJw/xbiPAZvrDilJGlw/Z/SXArszc09mTgLrgGXHGPde\n4MvAj2vMJ0mapX6KfiGwt2t9X7XtiIg4A7gqMz8HRH3xJEmzVdeLsf8EdM/dW/aSNCbm9jFmP3Bm\n1/qialu3S4B1ERHAy4AlETGZmZt6DzYxMXFkudVq0Wq1ZhhZksrWbrdpt9u1HS8y8/gDIk4BHgbe\nAjwBfBdYnpk7pxh/B/DvmfmVY+zL6f48HW3Fmq1HltesuGSESTTu3n7TbUeWN9xy/QiTqE4RQWYO\nPFMy7Rl9Zh6KiBuBLXSmetZm5s6IWNnZnat7f8ugYSRJ9etn6obMvAc4r2fb7VOMfXcNuSRJNfHK\nWEkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4i16SCmfRS1Lh+rpgSvXbvONJNm07wMHJ50cdRWNm\n473bWX/3Vg4+NznqKCqEZ/QjMkjJL5jnP9fJoI6SXzB/Xk1pVAKbY0QGKfmlF50xpDQaJ3WU/NVL\n/PA7/YpTN2PAT6TUVPwEStXBM3pJKpxFL0mFs+glqXAWvSQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6\nSSqcRS9JhbPoJalwFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkqnEUvSYWz6CWpcBa9JBXOopek\nwln0klQ4i16SCmfRS1Lh+ir6iFgcEbsi4pGIuPkY+6+NiO3V17ci4g/qjypJGsS0RR8Rc4BbgSuB\nC4DlEXF+z7DHgD/OzNcBHwH+pe6gkqTB9HNGfymwOzP3ZOYksA5Y1j0gM+/LzGer1fuAhfXGlCQN\nqp+iXwjs7Vrfx/GLfAVw92xCSZLqM7fOg0XEm4DrgDdMNWZiYuLIcqvVotVq1RlBkhqv3W7Tbrdr\nO14/Rb8fOLNrfVG17ddExGuB1cDizHxmqoN1F70k6Wi9J8GrVq2a1fH6mbq5HzgnIs6KiPnANcCm\n7gERcSawAXhHZv5oVokkSbWa9ow+Mw9FxI3AFjpPDGszc2dErOzsztXA3wIvBT4bEQFMZualwwwu\nSepPX3P0mXkPcF7Pttu7lt8DvKfeaJKkOnhlrCQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6SSqcRS9J\nhbPoJalwFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkqXK23EjzZbd7xJJu2HeDg5POjjqIR23jv\ndtbfvZWDz02OOorkGX2dBin5BfP8JyhRHSW/YP68mtLoZGfL1GiQkl960RlDSqNRqqPkr15ySU1p\ndLJz6mZI1qzwm1QdG265ftQRdJLzjF6SCmfRS1LhLHpJKpxFL0mFs+glqXAWvSQVzqKXpMJZ9JJU\nOItekgpn0UtS4Sx6SSqcRS9JhbPoJalwFr0kFc6il6TCWfSSVDhvPNIH7wV7cvK+rypFX2f0EbE4\nInZFxCMRcfMUYz4dEbsj4oGIuLDemKM105L3PrBl8L6vKsW0jRQRc4BbgSuBC4DlEXF+z5glwNmZ\n+SpgJXDbELKOzOGSf+Lh7007dpzvA9tut0cdYVZOdP667/va5Me/ydmh+flnq5+pm0uB3Zm5ByAi\n1gHLgF1dY5YBXwDIzO9ExGkRcXpmPlV34FF64uHv8Z+fWDnqGANrt9u0Wq1RxxjYKPPXcd/XJj/+\nTc4Ozc8/W/0U/UJgb9f6Pjrlf7wx+6ttY1v0zrufHJxnlwp7MXbFmq1DPf4pc2Kox9exvf2mzkzg\nQ/dt5cFnTvysoPPsarrIzOMPiLgMmMjMxdX6B4HMzI93jbkN+EZmrq/WdwF/0jt1ExHH/8MkSceU\nmQOfafZzRn8/cE5EnAU8AVwDLO8Zswm4AVhfPTH89Fjz87MJKkkazLRFn5mHIuJGYAudd+mszcyd\nEbGysztXZ+ZXI+KtEfEo8HPguuHGliT1a9qpG0lSsw31yp7qbZZfioidEfHDiHh9RLwkIrZExMMR\nsTkiThtmhkFExLkRsS0ivl/9+mxEvK8J2Q+LiL+KiB9ExI6I+GJEzG9Y/psi4sHq633VtrHNHxFr\nI+KpiNjRtW3KvBHxoeoCw50RccVoUv/KFPn/vPo/dCgiLu4Z34T8/1DleyAiNkTEi7r2NSH/30fE\n9qqD7omI3+3aN7P8mTm0L+DzwHXV8lzgNODjwAeqbTcDHxtmhhr+DnOAA8ArmpIdOAN4DJhfra8H\n3tmg/BcAO4AFwCl0pg3PHuf8wBuAC4EdXduOmRd4DbCt+p74feBRqp+uxyz/ecCrgHuBi7u2v7oh\n+S8H5lTLHwM+2rDH/7e6lt8LfG7Q/EM7o6+ePd+YmXcAZOYvM/NZOhdX3VkNuxO4algZanI58KPM\n3Euzsp8CnBoRc4HfpHNtQ1Pyvxr4TmYezMxDwDeBtwFLGdP8mfkt4JmezVM93kuBddX3xH8Duzn6\n2pQT6lj5M/PhzNwN9L6JYhnNyP/1zDx8ocx9wKJquSmP/8+6Vk8FDv9dZpx/mFM3rwSejog7qimQ\n1RHxAuDIFbOZ+STw8iFmqMPVwF3VciOyZ+YB4FPA43QK/tnM/DoNyQ/8AHhjNfXxAuCtdH6iakr+\nw14+Rd6pLjBsiibmfzfw1Wq5Mfkj4iMR8ThwLfB31eYZ5x9m0c8FLgY+k5kX03k3zgeB3ld/x/bV\n4IiYR+fZ80vVpkZkj4gX0znrOovONM6pEfEXNCR/Zu6iM+3xNTrfnNuAQ8caeiJz1aBpeYsQEX8N\nTGbmv446y0xl5t9k5pnAF+lM3wxkmEW/D9ibmYcvV91Ap/ifiojTAaoXF348xAyztQT4XmY+Xa03\nJfvlwGOZ+ZNq6uPfgD+iOfnJzDsy85LMbAE/BR6mQfkrU+XdT+cnlMMWVduaojH5I+JddH4ivLZr\nc2Pyd7mLzvQlDJB/aEVf/ci6NyLOrTa9BfghnYur3lVteyewcVgZarAc6D4LaEr2x4HLIuI3IiLo\nPPYP0Zz8RMTvVL+eCfwZnf/o454/+PX57KnybgKuqd4J9UrgHOC7JyrkcfTm7913WCPyR8Ri4P3A\n0sw82DWuKfnP6dp3Fb/6IMmZ5x/yK8mvo3Nl7QPAV+i86+alwNfpnKFtAV48yle7j5P9BcD/AC/s\n2taI7FXWDwM76bx75U5gXsPyf5POXP02oDXujz+dJ6IDwEE6T7TXAS+ZKi/wITrvltgJXDGm+a+i\nMxf8f3Suir+7Yfl3A3uA71dfn21Y/i8DD1b9uRH4vUHze8GUJBXOWyFJUuEsekkqnEUvSYWz6CWp\ncBa9JBXOopekwln0klQ4i16SCvf/zZTtwsZFlEMAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Cdf(marginals.Marginal(0).MakeCdf())\n", "thinkplot.Cdf(marginals.Marginal(1).MakeCdf());" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here are the posterior marginal distributions for `σ_a` and `σ_b`" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAECBJREFUeJzt3W2MXGd5xvHrcrxeFUojQoSpbZK2CSQiItipcZFA7QRa\nbOeDjYhUx5EQpHVkSzbkUxtUqcpaoiq0FNXU0MTYhaRqZLe4xVsJx47kjFCqvHiLY5vExg4prt8J\nNIEmiM3WvvthJmYy3peZ3TNndu79/6SVzznz7Dn3o2fn2mee2TN2RAgAkMusbhcAACge4Q4ACRHu\nAJAQ4Q4ACRHuAJAQ4Q4ACU0Y7ra32T5v+9A4bb5s+7jtZ2wvLLZEAEC7Wpm5f13S0rEetL1c0nUR\n8S5JayXdX1BtAIBJmjDcI+JxSS+N02SlpIfqbZ+SdKXtucWUBwCYjCLW3OdLOtmwf7p+DADQJbyh\nCgAJzS7gHKclvbNhf0H92GVs80E2ADAJEeF22rca7q5/jWZQ0npJO2x/QNLLEXF+nALbqa+nDAwM\naGBgoNtldAz9G9v6B7+r4ZGLhdXS3zdLX/nkLYWdb6aM3Z1/sk3Dr42Udt3+OX16+K//uOPXsdvK\ndUkthLvthyVVJL3N9n9Luk/SHEkREVsi4tu2b7P9vKRXJd3VdhVAjys62FcsmlfY+WaSsoN91fLF\npV2vXROGe0Tc2UKbDcWUA/S+rWum7xN+Jtm5aV23S+iqItbcUVepVLpdQkdl7N+eQ+c0eOCMhkcu\n6uyrv641W4e6XVJHdHvsdu07qB27hzo2s37x1M90+B5usWnkMtfAbUfmNXf0num+Vp5F1rXwsthu\n+w1V/hQSMxpr5eVgLbx8LMsAdayVl2Omr4WXhZk7ACREuANAQoQ7ACREuANAQoQ7ACREuANAQoQ7\nACREuANAQoQ7ACTEHaroSY0f+IWp6fSHeqE7mLmjJxUd7P19M/epUGaw98/pK+U6INzRo/jAr+KU\nGex8oFd5WJZBz+MDv4rDh3rlwcwdABIi3AEgIcIdABIi3AEgIcIdABIi3AEgIcIdABIi3AEgIcId\nABIi3AEgIcIdABIi3AEgIcIdABIi3AEgIcIdABIi3AEgIcIdABIi3AEgIcIdABJqKdxtL7N91PYx\n2/eO8viv2R60/Yztw7Y/VXilAICWTRjutmdJ2ixpqaSbJK22fWNTs/WSno2IhZJulfQ3tvnPtwGg\nS1qZuS+RdDwiTkTEiKTtklY2tQlJb6lvv0XSTyLi/4orEwDQjlbCfb6kkw37p+rHGm2W9B7bZyQd\nlHRPMeUBACajqKWTpZIORMSHbV8n6VHbN0fEK80NBwYGLm1XKhVVKpWCSgCAHKrVqqrV6pTO0Uq4\nn5Z0TcP+gvqxRndJ+ktJiogf2P4vSTdKGmo+WWO4AwAu1zzx3bhxY9vnaGVZZr+k621fa3uOpDsk\nDTa1OSHp9yXJ9lxJ75b0QtvVAAAKMeHMPSIu2N4gaa9qvwy2RcQR22trD8cWSZ+T9A3bh+rf9qcR\n8T8dqxoAMK6W1twj4hFJNzQde6Bh+6xq6+7AmPYcOqfBA2c0PHKx26X0hF37DmrH7iENvzbS7VLQ\ng7hDFaXpRLD39+X9ES472Pvn9JV2LXRe3mcGpp1OBPuKRfMKPed0Unawr1q+uLTrofO4ixRdsXUN\nQdKOnZvWdbsE9Bhm7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR\n7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQ\nEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAkR7gCQEOEOAAm1FO62l9k+avuY7XvHaFOx\nfcD292w/VmyZAIB2zJ6oge1ZkjZL+oikM5L2294VEUcb2lwp6SuSPhoRp21f3amCAQATa2XmvkTS\n8Yg4EREjkrZLWtnU5k5JOyPitCRFxI+LLRMA0I5Wwn2+pJMN+6fqxxq9W9JVth+zvd/2J4oqEADQ\nvgmXZdo4zy2SPizpzZKesP1ERDxf0PkBAG1oJdxPS7qmYX9B/VijU5J+HBG/kPQL29+R9D5Jl4X7\nwMDApe1KpaJKpdJexQCQXLVaVbVandI5Wgn3/ZKut32tpLOS7pC0uqnNLkl/Z/sKSf2SfkfSl0Y7\nWWO4AwAu1zzx3bhxY9vnmDDcI+KC7Q2S9qq2Rr8tIo7YXlt7OLZExFHbeyQdknRB0paIeK7tagAA\nhWhpzT0iHpF0Q9OxB5r2vyjpi8WVBgCYLO5QBYCECHcASIhwB4CECHcASIhwB4CECHcASIhwB4CE\nCHcASIhwB4CECHcASIhwB4CECHcASIhwB4CECHcASIhwB4CECHcASKio/yAbie05dE6DB85oeORi\nt0uZFnbtO6gdu4c0/NpIt0sBxsTMHRMqOtj7+3r7x67sYO+f01fatZBHbz/LUIqig33FonmFna8b\nyg72VcsXl3Y95MGyDNqydQ1B02jnpnXdLgEYFTN3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3\nAEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhFoKd9vLbB+1fcz2veO0\ne7/tEdsfL65EAEC7Jgx327MkbZa0VNJNklbbvnGMdp+XtKfoIgEA7Wll5r5E0vGIOBERI5K2S1o5\nSrtPS/qmpB8VWB8AYBJaCff5kk427J+qH7vE9jxJH4uIv5fk4soDAExGUW+o/q2kxrV4Ah4Aumh2\nC21OS7qmYX9B/VijxZK227akqyUttz0SEYPNJxsYGLi0XalUVKlU2iwZAHKrVquqVqtTOocjYvwG\n9hWSvi/pI5LOSnpa0uqIODJG+69L+veI+NdRHouJrofpZ83WoUvbW9cs7mIl08Pt99x/aXvnpnVd\nrAQzhW1FRFsrIhPO3CPigu0NkvaqtoyzLSKO2F5bezi2NH9LOwUAAIrXyrKMIuIRSTc0HXtgjLZ/\nVEBdAIAp4A5VAEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiI\ncAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeA\nhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhFoK\nd9vLbB+1fcz2vaM8fqftg/Wvx22/t/hSAQCtmjDcbc+StFnSUkk3SVpt+8amZi9I+t2IeJ+kz0n6\nWtGFAgBa18rMfYmk4xFxIiJGJG2XtLKxQUQ8GRE/re8+KWl+sWUCANrRSrjPl3SyYf+Uxg/vNZJ2\nT6UoAMDUzC7yZLZvlXSXpA+N1WZgYODSdqVSUaVSKbIEAOh51WpV1Wp1SudoJdxPS7qmYX9B/dgb\n2L5Z0hZJyyLipbFO1hjuAIDLNU98N27c2PY5Wgn3/ZKut32tpLOS7pC0urGB7Wsk7ZT0iYj4QdtV\noCP2HDqnwQNnNDxysdullGLXvoPasXtIw6+NdLsUoOsmDPeIuGB7g6S9qq3Rb4uII7bX1h6OLZL+\nXNJVkr5q25JGImJJJwvHxIoO9v6+6X1bRNnB3j+nr7RrAe1qac09Ih6RdEPTsQcatu+WdHexpWGq\nig72FYvmFXa+Tig72FctX1za9YB2FfqGKqavrWtmVhDt3LSu2yUAXTW9X2cDACaFcAeAhAh3AEiI\ncAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeA\nhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3\nAEiIcAeAhAh3AEiIcAeAhAh3AEhodrcLwOX2HDqnwQNnNDxysdulFGbXvoPasXtIw6+NdLsUYEZo\naeZue5nto7aP2b53jDZftn3c9jO2FxZb5sxSdLD393X/BVqZwd4/p6+U6wDT2YTPetuzJG2WtFTS\nTZJW276xqc1ySddFxLskrZV0fwdqnfaq1Woh5yk62FcsmlfIuabSvzKDfdXyxZP63qLGbzrK3Dcp\nf/8mo5VlmSWSjkfECUmyvV3SSklHG9qslPSQJEXEU7avtD03Is4XXfB0Vq1WValUCj3n1jWTC6pO\nKKp/Ozetm3oxHdCJ8ZsuMvdNyt+/yWjl9fp8SScb9k/Vj43X5vQobQAAJenZN1TXbB3q2rWHnj0x\n6vFzTzynb/10Z6HXuv2e7vWz2XNPDunwSzNyxQ3oOY6I8RvYH5A0EBHL6vuflRQR8YWGNvdLeiwi\ndtT3j0r6veZlGdvjXwwAMKqIcDvtW5m575d0ve1rJZ2VdIek1U1tBiWtl7Sj/svg5dHW29stDgAw\nOROGe0RcsL1B0l7V1ui3RcQR22trD8eWiPi27dtsPy/pVUl3dbZsAMB4JlyWAQD0ntLubrH9Q9sH\nbR+w/XRZ1+0U29tsn7d9qOHYW23vtf1923tsX9nNGqdijP7dZ/uU7e/Wv5Z1s8bJsr3A9j7bz9o+\nbPsz9eMpxm+U/n26fjzL+PXbfqqeJYdt31c/3vPjN07f2h670mbutl+Q9NsR8VIpF+ww2x+S9Iqk\nhyLi5vqxL0j6SUT8Vf1O3rdGxGe7WedkjdG/+yT9b0R8qavFTZHtd0h6R0Q8Y/tXJf2navdq3KUE\n4zdO/1YpwfhJku03RcTPbV8h6T8kfUbS7coxfqP1bbnaHLsy70t3ydfrqIh4XFLzL6qVkh6sbz8o\n6WOlFlWgMfon1caxp0XEuYh4pr79iqQjkhYoyfiN0b/X7zvp+fGTpIj4eX2zX7X3DkN5xm+0vklt\njl2ZYRuSHrW93/bdJV63TG9//a+EIuKcpLd3uZ5O2FD//KCtvfiyt5nt35C0UNKTkuZmG7+G/j1V\nP5Ri/GzPsn1A0jlJj0bEfiUZvzH6JrU5dmWG+wcj4hZJt0laX3/Zn122d6u/Kum3ImKhaj94Pf3y\nvr5k8U1J99RnuM3j1dPjN0r/0oxfRFyMiEWqveJaYvsmJRm/Ufr2Hk1i7EoL94g4W//3RUn/ptpn\n1mRz3vZc6dK654+6XE+hIuLF+OWbNF+T9P5u1jMVtmerFnz/GBG76ofTjN9o/cs0fq+LiJ9Jqkpa\npkTjJ72xb5MZu1LC3fab6rMI2X6zpI9K+l4Z1+4w643rYIOSPlXf/qSkXc3f0GPe0L/6E+Z1H1dv\nj+E/SHouIjY1HMs0fpf1L8v42b769WUJ278i6Q9Ue1+h58dvjL4dnczYlfLXMrZ/U7XZeqj2BsE/\nRcTnO37hDrL9sKSKpLdJOi/pPknfkvQvkt4p6YSkP4yIl7tV41SM0b9bVVu/vSjph5LW9uInf9r+\noKTvSDqs2s9kSPozSU9L+mf1+PiN0787lWP83qvaG6az6l87IuIvbF+lHh+/cfr2kNocO25iAoCE\n0vxpIgDglwh3AEiIcAeAhAh3AEiIcAeAhAh3AEiIcAeAhAh3AEjo/wF3OxFAlFVdUAAAAABJRU5E\nrkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Cdf(marginals.Marginal(2).MakeCdf())\n", "thinkplot.Cdf(marginals.Marginal(3).MakeCdf());" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, the posterior marginal distribution for the correlation coefficient, `ρ`" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADydJREFUeJzt3W+oXHedx/H3p6YprNRCFQqmtuxW26JYi8RsH5Tdsco2\nkYWIfVBb6GIXJWWt5Jl1H4j3QcEVXLBuVtMsQfGBJGAWknVNrdgOUmk1gTb1T2JTXWKTdiutuqDQ\n9hq+++COyXhN7kwnM3Nuf3m/YOCcM78558uPuZ/53d+ZcyZVhSSpTRd0XYAkaXYMeUlqmCEvSQ0z\n5CWpYYa8JDXMkJekho0M+SQ7kzyf5MkV2nwxydEkTyS5frolSpImNc5I/ivAzWd7Mskm4Kqqehuw\nBdg+pdokSedoZMhX1SPAb1Zoshn42qDtD4BLklw2nfIkSediGnPy64BnhtZPDLZJkjrmiVdJatia\nKezjBPCWofXLB9v+TBJvlCNJE6iqTPK6cUM+g8eZ7AM+DuxOcgPw26p6/mw78oZoSxYWFlhYWOi6\njFXBvjitlb64Zeu5f//ip499k7ff8PdTqGa+9tx319T3mUyU78AYIZ/k60APeGOSXwKfAdYCVVU7\nqupbST6Q5Gng98CdE1cjaar2PnSI3fsP8vIri53VMGnoLSz8LwsL0w/M883IkK+q28doc/d0ypE0\nTV0H/EVrL+zs2FriideO9Hq9rktYNeyL06bdF10H/K2b1k/8et8X05F5zpEnKefkpfkZnhufxVyx\n5iPJxCdeHclLUsMMeUlqmCEvSQ0z5CWpYYa8JDXMkJekhhnyktQwQ16SGmbIS1LDDHlJapghL0kN\nM+QlqWGGvCQ1zJCXpIYZ8pLUMENekhpmyEtSw0b+xqukc7Mafkxb5y9H8tKMrYaA9we1z1+GvDRj\nqyHgz+UHtfXa5nSNNEf+mLbmzZG8JDXMkJekhhnyktQwQ16SGmbIS1LDDHlJapghL0kNM+QlqWGG\nvCQ1zJCXpIYZ8pLUMENekhpmyEtSw8YK+SQbkxxJ8lSSe87w/BuS7EvyRJIfJfnI1CuVJL1qI0M+\nyQXANuBm4B3AbUmuXdbs48BPqup64L3AvybxNsaS1LFxRvIbgKNVdayqFoFdwOZlbQq4eLB8MfBi\nVf1hemVKkiYxTsivA54ZWj8+2DZsG/D2JM8Ch4Ct0ylPknQupjWlcjPweFXdlOQq4DtJrquq3y1v\nuLCwcGq51+vR6/WmVIIktaHf79Pv96eyr1TVyg2SG4CFqto4WP8UUFX1uaE23wQ+W1XfH6x/F7in\nqg4u21eNOp7Umlu2bj+17M//aRJJqKpM8tpxpmsOAG9NcmWStcCHgX3L2hwD3j8o5jLgauAXkxQk\nSZqekdM1VXUyyd3Agyx9KOysqsNJtiw9XTuAe4GvJnly8LJPVtWvZ1a1JGksY83JV9UDwDXLtt0/\ntPwcS/PykqRVxCteJalhhrwkNcyQl6SGGfKS1DBDXpIaZshLUsMMeUlqmCEvSQ0z5CWpYYa8JDXM\nkJekhhnyktQwQ16SGmbIS1LDDHlJapghL0kNM+QlqWGGvCQ1zJCXpIYZ8pLUMENekhq2pusCpHnZ\n+9Ahdu8/yMuvLHZdijQ3juR13ug64C9ae2Fnx9b5y5DXeaPrgL910/rOjq/zl9M1Oi/tue+urkuQ\n5sKRvCQ1zJCXpIYZ8pLUMENekhpmyEtSwwx5SWqYIS9JDTPkJalhhrwkNcyQl6SGjRXySTYmOZLk\nqST3nKVNL8njSX6c5OHplilJmsTIe9ckuQDYBrwPeBY4kGRvVR0ZanMJ8O/A31XViSRvmlXBkqTx\njTOS3wAcrapjVbUI7AI2L2tzO7Cnqk4AVNUL0y1TkjSJcUJ+HfDM0PrxwbZhVwOXJnk4yYEkd0yr\nQEnS5KZ1q+E1wLuBm4DXA48mebSqnp7S/iVJExgn5E8AVwytXz7YNuw48EJVvQS8lOR7wLuAPwv5\nhYWFU8u9Xo9er/fqKpakxvX7ffr9/lT2lapauUHyOuBnLJ14fQ74IXBbVR0eanMt8G/ARuAi4AfA\nrVX102X7qlHHk2bllq3bTy37oyF6LUlCVWWS144cyVfVySR3Aw+yNIe/s6oOJ9my9HTtqKojSb4N\nPAmcBHYsD3hJ0vyNNSdfVQ8A1yzbdv+y9c8Dn59eaZKkc+UVr5LUMENekhpmyEtSwwx5SWqYIS9J\nDTPkJalhhrwkNcyQl6SGGfKS1DBDXpIaZshLUsMMeUlqmCEvSQ0z5CWpYYa8JDXMkJekhhnyktQw\nQ16SGmbIS1LDDHlJapghL0kNM+QlqWGGvCQ1zJCXpIYZ8pLUMENekhpmyEtSwwx5SWqYIS9JDTPk\nJalhhrwkNcyQl6SGGfKS1DBDXpIatqbrAnT+2fvQIXbvP8jLryx2XYrUPEfymruuA/6itRd2dmxp\n3sYK+SQbkxxJ8lSSe1Zo954ki0k+NL0S1ZquA/7WTes7O740byOna5JcAGwD3gc8CxxIsreqjpyh\n3b8A355FoWrTnvvu6roEqWnjjOQ3AEer6lhVLQK7gM1naPcJ4BvAr6ZYnyTpHIwT8uuAZ4bWjw+2\nnZLkzcAHq+rLQKZXniTpXEzrxOsXgOG5eoNeklaBcb5CeQK4Ymj98sG2YeuBXUkCvAnYlGSxqvYt\n39nCwsKp5V6vR6/Xe5UlS1Lb+v0+/X5/KvtKVa3cIHkd8DOWTrw+B/wQuK2qDp+l/VeA/6qq/zzD\nczXqeGrfLVu3n1r2xKs0WhKqaqIZkpEj+ao6meRu4EGWpnd2VtXhJFuWnq4dy18ySSGSpOkb64rX\nqnoAuGbZtvvP0vYfp1CXJGkKvOJVkhpmyEtSwwx5SWqYIS9JDTPkJalhhrwkNcyQl6SGGfKS1DBD\nXpIaZshLUsMMeUlqmCEvSQ0z5CWpYYa8JDXMkJekhhnyktQwQ16SGmbIS1LDDHlJapghL0kNM+Ql\nqWGGvCQ1zJCXpIYZ8pLUMENekhpmyEtSwwx5SWqYIS9JDTPkJalhhrwkNcyQl6SGGfKS1DBDXpIa\nZshLUsMMeUlq2JquC1B39j50iN37D/LyK4tdlyJpRsYaySfZmORIkqeS3HOG529PcmjweCTJO6df\nqqat64C/aO2FnR1bOl+MDPkkFwDbgJuBdwC3Jbl2WbNfAH9TVe8C7gX+Y9qFavq6DvhbN63v7PjS\n+WKc6ZoNwNGqOgaQZBewGTjyxwZV9dhQ+8eAddMsUrO35767ui5B0gyMM12zDnhmaP04K4f4R4H9\n51KUJGk6pnriNcl7gTuBG8/WZmFh4dRyr9ej1+tNswRJes3r9/v0+/2p7CtVtXKD5AZgoao2DtY/\nBVRVfW5Zu+uAPcDGqvr5WfZVo46n+bll6/ZTy07XSKtXEqoqk7x2nOmaA8Bbk1yZZC3wYWDfsgKu\nYCng7zhbwEuS5m/kdE1VnUxyN/AgSx8KO6vqcJItS0/XDuDTwKXAl5IEWKyqDbMsXJI02lhz8lX1\nAHDNsm33Dy1/DPjYdEuTJJ0rb2sgSQ0z5CWpYYa8JDXMkJekhhnyktQwQ16SGmbIS1LDDHlJapgh\nL0kNM+QlqWGGvCQ1zJCXpIYZ8pLUMENekhpmyEtSwwx5SWqYIS9JDTPkJalhhrwkNWys33jVbO19\n6BC79x/k5VcWuy5FUmMcya8CXQf8RWsv7OzYkmbLkF8Fug74Wzet7+z4kmbL6ZpVZs99d3VdgqSG\nOJKXpIYZ8pLUMENekhpmyEtSwwx5SWqYIS9JDTPkJalhhrwkNcyQl6SGecXrEG8UJqk1juSHdB3w\n3ihM0rQZ8kO6DnhvFCZp2saarkmyEfgCSx8KO6vqc2do80VgE/B74CNV9cQ0C503bxQmqQUjR/JJ\nLgC2ATcD7wBuS3LtsjabgKuq6m3AFmD7DGptSr/f77qEVcO+OM2+OM2+mI5xRvIbgKNVdQwgyS5g\nM3BkqM1m4GsAVfWDJJckuayqnp+kqPPhBGi/36fX63VdxqpgX5xmX5xmX0zHOHPy64BnhtaPD7at\n1ObEGdqMreuA9wSopFasyhOvXQe8J0AltSJVtXKD5AZgoao2DtY/BdTwydck24GHq2r3YP0I8LfL\np2uSrHwwSdIZVVUmed04c/IHgLcmuRJ4DvgwcNuyNvuAjwO7Bx8Kvz3TfPykRUqSJjMy5KvqZJK7\ngQc5/RXKw0m2LD1dO6rqW0k+kORplr5Ceedsy5YkjWPkdI0k6bVrJidek2xMciTJU0nuOUubLyY5\nmuSJJNfPoo7VYFRfJLk9yaHB45Ek7+yiznkY530xaPeeJItJPjTP+uZpzL+RXpLHk/w4ycPzrnFe\nxvgbeUOSfYOs+FGSj3RQ5swl2Znk+SRPrtDm1edmVU31wdIHx9PAlcCFwBPAtcvabAL+e7D818Bj\n065jNTzG7IsbgEsGyxvP574Yavdd4JvAh7quu8P3xSXAT4B1g/U3dV13h33xz8Bn/9gPwIvAmq5r\nn0Ff3AhcDzx5lucnys1ZjORPXTxVVYvAHy+eGvYnF08BlyS5bAa1dG1kX1TVY1X1f4PVxziH6wtW\nuXHeFwCfAL4B/Gqexc3ZOH1xO7Cnqk4AVNULc65xXsbpiwIuHixfDLxYVX+YY41zUVWPAL9ZoclE\nuTmLkJ/7xVOr2Dh9MeyjwP6ZVtSdkX2R5M3AB6vqy0DL38Qa531xNXBpkoeTHEhyx9yqm69x+mIb\n8PYkzwKHgK1zqm21mSg3vZ/8KpHkvSx9K+nGrmvp0BeA4TnZloN+lDXAu4GbgNcDjyZ5tKqe7ras\nTtwMPF5VNyW5CvhOkuuq6nddF/ZaMIuQPwFcMbR++WDb8jZvGdGmBeP0BUmuA3YAG6tqpX/XXsvG\n6Yv1wK4kYWnudVOSxaraN6ca52WcvjgOvFBVLwEvJfke8C6W5q9bMk5f3Al8FqCqfp7kf4BrgYNz\nqXD1mCg3ZzFdc+riqSRrWbp4avkf6T7gH+DUFbVnvHiqASP7IskVwB7gjqr6eQc1zsvIvqiqvxo8\n/pKlefl/ajDgYby/kb3AjUlel+QvWDrRdnjOdc7DOH1xDHg/wGAO+mrgF3Otcn7C2f+DnSg3pz6S\nLy+eOmWcvgA+DVwKfGkwgl2sqg3dVT0bY/bFn7xk7kXOyZh/I0eSfBt4EjgJ7Kiqn3ZY9kyM+b64\nF/jq0FcLP1lVv+6o5JlJ8nWgB7wxyS+BzwBrOcfc9GIoSWrYqrwLpSRpOgx5SWqYIS9JDTPkJalh\nhrwkNcyQl6SGGfKS1DBDXpIa9v+21GBpMB+8dgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.Cdf(marginals.Marginal(4).MakeCdf());" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can ignore everything after this, which is my development code and some checks." ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "ename": "Exception", "evalue": "YouShallNotPass", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mException\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"YouShallNotPass\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mException\u001b[0m: YouShallNotPass" ] } ], "source": [ "raise Exception(\"YouShallNotPass\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def estimate(X):\n", " return X.mean(axis=1), np.cov(X)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "estimate(generate(x̄, S, n).transpose())" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def z_prime(r):\n", " return 0.5 * np.log((1+r) / (1-r))\n", "\n", "def sampling_distributions(stats, cov, n):\n", " sig1, sig2, _ = std_rho(cov)\n", " \n", " array = np.zeros((len(stats), 8))\n", " for i, (x̄, S) in enumerate(stats):\n", " array[i, 0:2] = x̄\n", " s1, s2, r = std_rho(S)\n", " array[i, 2] = s1\n", " array[i, 3] = s2\n", " array[i, 4] = r\n", " array[i, 5] = (n-1) * S[0, 0] / cov[0, 0]\n", " array[i, 6] = (n-1) * S[1, 1] / cov[1, 1]\n", " array[i, 7] = z_prime(r)\n", " return array\n", "\n", "dists = sampling_distributions(stats, cov, n)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cdf0 = thinkbayes2.Cdf(dists[:, 0])\n", "cdf1 = thinkbayes2.Cdf(dists[:, 1])\n", "thinkplot.Cdfs([cdf0, cdf1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cdf2 = thinkbayes2.Cdf(dists[:, 2])\n", "cdf3 = thinkbayes2.Cdf(dists[:, 3])\n", "thinkplot.Cdfs([cdf2, cdf3])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cdf4 = thinkbayes2.Cdf(dists[:, 4])\n", "thinkplot.Cdfs([cdf4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cdf5 = thinkbayes2.Cdf(dists[:, 5])\n", "cdf6 = thinkbayes2.Cdf(dists[:, 6])\n", "thinkplot.Cdfs([cdf5, cdf6])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cdf7 = thinkbayes2.Cdf(dists[:, 7])\n", "thinkplot.Cdfs([cdf7])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def sampling_dist_mean(i, mean, cov, cdf):\n", " sampling_dist = scipy.stats.norm(loc=mean[i], scale=np.sqrt(cov[i, i]/n))\n", " xs = cdf.xs\n", " ys = sampling_dist.cdf(xs)\n", " thinkplot.plot(xs, ys)\n", " thinkplot.Cdf(cdf)\n", " \n", "sampling_dist_mean(0, mean, cov, cdf0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sampling_dist_mean(1, mean, cov, cdf1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def sampling_dist_std(i, mean, cov, cdf):\n", " sampling_dist = scipy.stats.chi2(df=n)\n", " xs = cdf.xs\n", " ys = sampling_dist.cdf(xs)\n", " thinkplot.plot(xs, ys)\n", " thinkplot.Cdf(cdf)\n", " \n", "sampling_dist_std(5, mean, cov, cdf5)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sampling_dist_std(6, mean, cov, cdf6)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def sampling_dist_r(i, mean, cov, cdf):\n", " _, _, rho = std_rho(cov)\n", " sampling_dist = scipy.stats.norm(loc=z_prime(rho), scale=1/np.sqrt(n-3))\n", " xs = cdf.xs\n", " ys = sampling_dist.cdf(xs)\n", " thinkplot.plot(xs, ys)\n", " thinkplot.Cdf(cdf)\n", " \n", "sampling_dist_r(7, mean, cov, cdf7)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pdf_X = scipy.stats.multivariate_normal(mean, cov/n)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pdf_X.pdf(mean) - pdf_X.pdf(mean-0.1) " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def make_multi_norm_marginal(index, mean, cov, n):\n", " \n", " sigmas = std_rho(cov)\n", " width = 6 * sigmas[index] / np.sqrt(n)\n", " xs = np.linspace(mean[index]-width/2, mean[index]+width/2, 101)\n", " array = np.tile(mean, (len(xs), 1))\n", " array[:, index] = xs\n", "\n", " pdf_X = scipy.stats.multivariate_normal(mean, cov/n)\n", " ys = pdf_X.pdf(array)\n", " \n", " pmf = thinkbayes2.Pmf(dict(zip(xs, ys)))\n", " pmf.Normalize()\n", " return pmf\n", "\n", "pmf = make_multi_norm_marginal(0, mean, cov, n)\n", "thinkplot.Pdf(pmf)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pmf = make_multi_norm_marginal(1, mean, cov, n)\n", "thinkplot.Pdf(pmf)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def generate_statistics(mean, cov, n, iters):\n", " return [estimate(generate(mean, cov, n)) for _ in range(iters)]\n", "\n", "stats = generate_statistics(mean, cov, n, 1000)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "s0 = np.zeros(len(stats))\n", "s1 = np.zeros(len(stats))\n", "\n", "for i, (x̄, S) in enumerate(stats):\n", " sigmas = std_rho(S)\n", " s0[i] = sigmas[0]\n", " s1[i] = sigmas[1]\n", "\n", "thinkplot.Scatter(s0, s1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "s0 = np.zeros(len(stats))\n", "s1 = np.zeros(len(stats))\n", "\n", "for i, (x̄, S) in enumerate(stats):\n", " s0[i] = (n-1) * S[0][0]\n", " s1[i] = (n-1) * S[1][1]\n", "\n", "thinkplot.Scatter(s0, s1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pdf_S = wishart(df=n-1, scale=cov)\n", "stats = pdf_S.rvs(1000)\n", "\n", "s0 = np.zeros(len(stats))\n", "s1 = np.zeros(len(stats))\n", "\n", "for i, S in enumerate(stats):\n", " s0[i] = S[0][0]\n", " s1[i] = S[1][1]\n", "\n", "thinkplot.Scatter(s0, s1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sigmas = std_rho(cov)\n", "width = 6 * sigmas[0] / np.sqrt(2 * (n-1))\n", "X = np.linspace(sigmas[0]-width/2, sigmas[0]+width/2, 101)\n", "\n", "width = 6 * sigmas[1] / np.sqrt(2 * (n-1))\n", "Y = np.linspace(sigmas[1]-width/2, sigmas[1]+width/2, 101)\n", "Z = np.zeros((len(X), len(Y)))\n", "\n", "pdf_S = wishart(df=n-1, scale=cov)\n", "\n", "for i, x in enumerate(X):\n", " for j, y in enumerate(Y):\n", " S = cov.copy()\n", " S[0, 0] = x**2\n", " S[1, 1] = y**2\n", " try:\n", " density = pdf_S.pdf((n-1) * S)\n", " Z[i, j] = density\n", " except:\n", " Z[i, j] = np.nan" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "thinkplot.Scatter(s0, s1)\n", "plt.contour(X, Y, Z)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "pmf_0 = thinkbayes2.Pmf()\n", "\n", "for i, (x̄, S) in enumerate(stats):\n", " sig1, sig2, rho = std_rho(S)\n", " density = pdf_S.pdf((n-1) * S)\n", " pmf_0[sig1] += 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "thinkplot.Cdf(pmf_0.MakeCdf())\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "pdf_S = wishart(df=n-1, scale=cov)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pdf_S.pdf(cov)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }