From 92c245e86749a901a24377bb323d9961935f82df Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sun, 29 Apr 2018 13:22:49 -0400 Subject: [PATCH] More on 5.9.8 --- BDA 5.9.8.ipynb | 97 ++++++++++++++++++++++++++++++------------- Useful Formulae.ipynb | 2 +- 2 files changed, 68 insertions(+), 31 deletions(-) diff --git a/BDA 5.9.8.ipynb b/BDA 5.9.8.ipynb index d3ea631..5228ce3 100644 --- a/BDA 5.9.8.ipynb +++ b/BDA 5.9.8.ipynb @@ -46,70 +46,107 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In the special case under consideration, $p_1$ is normal with mean $-1$ and $\\sigma=.5$, $p_2$ is normal with mean $1$ and $\\sigma=.5$ and we can set $\\lambda_1=.1$ and $\\lambda_2=.9$. The $p_m(\\{y_{i}\\})$ can be calculated from the $t$ distribution. Drawing a sample of size $10$ from $p_1$ and getting a sample mean of $-.25$ and a sample variance of $1$ gives a $t$-statistics of $\\sqrt{10}(-.25+1)$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. " + "In the special case under consideration, $p_1$ is normal with mean $-1$ and $\\sigma=.5$, $p_2$ is normal with mean $1$ and $\\sigma=.5$ and we can set $\\lambda_1=.1$ and $\\lambda_2=.9$. The $p_m(\\{y_{i}\\})$ can be calculated from the $t$ distribution. Drawing a sample of size $10$ from $p_1$ and getting a sample mean of $-.25$ and a sample variance of $1$ gives a $t$-statistics of \n", + "$$\\frac{(\\overline{y}-\\mu)}{s/\\sqrt{N}}=\\frac{(-.25+1)}{1/\\sqrt{10}}=\\sqrt{10}(-.25+1)$$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. " ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy.stats import norm, t\n", + "import matplotlib.pyplot as plt\n", + "t_1=np.sqrt(10)*.75\n", + "t_2=np.sqrt(10)*(-1.25)\n", + "p1y=t.pdf(t_1,df=9)\n", + "p2y=t.pdf(t_2,df=9)\n", + "lambda1,lambda2=(.1,.9)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "wt1=lambda1*p1y/(lambda1*p1y+lambda2*p2y)\n", + "wt2=lambda2*p2y/(lambda1*p1y+lambda2*p2y)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "0.041664931082753924\n", - "0.0035119750957915393\n" + "0.600590082806086 0.3994099171939141\n" ] } ], "source": [ - "import numpy as np\n", - "from scipy.stats import norm, t\n", - "t_1=np.sqrt(9)*.75\n", - "t_2=np.sqrt(9)*1.25\n", - "print(t.pdf(t_1,df=9))\n", - "print(t.pdf(t_2,df=9))" + "print(wt1, wt2)" ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 26, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "0.5655172413793104" - ] - }, - "execution_count": 18, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], + "source": [ + "def posterior(prior_mean,prior_variance,sample_mean,pop_variance,n):\n", + " post_var=1/((1/prior_variance) + n/pop_variance)\n", + " post_mean=(prior_mean/prior_variance+sample_mean*n/pop_variance)/(1/post_var)\n", + " return post_mean, post_var" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], "source": [ - ".1*.041/(.1*.041+.9*.0035)" + "post_mean1,post_var1=posterior(-1,.25,-.25,1,10)\n", + "post_mean2,post_var2=posterior(1,.25,-.25,1,10)" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 28, "metadata": {}, "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1+PHPSQDD6gYiAhpUUMBElriwWFFcICCIIoJVwd1atOpXf6Kt2tL67dd9qUu1VbSWgiCCKCAILkgBJW6ALMoSMIKK7IgsCc/vj5OBEBIyJPfOvXfmvF+veU1mcufeM4Gc3Dn3ec4jzjmMMcYkl7SgAzDGGOM9S+7GGJOELLkbY0wSsuRujDFJyJK7McYkIUvuxhiThCy5G2NMErLkbowxSciSuzHGJKFqFW0gIi8BPYEfnXMnlfF9AZ4EcoGtwCDn3GcV7bd+/fouMzPzgAM2xphU9umnn/7knGtQ0XYVJnfgZeBp4F/lfL870Lz4dhrwXPH9fmVmZpKXlxfH4Y0xxsSIyIp4tquwLOOcmw6s288mvYF/OTUbOEREGsUXpjHGGD94UXNvDHxb4nFB8XPGGGMC4kVylzKeK7PVpIhcLyJ5IpK3Zs0aDw5tjDGmLPHU3CtSADQt8bgJsKqsDZ1zLwAvAOTk5FivYWOS2M6dOykoKGDbtm1BhxJJGRkZNGnShOrVq1fq9V4k9/HAYBEZiV5I3eicW+3Bfo0xEVZQUEDdunXJzMxEB9WZeDnnWLt2LQUFBTRr1qxS+4hnKOQIoAtQX0QKgPuB6sUB/B2YiA6DXIIOhbyqUpEYY5LKtm3bLLFXkohw+OGHU5XydYXJ3Tk3oILvO+C3lY7AGJO0LLFXXlV/dl6UZYwJxNy5MGkStGgBvXtDms23NmY3+3UwkfTKK9CuHQwZAhddBBdcANu3Bx2Viar77ruPqVOnBh2Gpyy5m8hZtAiuuw7OOgu+/x6eegomToTbbgs6MhNFRUVFDB06lHPOOeeAXhN2ltxN5Pzud1C7NgwfDg0bws03w513wnPPwXvvBR2dCZP8/HxOPPFEBg4cSHZ2Nn379mXr1q1kZmYydOhQOnfuzOjRoxk0aBCvv/46ANOmTaNt27ZkZWVx9dVXs734I2Hp14Sd1dxNpHzxBUyZAg89BEccsef5oUNh9Gi49VbdxurvIRP7h/FSmzbwxBMVbrZ48WJefPFFOnXqxNVXX82zzz4L6DjyGTNmAPDOO+8AOsJn0KBBTJs2jRYtWnDllVfy3HPPceutt+7zmrCzXwETKX/7G9SqBddeu/fzGRnwwAMwbx6MHx9MbCacmjZtSqdOnQC4/PLLdyfnSy+9dJ9tFy9eTLNmzWjRogUAAwcOZPr06bu/X9ZrwsrO3E1kbNumZ+eXXgqHHrrv9/v1gz/8AR58UEfP2Ci8EInjDNsvpYcUxh7Xrl17n211ZHf5ynpNWNmZu4mMKVNg82ZN7mWpVg3+539g9my9GQOwcuVKZs2aBcCIESPo3LlzudueeOKJ5Ofns2TJEgBeffVVzjzzzITE6TVL7iYyxozRM/azzy5/m4EDoU4d+Oc/ExeXCbeWLVvyyiuvkJ2dzbp16/jNb35T7rYZGRkMGzaMSy65hKysLNLS0rjxxhsTGK13rCxjIsE5PXPv1g3210epTh09sx8xAh5/HOrVS1yMJpzS0tL4+9//vtdz+fn5ez1++eWXd3/dtWtXPv/88332U/o1YWdn7iYSFi7UMe1du1a87bXXwtat8Npr/sdlTFhZcjeRMG2a3seT3E87DVq3hmHD/I3JhF9mZibz588POoxAWHI3kfDee9CsGcSzproIXH45zJoFEfskbYxnLLmb0Nu1Cz78cP8XUkvr31/vR470JyZjws6Suwm9JUtg/Xro2DH+12RmQocOemHVmFRkyd2E3ief6P2ppx7Y6wYM0LbACxZ4H5MxYWfJ3YTeJ59oo7CWLQ/sdf36aY8ZO3s3VTFu3DgWVOIMYfz48fzf//2fDxHFx5K7Cb1PPoGcHEhPP7DXNWyodfoRI3ScvDGVUZnkXlhYSK9evRgyZMgBvcZLltxNqO3YAZ9/fuAlmZjLLoOlS2HOHG/jMtFQXsvf8tr6DhkyhFatWpGdnc0dd9zBzJkzGT9+PHfeeSdt2rRh6dKlLF26lG7dutG+fXvOOOMMFi1aBMCgQYO4/fbbOeuss7jrrrt4+eWXGTx4MAArVqyga9euZGdn07VrV1auXFnma7xkM1RNqM2frwk+J6dyr+/TB268Uc/eK/sHwlRdgB1/92n5+9hjj/H888/v09b3yiuvZOzYsSxatAgRYcOGDRxyyCH06tWLnj170rdvX0BnsP7973+nefPmfPzxx9x00028V7yQwNdff83UqVNJT0/fa9br4MGDufLKKxk4cCAvvfQSt9xyC+PGjdvnNV6yM3cTavPm6f3JJ1fu9YccAj166JDICCyeY3xQuuXvtGnTymzrW69ePTIyMrj22mt54403qFWr1j772rJlCzNnzuSSSy6hTZs23HDDDaxevXr39y+55JIyk/SsWbO47LLLALjiiiv26glf3muqys7cTajNm6e92o8/vvL7uOwyGDsWPvggvhmuxnsBdvzdp+VveapVq8Ynn3zCtGnTGDlyJE8//fTuM/KYXbt2ccghh/BFOR9D4m0JXDImv9oI25m7CbW5c6FVqwO/mFpSjx5Qty785z/exWWio3TL33POOafMtr5btmxh48aN5Obm8sQTT+xO4HXr1mXz5s0A1KtXj2bNmu1eZs85x5dffllhDB07dmRk8Yy64cOH77ftsFcsuZtQmzcPsrKqto+aNbX2PmYMFF83MymkdMvf2267rcy2vps3b6Znz55kZ2dz5pln8vjjjwPQv39/Hn74Ydq2bcvSpUsZPnw4L774IieffDKtW7fmzTffrDCGp556imHDhpGdnc2rr77Kk08+6ffbRipaecQvOTk5Li8vL5Bjm2j46Sdo0AAeeUQX4aiKyZO1XfDYsXDhhd7EZ/Zv4cKFtDzQyQkey8/Pp2fPnpFtHlbWz1BEPnXOVTjEwM7cTWjFLqZW9cwdtNbeoIGVZkzqsORuQiuW3LOzq76vatV0EY+33oJNm6q+PxMN1vLXmBCaOxfq19eZpl4YMEAX2R471pv9mYoFVfZNBlX92VlyN6G1aJGOlIlzJFuFOnSAY4+FEnNLjI8yMjJYu3atJfhKcM6xdu1aMjIyKr0PG+duQmvxYh3l4hURuO46uPtu+PprKJ7DYnzSpEkTCgoKWLNmTdChRFJGRgZNmjSp9OstuZtQWrdOR8t4nYAHDYJ774V//AMeftjbfZu9Va9enWbNmgUdRsqysowJpa+/1vsTTvB2v0ceCb17a2nGxrybZBZXcheRbiKyWESWiMg+PSxF5GgReV9EPheRuSKS632oJpXEkrsfpZPrr9dPBcV9m4xJShUmdxFJB54BugOtgAEi0qrUZn8ARjnn2gL9gWe9DtSklsWLteXAscd6v+9zztHFtp9+2vt9GxMW8Zy5nwoscc4tc87tAEYCvUtt44B6xV8fDKzyLkSTir7+WhN79ere7zstDW65BWbMgI8/9n7/xoRBPMm9MfBticcFxc+V9EfgchEpACYCN3sSnUlZixd7X28v6ZprtB3wI4/4dwxjghRPci9rlHHpgasDgJedc02AXOBVEdln3yJyvYjkiUieDY8y5dm1C775xt/kXreuLuLxxhu6UpMxySae5F4ANC3xuAn7ll2uAUYBOOdmARlA/dI7cs694JzLcc7lNGjQoHIRm6RXUKAzSf0eh37LLdqW4MEH/T2OMUGIJ7nPAZqLSDMRqYFeMB1fapuVQFcAEWmJJnc7NTeVsnix3vt55g7QqJGOnBk2zM7eTfKpMLk75wqBwcBkYCE6KuYrERkqIr2KN/sf4DoR+RIYAQxyNufYVJKfwyBLu+cePXsfOtT/YxmTSHHNUHXOTUQvlJZ87r4SXy8AOnkbmklVS5ZA7do64chvjRrBb38Ljz+ubQlOPNH/YxqTCDZD1YTOihWQmeldw7CK3HWXrtZ0//2JOZ4xiWDJ3YROfj4cc0zijtegAdx2G4waBZ99lrjjGuMnS+4mdPLz9cw9ke64Aw47TGvwxiQDS+4mVDZtgvXrE5/cDz5Ya+6TJ8MHHyT22Mb4wZK7CZUVK/Q+kWWZmN/+Fho31iRvY71M1FlyN6GSn6/3iT5zhz0XVWfP1rVWjYkyS+4mVGJn7kEkd4CrrtLx9ffcA0VFwcRgjBcsuZtQyc/XM+igulNUqwZ/+Qt89RX85z/BxGCMFyy5m1CJDYNM1Bj3slx8MbRrB/fdZ6s1meiy5G5CJYhhkKWlpcEDD2gsI0YEG4sxlWXJ3YTKihXBjJQp7fzzISsLHnssoiNnnNO/TgsWwC+/BB2NCYAldxMaW7bo2qZBn7mDloVuvx3mzYOpU4OO5gDs3AkPPQRNm+pagq1b66okAwbsGYpkUoIldxMaQY+UKW3AAGjYEJ54IuhI4lRQAB06aLOcVq3guedg+HBdleStt/S5114LOkqTIHF1hTQmEYKcwFSWgw6Ca6+Fv/4VvvtOJziF1ooVcNZZ+tFnzBi46KI937vsMu2vcNll0L+/TgO+7rrgYjUJYWfuJjSCnMBUnquu0mX/Xnkl6Ej2Y/NmuOACWLcOpk3bO7HHNG0K774L3bvrmfykSYmP0ySUJXcTGvn5erbcsGHQkexx3HHQpQu89FKIL6zecINeOB09Gk45pfztMjK09WV2Nvz611rGMUnLkrsJjdhImbSQ/a+85hpdhu+//w06kjK88YaO17z/fjj33Iq3r1NHE/yOHTBoUIj/YpmqCtmvkUllie7jHq8LL9ST3tGjg46klHXr4De/0RlXQ4bE/7rmzeHRR7WEYwP5k5YldxMaYZjAVJY6dbRU/frrWn8Pjb/8RS+gvvQSVK9+YK+97jot4dxxh9bsTdKx5G5C4Zdf4Mcfw5ncAfr2hVWrtGNkKCxdCk8/rVd8Tz75wF+flqavX70aHnzQ+/hM4Cy5m1AI2zDI0nr2hBo19Ow9FO69V8/Whw6t/D5OPRX69YMnn9RPACapWHI3oRDGYZAl1aun1yvHjQvBNchvvtHJSIMHw1FHVW1f998PP/+sNXiTVCy5m1AIe3IHyM2F5cs1twbqoYf0rP2226q+r1atdGLT3/6m6xuapGHJ3YTCihWarxo1CjqS8nXvrveBzv/57judUXX11XDkkd7s86679Oz9H//wZn8mFCy5m1DIz4ejjw7fGPeSmjWDE06Ad94JMIinntIhO3fe6d0+Tz4Zzj5bz9537vRuvyZQIf5VMqkkrMMgS+vWDT74IKAuutu2wYsvQu/e+pfGS7fdpjNWx4zxdr8mMJbcTSiEpY97Rbp31xw7fXoAB3/9dVi7VicueS03Vyc3Pfmk9/s2gbDkbgK3bZsOt47CmXvnzrrO6gcfBHDw557TBHz22d7vOy1N/2jMng3z53u/f5NwltxN4Fau1PsoJPfatXV4eMKT+9y5MHOmdnT068LEFVfoVe0XX/Rn/yahLLmbwIV9AlNpZ54JeXm6clTCPP+8NrgZNMi/Y9Svr/X8V1+1lcGTgCV3E7gojHEv6cwzobBQT6QTYscObfDVpw8cdpi/x7rmGq3rjx/v73GM7yy5m8Dl52sdu6qTLROlUydIT4cPP0zQASdO1AlGV1zh/7HOPVcX9rDSTOTFldxFpJuILBaRJSJSZm9REeknIgtE5CsR+Y+3YZpklp+v+aRaRBZ9rFMHcnISmNxffRWOOCK+fu1VlZ4OAwfClCl6ldtEVoXJXUTSgWeA7kArYICItCq1TXPgbqCTc641cKsPsZokFZVhkCV17qx19x07fD7Q+vXw9tu6Wnei/vpddpk20Bk1KjHHM76I58z9VGCJc26Zc24HMBLoXWqb64BnnHPrAZxzP3obpklmUZnAVNLpp+s1xy+/9PlAo0frX5BElGRiWrbUWasjRybumMZz8ST3xsC3JR4XFD9XUgughYj8V0Rmi0g3rwI0yW3HDu2THsXkDgno7/7vf8OJJ+pqS4nUv7++ueXLE3tc45l4kruU8VzppqfVgOZAF2AA8E8ROWSfHYlcLyJ5IpK3Zs2aA43VJKFvv9UKQNTKMk2aQOPGPif3lSvho490MWsp69fQR/376/1rryX2uMYz8ST3AqBpicdNgFVlbPOmc26nc245sBhN9ntxzr3gnMtxzuU0aNCgsjGbJBK1YZAlnX66z8n9jTf0/tJLfTxIOTIzoUMHW2M1wuJJ7nOA5iLSTERqAP2B0oNgxwFnAYhIfbRMs8zLQE1yinpyX7ZMlwf0xeuvQ3a2thwIQv/+OjN2wYJgjm+qpMLk7pwrBAYDk4GFwCjn3FciMlREehVvNhlYKyILgPeBO51za/0K2iSPFSt0Nn3j0ldxIiBWd//4Yx92vnq1zpK6+GIfdh6nfv20HDR6dHAxmEqLa5y7c26ic66Fc+4459wDxc/d55wbX/y1c87d7pxr5ZzLcs7ZZXYTl/x8rV9Xrx50JAeuXTsdnThrlg87HztWL0YEmdyPPBI6dtRYTOTYDFUTqCgOg4ypVUurJnPm+LDz11/XUTKtWlW8rZ/69NHxnjZqJnIsuZtARXECU0nt2sFnn3m8aPaaNTr99eKLEz9KprQ+ffR+3Lhg4zAHzJK7CczOnbr4T1TP3EGT+7p1e9oWe+LNN3Upvb59PdxpJR17rH48sdJM5FhyN4EpKNAcFuXk3r693n/2mYc7HTNGk+rJJ3u40yro0wdmzPBxWJDxgyV3E5jYMMgol2WysrTXlmfJfcsWeO89uPDC4EsyMX36aN3J2gBHiiV3E5jYIh1RPnOvWVOveXqW3KdO1Z4MPXt6tEMPZGfrgtxWmokUS+4mMPn5enLatGmFm4Za7KKqJyZMgHr1tO1kWIjo2fvUqbBpU9DRmDhZcjeByc/XBTpq1Ag6kqpp1w6+/96D9ufO6cIc550XvoH/vXvrJ4p33w06EhMnS+4mMCtWRLskExNr2Fjls/cvvtAWmT16VDkmz3XsCIccop8sTCRYcjeBifIEppLatNHKxaefVnFHb7+tO+re3ZO4PFWtGnTrpsl9166gozFxsORuAlFYqO1+kyG516kDJ5zgQXKfMAFOOQUaNvQkLs/17KnDIfPygo7ExMGSuwnEqlVQVBTtYZAltWmjDRQrbc0a+OSTcJZkYrp10y5vVpqJBEvuJhBRbvVbluxsfU8bN1ZyB5Mm6QXVMCf3ww/XHu9vvx10JCYOltxNIJIxuQPMn1/JHUyYoF0Y27b1LCZf9OypV45XlV6vx4SNJXcTiNgEpqiPcY+JJfdKlWZ27oTJkyE3V8seYRb7ZDFxYrBxmAqF/H+SSVb5+dCoEWRkBB2JN5o00ZGClUruM2dqPSfMJZmYk06Co4+20kwEWHI3gUiWYZAxInr2XqnkPmGCTlo691zP4/KciP4Revdd2LYt6GjMflhyN4HIz0+ekTIx2dkwb14leru//TaceSbUretLXJ7r2RO2btWe8ya0LLmbhCsqSp4x7iVlZcHmzXuuJ8Rl+XJYuDAaJZmYs87SjmlWmgk1S+4m4Vav1muIyZbcK3VRNTZmPErJvWZN6NpVk7unS1AZL1lyNwmXDH3cy3LSSXp/wMm9eXO9RUmPHvoPuWhR0JGYclhyNwmXDH3cy1KnDhx33AEk959/hvffj9ZZe0xurt7bbNXQsuRuEm75cr0/+uhg4/DDAY2Yee892L49msn96KP1o4qNdw8tS+4m4fLzdTJmrVpBR+K97Gz45hsdTFKhCRP0dP9Xv/I9Ll/06AEffVSFngvGT5bcTcItX558JZmY7GztiLtgQQUbOqfJ/bzzortaSW6utvecOjXoSEwZLLmbhFu+XJfkTEZZWXpfYWlm7lwoKIhmSSamY0c4+GCru4eUJXeTULE+7sma3I89VstNFSb3WEKMXZiMomrV4Pzzte5uC3iEjiV3k1DffacJPlnLMunp0Lq1zlTdrwkToH17vfgQZT16wA8/wOefBx2JKcWSu0mo2EiZZD1zhz0jZsqd37N2LcyeHe2STEy3btpvxkozoWPJ3SRUKiT3rCz46Sc9oS3TO+9oGSMZkvsRR+jSgDYkMnQsuZuEys/XE71k6eNellgbgnJLMxMmaFLMyUlYTL7q0UOXCFyzJuhITAmW3E1CLV+uvc+jOvovHrERM2Um98JCXVIvCgtzxCs3V2tQkyYFHYkpIa7/XSLSTUQWi8gSERmyn+36iogTkSQ5JTFeS+ZhkDH16+t10jJHzMyaBRs2JEdJJqZdO2jY0EozIVNhcheRdOAZoDvQChggIq3K2K4ucAvwsddBmuSRbIt0lCfW230fEyboEMIoLMwRr7Q0PXufPFk/mZhQiOfM/VRgiXNumXNuBzAS6F3Gdn8GHgJseRZTpu3bdShksp+5g5ZmFiwoI9dNmABnnKGTf5JJbq5+Ipk1K+hITLF4kntj4NsSjwuKn9tNRNoCTZ1z++3eLyLXi0ieiOStsYsvKWflSi3Npkpy37YNliwp8eTKlTB/fnKVZGLOPVc/kdiQyNCIJ7lLGc/tHsErImnA48D/VLQj59wLzrkc51xOgwYN4o/SJIXYMMhUKctAqdJMFBfmiNfBB+snEqu7h0Y8yb0AKDlwrQmwqsTjusBJwAcikg+cDoy3i6qmtNgiHalw5t6ypc5W3Su5v/02HH88nHBCYHH5KjdX3/DKlUFHYogvuc8BmotIMxGpAfQHxse+6Zzb6Jyr75zLdM5lArOBXs65PF8iNpG1fLl+cm/cuOJtoy4jQxdX2j1i5uefYdo0XVxayvownARin0js7D0UKkzuzrlCYDAwGVgIjHLOfSUiQ0Wkl98BmuSxfLmu8ZCeHnQkibHXiJnYwhw9ewYak69OPFFrbpbcQ6FaPBs55yYCE0s9d18523apelgmGeXnp0ZJJiYrC0aNgi1boM7bb0PdulqXTlYievY+bJheTc7ICDqilJYkU+RMFCxblnrJHWD+PKf19vPPT+6puaDJfetW+PDDoCNJeZbcTUJs3KitR5o3DzqSxNk9YmbCSli1KrlLMjFdukDNmjYkMgQsuZuEWLpU748/Ptg4EumYY3SJ1Lnv/qAli+7dgw7JfzVrwtlna3Ivt+exSQRL7iYhvvlG71MpuaelaWlm3sJ0OO007QSZCnJztQb39ddBR5LSLLmbhIjN1DzuuGDjSLSs47Yyb3MmrkcKlGRiYkMirTQTKEvuJiGWLIGjjoLatYOOJLGydn3JOg5n1akXBh1K4hxzjK41aEMiA2XJ3STEkiWpVZKJyf5Wz17nFe3TSDW55ebC9OmweXPQkaQsS+4mIVIyuW/fTtanLwMwb36SzkotT48esHMnTJ0adCQpy5K78d2WLfD99ymY3N9/n0O3fkfjw38pe+GOZNaxozYTs7p7YCy5G9+l4jBIAMaOhTp1yM6pUf56qsmqenU47zytu9uQyEBYcje+i42USankXlQEb74J3buTdXI6CxdqlSKl9OgBq1fDF18EHUlKsuRufJeKY9yZPRt++AH69CErC3bsSMFh39266f3b+13Dx/jEkrvx3ZIlun5y3bpBR5JAY8dqaSI3t+yFO1JBw4Zw+ukwblzQkaQkS+7Gdyk3UsY5Te5du8LBB3PiidrHPuWSO8BFF8Fnn8GKFUFHknIsuRvfpVxynzdPp9/36QNoI8gTTiD1RszA7p8BY8cGG0cKsuRufPXzz/DddymW3MeO1UZhvXvvfmqvhTtSyfHH65t/442gI0k5ltyNrxYv1vuWLYONI6HGjoVOnbTmXCwrSysTGzcGGFdQLroIZszQC8wmYSy5G18tXKj3J54YbBwJs2wZfPnlnnJEsdhF1fnzA4gpaH366HWIN98MOpKUYsnd+GrhQl0zNWUW6Xj9db0vldxjqzKlZGkmK0vbgVrdPaEsuRtfLVqkv9fJvrrcbiNHau/2UusJNm2qs/FT8qKqiJZmpk2DDRuCjiZlWHI3vlq4MIVKMl9/DZ9/Dpdeus+3RKBNGx0VmJIuukin6FqvmYSx5G58U1ios1NT5mLqa69pFu/Xr8xvt2+v5fjCwgTHFQannqoN/W3UTMJYcje+WbpUT9ZSJrmPHAlnnAGNG5f57XbtYNs2WLAgwXGFQVqaXoeYNMl6vCeIJXfjm0WL9D4lyjLz52vW7t+/3E3at9f7lC3N9O8Pv/wC48cHHUlKsORufJNSwyBHjtSz04svLneTFi2gTh349NMExhUmHTvqleURI4KOJCVYcje+WbhQy6wHHxx0JD5zTpN7165wxBHlbpaWBm3bpnByT0vTs/fJk2Ht2qCjSXqW3I1vUmakzMcf6wWGAQMq3LR9e21vnpIXVUF/RoWFe+YDGN9Ycje+2LULvvoKTjop6EgS4OWXoVYt6Nu3wk3btdOyc+x6RMpp00a7qFlpxneW3I0vli2DrVv3TLtPWr/8oiWZiy+Oq2F97KJqypZmRPTsffp07ShnfGPJ3fgiNhPz5JODjcN3b76p3cAGDYpr8xNOgNq1U3jEDGhyd07nBRjfWHI3vpg7V6+ftWoVdCQ+e+UVOPpo6NIlrs3T07UyMWeOv2GFWosWkJOjPztbPNs3cSV3EekmIotFZImIDCnj+7eLyAIRmSsi00TkGO9DNVEyd642C6tVK+hIfPTddzBlClx5pf4li9Ppp+uZ+44dPsYWdlddpf9JPv886EiSVoX/I0UkHXgG6A60AgaISOnzsc+BHOdcNvA68JDXgZpomTs3Bert//qXXjkeOPCAXnb66bB9u46aSVkDBsBBB8GwYUFHkrTiOd04FVjinFvmnNsBjAR6l9zAOfe+c25r8cPZQBNvwzRRsmWLjgxM6uReVATPPw9nn33Ay0x16KD3s2b5EFdUHHqoNhMbPlx7MhjPxZPcGwPflnhcUPxcea4BJlUlKBNtsQUpkjq5T5qkSyvddNMBv7RxY2jSJMWTO2hpZv16W8TDJ/EkdynjuTKvgojI5UAO8HA5379eRPJEJG/NmjVNjywxAAAR/UlEQVTxR2kiJVZGTeqRMs8+C40aQa9elXp5hw4we7bHMUXN2WfrxWgrzfginuReADQt8bgJsKr0RiJyDvB7oJdzbntZO3LOveCcy3HO5TRo0KAy8ZoImDMH6tfX39uktGwZvPMOXH89VK9eqV106KAn/qtXexxblKSn6xDSKVPg228r3NwcmHiS+xyguYg0E5EaQH9gr7ZuItIWeB5N7D96H6aJkrw8OOUUna+SlJ5/XkfHXHddpXdx+ul6n/Klmdj8gH/+M9AwklGFyd05VwgMBiYDC4FRzrmvRGSoiMQ+kz4M1AFGi8gXImI9PVPUzz9r24GcnKAj8cnPP2si6t273L7t8WjXTpceTPnk3qwZdO8OL7yQ4mNDvVctno2ccxOBiaWeu6/E1+d4HJeJqC++0NGBp5wSdCQ+GTYM1q2D22+v0m4OOkgXJ/rwQ4/iirKbb9YEP2ZMXM3XTHxshqrxVGzmZVKeuRcWwqOPal/yTp2qvLsuXbTHzKZNVQ8t0s47T2e8/e1vQUeSVCy5G0/NmaPVikaNgo7EB2PGQH4+/L//58nuunTRTzkzZniyu+hKS4Pf/lZrVCnbUc17ltyNp+bMSdKzdufgoYe089cFF3iyyw4dtO7+wQee7C7aBg3SjmpPPx10JEnDkrvxzA8/wDffaNUi6UyYoA1h7rzzgPrI7E+tWnDaaZbcAV2ua9AgnbFqrYA9YcndeCZWXjjjjGDj8NyuXXDvvXDccdokzEOxuvvGjZ7uNpruuEN/1o8+GnQkScGSu/HMRx9BzZp7FqRIGmPH6jCg+++v9KSl8sTq7h995OluoykzU0fLvPCCrbHqAUvuxjMffaRlhho1go7EQ0VFcN99uhjsZZd5vvsOHfQP4pQpnu86moYM0bkENnKmyiy5G09s2qQnt0lXknnpJViwAIYO1enyHqtZU1usTJhg61YA0Lq19ut56inYvDnoaCLNkrvxxKxZWl5IquS+fj3cc4++qTgWv66s3FxtV/PNN74dIlruuUd/9nb2XiWW3I0npk3TcnSsV3lS+NOfdDbqU0/52ignN1fvJ07c/3Yp47TT9Oz9wQet9l4FltyNJyZPhs6doU6doCPxyLx5Oub6hht00VMfZWbqWrMTJvh6mGh54AEtyzz4YNCRRJYld1Nlq1bpsnrdugUdiUd27tSFJA47DP7854QcMjdX+8xs2ZKQw4XfSSfBFVdoaaagIOhoIsmSu6myyZP1PmmS+0MP6eDz556Dww9PyCF79tS/KZNsDbM9/vQnvZAzZEjQkUSSJXdTZZMnay+ZrKygI/HAp59qUunXDy6+OGGH7dwZjjwSXnstYYcMv8xM7eMzfDhMnx50NJFjyd1Uyc6dOkb7/POTYHGO9et1VMyRR+oyegmUnq6HnjDBRgDu5e67dUmvwYO1K6eJmyV3UyXvv685sU+foCOpoqIiGDhQ+5qMGpWwckxJl14K27bBW28l/NDhVasWPP64XuC2oZEHxJK7qZLXX9cRMuedF3QkVeCcLr7x1lvw2GN71sBLsI4dtV3yiBGBHD68+vSBHj3g97+Hr78OOprIsORuKq2wUNuuXHABZGQEHU0VPPKIjmW/7Tb9+B+QtDQdIDJpkjVG3IuI9pvJyNDOkUVFQUcUCZbcTaVNeWcXP/0E/c7bAMuX6zTL1athw4bo1Ef/+le9aHfppZrkA3bttZq7hg0LOpKQOeooLcvMmgUPPxx0NJEgLqCGFjk5OS4vLy+QY5sDtHatNo758ksd0L50KXz7LX1XPsp0dwYFNKEGO/d+TVoaNGyodYbYLJ2TTtLbCSd41hO90oqKdIjdI49oQ7BXXoFqcS0p7Ltzz9VWBEuX+tLOJrqc0z/Cb7wB770Hv/pV0BEFQkQ+dc5VuCSOJXezr3Xr9Epp7LZgwZ7vNWoELVqwpkErGr/xFDd3/oJH+8/RDlgi8Msvelu/XmsL332nWWrp0j2dsQ4+WOvaHTvqrUMHXYUnUdasgV//Gt59V5d3e/LJUGXR0aN1JOb48Z4t+pQ8Nm3S1dc3bdLFU5JyPcf9s+RuDkxBAYwbp0X0Dz/UM9vatXUAdpcuunbeySdDgwaAzgofMkQHMZx0Uhz7/+UXWLRIz/5nz4aZM2H+fE341appE/gzz9SzsU6d4JBDvH+PzulA8ptv1uTw7LNwzTXeH6eKdu7U9aKPOgr++98kGGLqtfnztf9MmzYwdaqeWKSQeJM7zrlAbu3bt3cmYGvXOvfMM86ddppzmvqca9nSuXvuce6//3Vux44yX7Z9u3NHHeVc165VPP6GDc5NmuTc3Xc716mTc9WrawwizrVt69zvfufcmDHOrVlTteMUFTn3zjvOnX667v+UU5ybO7eKwfvrmWc01A8+CDqSkBo9Wv+f9OnjXGFh0NEkFJDn4sixltxTzc6dzk2c6Fy/fs7VqKH/BbKznfvf/3Vu4cK4dvHKK/qyiRM9jm3rVufee8+5P/7RubPPdq5mzT1/dFq1cu7GG50bMcK5ggLndu3a/7527HBuxgzn7r3XuWOP1X00bercCy9EIhls3ercEUfoH9CK3mrKeuIJ/Xf9zW9S6ocUb3K3skyqWLQIXn4ZXn1VO30dfrjWnQcNgrZt497Nzp16bbRmTa2w+Foy2LED8vK0TDR9ui7SGuusVaeOXqht0kS/rlkTtm/Xckt+PixZoiN2RLSsdPXVcMklcNBBPgbsrSefhFtv1eH3PXsGHU1I3XWX9gK68UZ45pngL9QngNXcjQ5JfO01TeqzZ+tFw9xcTeg9e1ZqPbznnoObboK339Z5JQlVWKijdmbO1GGXy5frH6qtW/V20EFQty40barL4rVrp8scHXZYggP1xs6d2q9n1y4tMyfV8oVecU5bFDz4oM4w/sc/PF/nNmys5p6qCgudmzzZuQEDnMvI0I+trVs798gjzq1eXaVdr1njXP36zv3qVyn1KThQEyfqP+H99wcdSYjt2uXc0KH6g+rSperXaEKOOMsy4RjYa6pu8WIdq/2vf+nww0MP1ZEggwbpSBQP6ie33QYbN+ogExvBkRjdu8Pll8Nf/qItlQPqjBBuInDvvdCsmc4CO/VU7eFw2mlBRxao5C9QJbPvvtNp8x06aBnioYd0eNjo0TpT9OmndQijB5n41Vfh3//WT8CtW3sQu4nb00/rpYX+/eGHH4KOJsQuv3zPMN5OneD++/W6TYqymnvUrFoFY8ZoAp8xQ2uOWVlw5ZV6gdSHSR0ffwxnnaUnQu++G5qJnCklL0+nAbRsCR98kETLGfph40a45Rb9FNu8OTz6qF5jSpKPm/HW3O3MPeyKirSfxr33anmlcWP9j7thAwwdqqNg5s6FO+7wJbF/+qmWBho1gpEjLbEHJSdHr41//jl07Qo//RR0RCF28MFaopw4UUfP9OqlF9anTNkzSzoF2Jl72BQW6hjDGTPgo490+v+6dTrSpUMHzbR9+ugpnM9GjdKSfYMG+mk3M9P3Q5oKjBsHAwbo7NXhw60GX6GdO/Ui0YMPaqmyTRuty/fvH0jPfi94OloG6AYsBpYAQ8r4/kHAa8Xf/xjIrGifNlrG6ciWhQt1Ys5ddzl3zjnO1amzZ+JOZqZzAwc699przq1bl7CwVq7UOU7gXIcOzn3/fcIObeIwc6ZzxxzjXHq6czff7NyPPwYdUQRs2+bciy/qhD3Q2dA9ezr37LPOLVsWdHQHBK8mMYlIOvA1cC5QAMwBBjjnFpTY5iYg2zl3o4j0B/o45y7d335T5sx92zb4/ntYuVIn1ixdqvdLlsDChdpzBXRsbuvW2kjrjDO0p0uTJgkLs7BQPyT8+9860CAtDf7wB+2Ga+Orw2fjRp2/889/6n+dfv20F3znzhHvrZ8IX36pIwTGjNEJbwDHHKMNyU45RedHtGihv38hnBTl2SQmEekA/NE5d37x47sBnHN/LbHN5OJtZolINeB7oIHbz85DndyLinS24/5umzbpb9iGDXofu23YoMn8++/1Y+D69Xvvu1o1rW8cf7yOcGnbVhtytWzpaxbdtQt+/lnX59y4UfuE5efrXKC8PL1ounmzzgEaOFBL+Mcc41s4xiOLFulM1uHD9d+vZk0t1WRl6blC06a6JGzDhvpvW6tWqBpgBss57a38zjvaoW3OHJ0YF5ORAccdp2u4Nmy453bEEVCvnl7VrltX7+vU0UZ71avvufn0h8HL5N4X6Oacu7b48RXAac65wSW2mV+8TUHx46XF25R72aeyyf2lq6bzyH+OAsA5Ady+XyPgYo8o8X39buz7UNY2Jfax+xVS7nN7fZ2WhpM0/e1JS8cV35OehktLh2rV9D62vxIHr+zX8WxbVKSJvax/6mrVtKtjhw7aR7x7dzvzi6Kff9ZPXu++q9ffFyzQ58pSo4b+EcjI0Pyzv5tI0gwyiU9RIWzbrkMoY7fCnfrRtrCIvTNGHESA4h+isPuHef8NP3Dp02dUKsR4k3s8Yx/K+qct/Q7j2QYRuR64HuDoo4+O49D7qn/UQZzUcM3uQ8Z+diBlfi27NyqxrUjZzyNIGruTs6QXJ+r0NP06LR2ppvfUqI7UqKG/KTVqINXS997f7ve879cVff9Av47n+3Xr6q1ePb1v3FjnfBx1lI2ASQa1a+tov1gPml274NtvdeTsDz/oB8ktW7QKGOvWsG2b/sHftWv/t9RSrfhW1voCDnbshO3bNNnvLCxO+iVubhfsKv6hlvV18RW1Q4+pl5B3UpECoGmJx02AVeVsU1BcljkYWFd6R865F4AXQM/cKxNwrwdOo9cDlXmlMakjLU3LalZa85IANYpv4RdPUWgO0FxEmolIDaA/ML7UNuOBgcVf9wXe21+93RhjjL8qPHN3zhWKyGBgMpAOvOSc+0pEhqJDcsYDLwKvisgS9Iy9v59BG2OM2b+4qq3OuYnAxFLP3Vfi623AJd6GZowxprLCN4jTGGNMlVlyN8aYJGTJ3RhjkpAld2OMSUKW3I0xJgkF1vJXRNYAKyr58vpAsnS0tvcSPsnyPsDeS1hV5b0c45xrUNFGgSX3qhCRvHh6K0SBvZfwSZb3AfZewioR78XKMsYYk4QsuRtjTBKKanJ/IegAPGTvJXyS5X2AvZew8v29RLLmbowxZv+ieuZujDFmPyKb3EXkzyIyV0S+EJEpInJU0DFVlog8LCKLit/PWBE5JOiYKkNELhGRr0Rkl4hEclSDiHQTkcUiskREhgQdT2WJyEsi8mPxKmmRJSJNReR9EVlY/H/rd0HHVFkikiEin4jIl8Xv5U++Hi+qZRkRqeec21T89S1AK+fcjQGHVSkich7aA79QRB4EcM7dFXBYB0xEWgK7gOeBO5xzIV0kt2zxLAYfFSLyK2AL8C/n3ElBx1NZItIIaOSc+0xE6gKfAhdG9N9EgNrOuS0iUh2YAfzOOTfbj+NF9sw9ltiL1eaAFzcMD+fcFOdcYfHD2ehqV5HjnFvonFscdBxVcCqwxDm3zDm3AxgJ9A44pkpxzk2njNXQosY5t9o591nx15uBhUDjYKOqHKe2FD+sXnzzLW9FNrkDiMgDIvIt8Gvgvoq2j4irgUlBB5GiGgPflnhcQEQTSTISkUygLfBxsJFUnoiki8gXwI/Au845395LqJO7iEwVkfll3HoDOOd+75xrCgwHBgcb7f5V9F6Kt/k9UIi+n1CK531EWFwLvZvEE5E6wBjg1lKf2iPFOVfknGuDfjo/VUR8K5mFet1759w5cW76H2ACcL+P4VRJRe9FRAYCPYGuYV5/9gD+TaIonsXgTYIV16fHAMOdc28EHY8XnHMbROQDoBvgy0XvUJ+574+INC/xsBewKKhYqkpEugF3Ab2cc1uDjieFxbMYvEmg4ouQLwILnXOPBR1PVYhIg9hIOBGpCZyDj3kryqNlxgAnoKMzVgA3Oue+CzaqyileWPwgYG3xU7OjOPJHRPoAfwMaABuAL5xz5wcb1YERkVzgCfYsBv9AwCFVioiMALqg3Qd/AO53zr0YaFCVICKdgY+AeejvOsA9xes6R4qIZAOvoP+30oBRzrmhvh0vqsndGGNM+SJbljHGGFM+S+7GGJOELLkbY0wSsuRujDFJyJK7McYkIUvuxhiThCy5G2NMErLkbowxSej/A2QntWHSphnaAAAAAElFTkSuQmCC\n", "text/plain": [ - "0.43448275862068964" + "
" ] }, - "execution_count": 17, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - ".9*.0035/(.1*.041+.9*.0035)" + "x=np.linspace(-3,3,1000)\n", + "y1=lambda1*norm.pdf(x,-1,.5)+lambda2*norm.pdf(x,1,.5)\n", + "y2=wt1*norm.pdf(x,post_mean1,np.sqrt(post_var1))+wt2*norm.pdf(x,post_mean2,np.sqrt(post_var2))\n", + "fig,ax=plt.subplots(1)\n", + "ax.plot(x,y1,color='red',label='prior')\n", + "ax.plot(x,y2,color='blue',label='posterior')\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This isn't consistent with the solutions, which I don't understand. The key issue is the meaning of the statement that \"the variance of each observation is known to be 1\". How does one compute $p_{m}(\\{y_{i}\\})$?" ] }, { @@ -136,7 +173,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.6.5" } }, "nbformat": 4, diff --git a/Useful Formulae.ipynb b/Useful Formulae.ipynb index 3fc65f4..5a3d891 100644 --- a/Useful Formulae.ipynb +++ b/Useful Formulae.ipynb @@ -136,7 +136,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.6.5" } }, "nbformat": 4,