Skip to content

Instantly share code, notes, and snippets.

@ricklupton

ricklupton/Uncertainty.ipynb

Last active November 14, 2017 15:44
Show Gist options
  • Save ricklupton/2a90b0a92241bec768e7f469ce832292 to your computer and use it in GitHub Desktop.
Save ricklupton/2a90b0a92241bec768e7f469ce832292 to your computer and use it in GitHub Desktop.
Download ZIP
Dirichlet allocation uncertainty
Display the source blob
Display the rendered blob
Raw
Uncertainty.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"from numpy.random import dirichlet\n",
"import pymc3 as pm\n",
"\n",
"def dirichlet_params(shares, i, stddev, clip_stddev=False):\n",
" # normalise\n",
" factor = sum(shares)\n",
" shares = np.array(shares) / factor\n",
" stddev /= factor\n",
" \n",
" # check stddev is within range\n",
" mi = shares[i]\n",
" limit = np.sqrt(mi * (1 - mi) / (1 + len(shares)))\n",
" if stddev > limit:\n",
" if clip_stddev:\n",
" print('Warning: clipped std dev from %.2g to %.2g' % (stddev * factor, limit * factor))\n",
" stddev = limit\n",
" else:\n",
" raise ValueError('stddev is too high (%.2g > %.2g)' % (stddev * factor, limit * factor))\n",
" concentration = mi * (1 - mi) / stddev**2 - 1\n",
" if not np.isfinite(concentration):\n",
" concentration = 1e10\n",
" \n",
" if clip_stddev:\n",
" return (concentration * shares, stddev * factor)\n",
" else:\n",
" return concentration * shares"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Describing the uncertainty of an allocation is, perhaps surprisingly, rather difficult conceptually. For example, assume we have 100 MJ of energy split into three parts of 80%, 17% and 3%. If the \"uncertainty\" is \"10%\", what could that mean? Intuitively, I'd expect the biggest share to be about 80 MJ and vary by about 10%, i.e. from 72 to 88 MJ. Similarly, I'd expect the middle share to vary from about 15 to 19 MJ and the smallest share from about 2.7 to 3.3 MJ. \n",
"\n",
"However, this is not possible. To see why, consider that the smallest share reduces to the minimum expected, from 3 MJ to 2.7 MJ. The other two shares must increase by a total of 0.3 MJ. Even if all of this increase went to the 80 MJ share, that would still be only a variation of 0.3 / 80 = 0.4%: much less than the expected variation. Conversely, if the largest share were to increase by 10% the corresponding change in the smaller shares would be bigger than expected.\n",
"\n",
"Because it is not possible to have \"10% uncertainty\" on all parts of the allocation, we must make a choice about how to interpret the uncertainty. At the same time, we have no information about the specific probability distributions describing the allocations, or their covariance (i.e. which part would increase if another part descreases). In this paper we choose to use the Dirichlet distribution to describe uncertain allocations because it naturally represents allocations that add up to 100%, and because it has a simple parameterisation in terms of the mean shares and a \"concentration parameter\" determining the level of uncertainty. We consider the following rules for the interpretation of a \"10% uncertainty\":\n",
"- 10 percentage points (pp) on the biggest (80%, case **1a**) or medium (17%, case **2a**) share, but not on the smallest (3%) as this would give a negative percentage.\n",
"- 10% of the percentage, i.e. 8pp on 80% (**1b**), 1.7% on 17% (**2b**) or 0.3% on 3% (**3b**).\n",
"- some kind of \"average\" variation across all 3 shares? not sure how to do this\n",
"\n",
"The figure & table below compare the results of these rules to the expected outcome -- 10% variations in the size of the output parts -- above."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"shares = [0.8, 0.17, 0.03]\n",
"\n",
"expected = np.array([\n",
" [80 * 0.9, 80 * 1.1],\n",
" [17 * 0.9, 17 * 1.1],\n",
" [3 * 0.9, 3 * 1.1],\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"inp = 100.0\n",
"N = 10000\n",
"# The 0.5 is to convert uncertainty ranges (95% confidence) to std-dev\n",
"out1a = inp * dirichlet(dirichlet_params(shares, 0, 0.5 * 0.10), N)\n",
"out2a = inp * dirichlet(dirichlet_params(shares, 1, 0.5 * 0.10), N)\n",
"out1b = inp * dirichlet(dirichlet_params(shares, 0, 0.5 * 0.10 * 0.80), N)\n",
"out2b = inp * dirichlet(dirichlet_params(shares, 1, 0.5 * 0.10 * 0.17), N)\n",
"out3b = inp * dirichlet(dirichlet_params(shares, 1, 0.5 * 0.10 * 0.03), N)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"labels = ['1a', '2a', '1b', '2b', '3b']\n",
"data = [out1a, out2a, out1b, out2b, out3b]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAJCCAYAAABdxXsqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGWxJREFUeJzt3X+MZWV9x/HPd2DH6TIo0CFkuyDXVsfGoAtk/RVMKram\nakxbkm4EDbVK6h+2qTYkVm1MY+I/TRqrDYaUFJbaCLQX2VZIp4ZajLGYtTMru4ssjphlV5a17AYF\nR0JXOt/+cc+wd8f5ce/cc+Y5z/N9v5LJzL1358z3ybn57Pc5zznnmrsLACIYS10AAGwWAg9AGAQe\ngDAIPABhEHgAwiDwAIRB4AEIg8ADEAaBByCMs5v+A1NTU97pdGrZ1qlTp7S4uKixsTGNj4/Xss1R\n6lhL6hqRl0HeU23Qpvf13NzcSXe/cJjfaTzwOp2OZmdna9nW/Py8JicntbCwoOnp6Vq2OUoda0ld\nI/IyyHuqDdr0vjazI8P+DlNaAGEQeADCIPAAhNH6wOvOd1OXAKAQrQ88AKgLgQcgDAIPQBgEHoAw\nCDwAYTR+pUUdllZqd2hH4koA5IwOD0AYBB6AMFo7peWEYwB1o8MDEAaBByAMAg9AGAQegDAIPABh\nEHgAwiDwAIRB4AEIg8ADEEZrr7RYyczhGU1sndDVF12duhQAGaLDAxAGgQcgDAIPQBgEHoAwCDwA\nYRB4AMLI6rSUpGZ3n/75pVelqwPAhtHhAQiDwAMQBoEHIAwCD0AYLFqMav9dp3/ecW26OgCsi8Db\niIN3S1snUlcBDGzPvmMv/nzNldsTVpLWuoFnZrdJerekp9z9suZLylh/t/dr705XB8Ii2NY2SId3\nu6SbJH2x2VIGN3N4RtPT083+kf7z7gAUYd3Ac/dvmFmn/zkz+yNJH5I0LukxSde7+3NNFAhgY/q7\nPfRsdJX2Hnd/vbvvkHRI0g011gQAjdjoosVlZvYZSedJmpT01fpKAoBmbLTDu13Sn7j7ayV9WhJL\nlgBab6Md3rmSjpvZFknvk8TBAiATkVdy1+3wzOxOSd+S9Goze8LMbpD0KUl7Jf2XpEebLREA6jHI\nKu11q7x0c821ANhk0bo9rrQAMsfpJ4Mj8Jpy8G7p2W2nH+/8QLpaAEhqYeB157upSwBQqNYFXlJc\nTgYUjfvhAQiDwAMQRraB153vcrwPwFA4hgdkhtNQNi7bDg8AhkWHt1n6V4A5Jw8ttLxzLPHKCzo8\nAGEQeADCIPAAhJF94HF6CoBBZR94ADAoVmmBDHDuXT3o8ACEQYfHHVKAMOjwAIRB4AEIo5gpbXe+\nq13Tu1KXMRguM8MAWKioXzGBB6BeJX6iGVNaAGEQeADCIPAAhBHzGB7n3gEhxQw8oKVYmW0WgZca\np6ggA0tB/PxzP9ON09OJq9k4juEBCKOowOPeeADWwpS2TZjeAo0qqsMDgLXQ4QGJ3bH3qI4feVIT\nW89JXUrxiuzwOJYHYCVFBh4ArIQpbVuxgFG0O/YeTV1CSHECj8vJgFr0h/V73/jyhJUMjyktgDCK\n7vCyugvyWpZ3p0xxgQ0pOvCANuG4XXrFB97S6SkF9HlA6+R2PK/swHv8wdM/X/C6dHUAaIWyA69U\nnLICbEiYwOs+fUCStItOD5uI43btEu60lKXgK8bs7tNfANZUVofXf8xuDctDr5iub7XQY9rbKLq4\nnhwWMMoKvA1aresrMggJvw0j2AbX1vAzd2/2D5idkHRkiF+ZknSyoXJSYUx5KG1MpY1HOnNMl7r7\nhcP8cuOBNywzm3X3nanrqBNjykNpYyptPNLoYwq3aAEgLgIPQBhtDLxbUhfQAMaUh9LGVNp4pBHH\n1LpjeADQlDZ2eADQCAIPQBgEHoAwCDwAYRB4AMIg8ACEQeABCIPAAxAGgQcgDAIPQBgEHoAwCDwA\nYRB4AMJo/DMtpqamvNPp1LKtU6dOaXFxUWNjYxofH69lm6PUsZbUNSIvg7yn2qBN7+u5ubmTw97i\nvfHA63Q6mp2drWVb8/Pzmpyc1MLCgqanp2vZ5ih1rCV1jcjLIO+pNmjT+9rMhvmsHElMaQEEQuAB\nCKO1gded76YuAUBhWht4AFA3Ag9AGAQegDAIPABhEHgAwmh14HXnu6zWAqhNqwMPAOpE4AEIg8AD\nEEbjNw+ow9JxvB3akbgSADmjwwMQBoEHIAwCD0AYBB6AMAg8AGEQeADCIPAAhEHgAQiDwAMQBoEH\nIAwCD0AYWQXezOGZ1CUAyFhWgQcAoyDwAIRB4AEII4v74QEYzZ59x854fM2V2xNVkhYdHoAw6PA2\nav9dp3/ecW26OgAMjMAbxuzu3vcjx6WtE2lrAdaxfBoLAg8IqT8MIx3PI/DWstTRARmgo1tfdoF3\n7w/u1fPPPa/p6enUpQDIzLqBZ2aXSPqipIskuaRb3P3zTReWFRYwgCwM0uG9IOlGd99nZudKmjOz\n+939kSYKWvoMWgCo27qB5+7HJR2vfv6pmR2StN3MrpL0IUnjkh6TdL27P9dksQAwiqFOPDazjqQr\nJO2VdI+7v97dd0g6JOmG2qsDgBoNvGhhZpOSvizpo+7+rJn9hpl9RtJ5kiYlfbWhGgGsgFXZ4Q0U\neGa2Rb2w+5K731M9fbuk33P3/Wb2h5Le2kSBAJoV6TrbQVZpTdKtkg65+2f7XjpX0vEqDN8nif9u\npNMrts89L01/Mm0tAM4wSId3laTrJR00s4eq5z4p6VPqHcs7UX0/t5EKNxsnGwPFGmSV9puSbJWX\nb663HABoDreHAhAGgQcgjOyupQUi41SU0dDhAQgj28Drzne57hbAULINPAAYFsfwmrT8nL6dH0hT\nBwBJdHgAAqHD20z9HR/dHrDp6PAAhEGHB7QY593VK/vA6853tWt612gb4YYBwItKvl0UU1oAYRB4\nAMIg8ACEQeABCIPAAxAGgQcgjOxPSwHQrP7TVN4+/bKElYyOwEuFy8ywCk42bg5TWgBhFBF43AwU\nwCDiTmm5nAwIp4gODwAGQeABCCPulBZoiTv2HtXxI09qYus5qUspHoHXBpyiAmyKoqa0rNQCWEtR\ngQcAa4k1peVUFGAk9z70pLb9eOLFx+9948sTVjO8WIGXA47nAY1hSgsgjOICj8vMAKymuMADgNVw\nDA9I4I69R1OXEFKxgbc0rR35M2tTWr6qzCIGMBKmtADCKLbDkyQ9/qD07ELqKgCmsC1BhwcgjLI7\nPEndpw9o1wWvS11GPTgpGS3T37nmcNUFHR6AMEIEXvfpA+o+fSB1GQASK35KWyymt63GIkU7lRd4\njz+46ktLXV4xx/TQKtFDbvn423hML8SUdrniprezu09/AVhVGR3eGl3daort9pjqbproHV2Oygi8\nEazU7RUTglyaVjtCbnBtnOKauzf7B8xOSDoyxK9MSTrZUDmpMKY8lDam0sYjnTmmS939wmF+ufHA\nG5aZzbr7ztR11Ikx5aG0MZU2Hmn0MYVctAAQE4EHIIw2Bt4tqQtoAGPKQ2ljKm080ohjat0xPABo\nShs7PABoBIEHIAwCD0AYBB6AMAg8AGEQeADCIPAAhEHgAQiDwAMQBoEHIAwCD0AYBB6AMAg8AGE0\n/pkWU1NT3ul0atnWqVOntLi4qLGxMY2Pj9eyzVHqWEvqGpGXQd5TbdCm9/Xc3NzJYW/x3njgdTod\nzc7O1rKt+fl5TU5OamFhQdPT07Vsc5Q61pK6RuRlkPdUG7TpfW1mw3xWjiSmtAACIfAAhEHgAQgj\ni8DrzndTlwCgAFkEHgDUgcADEAaBByAMAg9AGK0PPBYsANSl9YEHAHUh8ACEkU3gMbUFMKpsAg8A\nRkXgAUHs2XdMe/YdS11GUgQegDAavx9ekQ7eLW2dkHZcm7oS4BdE7+LWQuANanZ36goAjKjVgbd8\nZXbm8Iwmtk7o6ouuTlTRMvvvWvl5Oj+glTiGByAMAg9AGAQegDAIPABhtHrRAkD9+k9buebK7Qkr\n2XwEHlAAzr0bDIG3Fs69A4rCMTwAYWQZeDOHZ1KXsLb9d/UuPwPQKlkGHgBsBIEHIAwWLYBMsTI7\nPDo8AGEQeADCIPAAhEHgAQiDwAMQxrqBZ2a3mdlTZvZw33NfN7OdzZYGAPUa5LSU2yXdJOmLzZYC\nYLNFu3PKuh2eu39D0tMrvHS9mT1kZg+b2RvqLw0A6jXKMbyt7n65pA9Luq2megCgMaNcaXGn1OsA\nzeylZnaeu/+kprrK0H97qZ0fSFcHAEmjdXi+zuNGLf8IRwBYzyiB9x5JMrO3SHrG3Z+ppyQAaMa6\nU1ozu1PSWyVNmdkTkv6yeul5M/uOpC2SPthYhQBQk3UDz92vW+HpWxuopT24tTtQJK60ABAG98MD\nMsI98EbT2g6PVVgAdWtt4AFA3Qg8AGFkHXhMe4H67Nl37MWvUmUdeAAwDAIPQBgEHoAwCLzNMrub\nKziAxAg8AGEQeADCIPAAhMG1tEs4voYWK/ncuM1EhwcgjOwDj6stAAyKKe1m44N9gGSy7/AAYFAE\nHoAwmNIC+AX9q8LXXLk9YSX1osMDEAaBByAMprRAS3Gycf2K6PA4Fw/AIIoIPAAYBIEHIAyO4aXE\nVRfApoodeNwhBQilmCltd77L4gWANRUTeACwnthTWqBlOPeuWXR4AMIg8NqCj3EEGkfgAVjTnn3H\niplqt+4Y3qgrrd35rnZN76qpGqB5d+w9quNHntTE1nNSl1I8OjwAYbSuw9sUHCsDQqLDAxBGzA6v\nzbi+FmhMkYHHwgXa7o69R1OXEFKRgVcMuj20yJ59x/T8cz/Tth9P6L1vfHnqcjak2GN43EgAbXTH\n3qN0dwkVG3gSoQfgTHGmtJyKAtSmv0vNaXpbdIcnFdTlca0tMLKyO7zHH5Qkdavvuy54Xcpq6sFC\nBrBhZQfeCrpPHygj+NBqkRYmcprelhl4VUe3XPfpA2d8lzLv+laa4tL1JRMp5FbT9vAzd2/2D5id\nkHRkiF+ZknSyoXJSYUx5KG1MpY1HOnNMl7r7hcP8cuOBNywzm3X3nanrqBNjykNpYyptPNLoYyp+\nlRYAlhB4AMJoY+DdkrqABjCmPJQ2ptLGI404ptYdwwOAprSxwwOARhB4AMIg8ACEQeABCIPAAxAG\ngQcgDAIPQBgEHoAwCDwAYRB4AMIg8ACEQeABCIPAAxBG459pMTU15Z1Op5ZtnTp1SouLixobG9P4\n+Hgt2xyljrWkrhF5GeQ91QZtel/Pzc2dHPYW740HXqfT0ezsbC3bmp+f1+TkpBYWFjQ9PV3LNkep\nYy2pa0ReBnlPtUGb3tdmNsxn5UhiSgsgEAIPQBgEHoAwCDwAYRB4o9p/V+8LaJk9+46lLqF1CDwA\nYRB4AMIg8ICCMa09E4EHIAwCD0AYBN4oWJ1FSzGVXRmBByCMxm8eUKeZwzOa2Dqhqy+6OnUpADJE\nhwcgDAIPQBgEHoAwCDwAYRB4AMIg8ACEQeABCIPAqwtXXQCtR+ABCIPAAxAGgQcgDAIPQBgEHoAw\nCDwAYRB4QGG4+efqCDwAYRB4AMIg8ACEQeABCIPAAxAGgQcUjlXb0wg8AGGsG3hmdomZPWBmj5jZ\nd83sI9XzXzeznc2XCAD1GORzaV+QdKO77zOzcyXNmdn9DdcFALVbt8Nz9+Puvq/6+aeSDknaXr18\nvZk9ZGYPm9kbGqwTAEY21DE8M+tIukLS3uqpre5+uaQPS7qt1soAoGYDB56ZTUr6sqSPuvuz1dN3\nSpK7f0PSS83svPpLBIB6DBR4ZrZFvbD7krvf0/eSL/unyx8DQGsMskprkm6VdMjdP7vs5fdU/+Yt\nkp5x92fqL7Gl+NAeIDuDrNJeJel6SQfN7KHquU9W3583s+9I2iLpgw3UBwC1WTfw3P2bkmyFl/6t\n/nIAoDlcaQEgDAKvThzXA1qNwAMQBoEHIAwCD0AYBB6AMAg8AGEQeEBBuLvx2gg8AGEQeADCIPAA\nhEHgAQiDwAMQBoEHIIwsA2/m8EzqEoCscLpKT5aBBwAbQeABCIPAAxAGgQcgjGwCrzvfTV0CgMxl\nE3gAMCoCDygEp56sj8ADEAaBtxEH7179NT65DGgtAg9AGAQeEATH+Ag8AIEQeADCIPAAhEHgAQiD\nwAMQBoEHIAwCrwmcfAy0EoEHIAwCDygAJxUPhsADEAaBByAMAg9AGATesGZ3p64AwAYReADCyDbw\n+FAfAMPKNvBab627IgOJRD99hcAbBsfvgKwReEDmondtwyDwAISRReCxQAGsbCPdXeSOMIvAW82m\nBuFGjt9xzA8tFTX0zN2b/QNmJyQdGeJXpiSdbKicVBhTHkobU2njkc4c06XufuEwv9x44A3LzGbd\nfWfqOurEmPJQ2phKG480+piyntICwDAIPABhtDHwbkldQAMYUx5KG1Np45FGHFPrjuEBQFPa2OEB\nQCMIPABhEHgAwiDwAIRB4AEIg8ADEAaBByAMAg9AGAQegDAIPABhEHgAwiDwAIRB4AEI4+ym/8DU\n1JR3Op1atnXq1CktLi5qbGxM4+PjtWxzlDrWkrpG5GWQ91QbtOl9PTc3d3LYW7w3HnidTkezs7O1\nbGt+fl6Tk5NaWFjQ9PR0LdscpY61pK4ReRnkPdUGbXpfm9kwn5UjiSktgEAIPABhEHgAwiDwAIRB\n4AEIg8ADEAaBByAMAg9AGAQegDAIPABhEHgAwiDwAIRB4AEIg8ADEAaBByAMAg9AGAQegDAIPABh\nEHgAwiDwgELt2XcsdQmtQ+ABCIPAAxAGgQcgDAIPQBgEHoAwCDwAYRB4AMIg8ACEQeABCIPAAxAG\ngQcgDAIPQBgEHoAwCDwAYRB4AMJYN/DMbMLMvm1m+83su2b26er5x81sqvkSAaAeZw/wb/5X0tvc\nfcHMtkj6ppnNNFwXANRu3Q7Pexaqh1uqL68ef8zMDlYd4CubKhIA6jDQMTwzO8vMHpL0lKT73X1v\n9dIz7v5aSTdJ+lxDNQJALQYKPHf/P3e/XNLFkt5gZpdVL93Z9/3NDdQHALUZapXW3X8i6QFJ71h6\nqv/luooCgCYMskp7oZmdV/38S5LeLunR6uX39H3/ViMVAkBNBlml3SbpH8zsLPUC8p/d/T4zu0nS\n+WZ2QL2V3OsarBMARrZu4Ln7AUlXrPB8p/rxz2uuCQAawZUWAMIg8ACEQeABCIPAAxAGgQcgDAIP\nQBgEHoAwCDwAYRB4ddh/V+oKAAyAwAMQBoEHIAwCD0AYBB6AMAg8AGEQeADCIPAAhEHgAQiDwAMQ\nBoEHIAwCD0AYBB6AMAi8UXHjACAbBB6AMAg8AGEQeADCIPAAhEHgAQiDwAMQBoEHIAwCD0AYBB6A\nMAg8AGEQeADCIPAAhEHgAQiDwAMQBoEHIAwCDyjQnn3HUpfQSgQegDAIPABhEHgAwiDwAISRXeDd\n+4N7U5cAIFPZBR4AbBSBByAMAg9AGAQeUDBOQD5TVoE3c3gmdQkAMmbu3uwfMDsh6cgQvzIl6WRD\n5aTCmPJQ2phKG4905pgudfcLh/nlxgNvWGY26+47U9dRJ8aUh9LGVNp4pNHHlNWUFgBGQeABCKON\ngXdL6gIawJjyUNqYShuPNOKYWncMDwCa0sYODwAa0arAM7N3mNn3zOwxM/t46nqGZWaXmNkDZvaI\nmX3XzD5SPX+Bmd1vZt+vvp+futZhmdlZZvYdM7uvevwKM9tb7at/MrPx1DUOw8zOM7O7zexRMztk\nZm/OfT+Z2Z9V77uHzexOM5vIbT+Z2W1m9pSZPdz33Ir7xXr+thrbATO7cr3ttybwzOwsSV+Q9E5J\nr5F0nZm9Jm1VQ3tB0o3u/hpJb5L0x9UYPi7pa+7+Kklfqx7n5iOSDvU9/itJf+Pur5T0Y0k3JKlq\n4z4v6d/d/dcl7VBvbNnuJzPbLulPJe1098sknSXpWuW3n26X9I5lz622X94p6VXV14ck3bzu1t29\nFV+S3izpq32PPyHpE6nrGnFM/yrp7ZK+J2lb9dw2Sd9LXduQ47i4eqO9TdJ9kky9kz/PXmnftf1L\n0sskHVZ1DLvv+Wz3k6Ttkn4o6QJJZ1f76bdz3E+SOpIeXm+/SPo7Sdet9O9W+2pNh6fTO2zJE9Vz\nWTKzjqQrJO2VdJG7H69e+pGkixKVtVGfk/QxSYvV41+W9BN3f6F6nNu+eoWkE5J2V9P0vzezc5Tx\nfnL3Y5L+WtJRScclPSNpTnnvpyWr7ZehM6NNgVcMM5uU9GVJH3X3Z/tf895/RdksjZvZuyU95e5z\nqWup0dmSrpR0s7tfIelnWjZ9zXA/nS/pd9UL81+RdI5+cWqYvVH3S5sC75ikS/oeX1w9lxUz26Je\n2H3J3e+pnv4fM9tWvb5N0lOp6tuAqyT9jpk9Luku9aa1n5d0npmdXf2b3PbVE5KecPe91eO71QvA\nnPfTb0k67O4n3P3nku5Rb9/lvJ+WrLZfhs6MNgXef0t6VbWqNK7eAdevJK5pKGZmkm6VdMjdP9v3\n0lckvb/6+f3qHdvLgrt/wt0vdveOevvkP939fZIekPT71T/LbUw/kvRDM3t19dRvSnpEGe8n9aay\nbzKzrdX7cGlM2e6nPqvtl69I+oNqtfZNkp7pm/quLPUBymUHK98laV7SDyT9Rep6NlD/W9Rrtw9I\neqj6epd6x7y+Jun7kv5D0gWpa93g+N4q6b7q51+V9G1Jj0nqSnpJ6vqGHMvlkmarffUvks7PfT9J\n+rSkRyU9LOkfJb0kt/0k6U71jkH+XL1O/IbV9ot6i2dfqPLioHor1GtunystAITRpiktADSKwAMQ\nBoEHIAwCD0AYBB6AMAg8AGEQeADCIPAAhPH/1Sasup5jpYQAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fa7d47154e0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(5, 1, sharex=True, figsize=(5, 10))\n",
"for a, x, l in zip(ax, data, labels):\n",
" for i in range(3):\n",
" a.hist(x[:, i], alpha=0.4, bins=30);\n",
" #x0, x1 = pm.hpd(x[:, i])\n",
" #a.axvspan(x0, x1, alpha=0.1, color='k')\n",
" a.axvspan(expected[i, 0], expected[i, 1], alpha=0.1, color='k')\n",
" a.set_ylabel(l, rotation=0, labelpad=10)\n",
" a.set_yticks([])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr>\n",
" <th></th>\n",
" <th colspan=\"2\" halign=\"left\">0</th>\n",
" <th colspan=\"2\" halign=\"left\">1</th>\n",
" <th colspan=\"2\" halign=\"left\">2</th>\n",
" </tr>\n",
" <tr>\n",
" <th></th>\n",
" <th>min</th>\n",
" <th>max</th>\n",
" <th>min</th>\n",
" <th>max</th>\n",
" <th>min</th>\n",
" <th>max</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>1a</th>\n",
" <td>70.2</td>\n",
" <td>89.4</td>\n",
" <td>8.1</td>\n",
" <td>26.2</td>\n",
" <td>0.1</td>\n",
" <td>7.1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2a</th>\n",
" <td>69.4</td>\n",
" <td>89.7</td>\n",
" <td>7.8</td>\n",
" <td>26.9</td>\n",
" <td>0.0</td>\n",
" <td>7.4</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1b</th>\n",
" <td>71.9</td>\n",
" <td>87.6</td>\n",
" <td>9.5</td>\n",
" <td>24.3</td>\n",
" <td>0.3</td>\n",
" <td>6.3</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2b</th>\n",
" <td>78.2</td>\n",
" <td>81.7</td>\n",
" <td>15.3</td>\n",
" <td>18.7</td>\n",
" <td>2.3</td>\n",
" <td>3.8</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3b</th>\n",
" <td>79.7</td>\n",
" <td>80.3</td>\n",
" <td>16.7</td>\n",
" <td>17.3</td>\n",
" <td>2.9</td>\n",
" <td>3.1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>expec</th>\n",
" <td>72.0</td>\n",
" <td>88.0</td>\n",
" <td>15.3</td>\n",
" <td>18.7</td>\n",
" <td>2.7</td>\n",
" <td>3.3</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" 0 1 2 \n",
" min max min max min max\n",
"1a 70.2 89.4 8.1 26.2 0.1 7.1\n",
"2a 69.4 89.7 7.8 26.9 0.0 7.4\n",
"1b 71.9 87.6 9.5 24.3 0.3 6.3\n",
"2b 78.2 81.7 15.3 18.7 2.3 3.8\n",
"3b 79.7 80.3 16.7 17.3 2.9 3.1\n",
"expec 72.0 88.0 15.3 18.7 2.7 3.3"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pandas as pd\n",
"df = pd.DataFrame.from_items(\n",
" [(l, [z for i in range(3) for z in pm.hpd(d[:, i])]) \n",
" for l, d in zip(labels, data)] + [('expec', expected.flatten())],\n",
" columns=pd.MultiIndex.from_tuples([(i, z) for i in range(3) for z in ['min', 'max']]),\n",
" orient='index')\n",
"df.applymap(lambda x: '%.1f' % x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The **a** cases give too big an uncertainty range for all the parts.\n",
"\n",
"The **b** cases, unsurprisingly, give the expected uncertainty range for the *reference* part whose uncertainty was specified: **1b** has the expected range for the first part, **2b** for the second part, etc. But this example shows that parts smaller than the reference part have greater than expected uncertainty, while parts bigger than the reference part have smaller than expected uncertainty.\n",
"\n",
"Because bigger parts have a bigger influence on the overall results, it is better to exaggerate the uncertainty of small parts than it is to deny the uncertainty of large parts. We therefore use rule **1b**: adjust the distribution to match the uncertainty in the largest part to the specified value."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment