From 997683dcbc49b686b1041d377b7e5b9bb9ba128b Mon Sep 17 00:00:00 2001 From: Mauro Donega Date: Tue, 28 May 2024 17:43:11 +0200 Subject: [PATCH] Added Bayesian limits near boundaries --- book/_toc.yml | 4 +- .../LimitNearBoundariesBayesianApproach.ipynb | 218 ++++++++++++++++++ ...earBoundariesBayesianApproachPoisson.ipynb | 209 +++++++++++++++++ 3 files changed, 430 insertions(+), 1 deletion(-) create mode 100644 book/interactive-nbs/LimitNearBoundariesBayesianApproach.ipynb create mode 100644 book/interactive-nbs/LimitNearBoundariesBayesianApproachPoisson.ipynb diff --git a/book/_toc.yml b/book/_toc.yml index 7af77d6..314af09 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -28,8 +28,10 @@ sections: - file: confidenceIntervals.ipynb sections: - file: interactive-nbs/Poisson_CI.ipynb - - file: interactive-nbs/UpperLimit.ipynb - file: interactive-nbs/Gaussian_CI.ipynb + - file: interactive-nbs/LimitNearBoundariesBayesianApproach.ipynb + - file: interactive-nbs/LimitNearBoundariesBayesianApproachPoisson.ipynb + - file: interactive-nbs/UpperLimit.ipynb - file: mva.ipynb - file: appendix.ipynb sections: diff --git a/book/interactive-nbs/LimitNearBoundariesBayesianApproach.ipynb b/book/interactive-nbs/LimitNearBoundariesBayesianApproach.ipynb new file mode 100644 index 0000000..d5f37db --- /dev/null +++ b/book/interactive-nbs/LimitNearBoundariesBayesianApproach.ipynb @@ -0,0 +1,218 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Compute the bayesian upper limit for a gaussian near the physical boundary" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "theta_obs numerator denominator Ratio limit\n", + "-4.000000 0.000030 0.000032 0.950079 0.660000\n", + "-3.733333 0.000090 0.000094 0.950208 0.697000\n", + "-3.466667 0.000250 0.000263 0.950162 0.737000\n", + "-3.200000 0.000653 0.000687 0.950064 0.781000\n", + "-2.933333 0.001593 0.001677 0.950002 0.830000\n", + "-2.666667 0.003639 0.003830 0.950031 0.885000\n", + "-2.400000 0.007789 0.008198 0.950173 0.947000\n", + "-2.133333 0.015628 0.016449 0.950086 1.015000\n", + "-1.866667 0.029429 0.030974 0.950124 1.092000\n", + "-1.600000 0.052065 0.054799 0.950096 1.178000\n", + "-1.333333 0.086662 0.091211 0.950124 1.275000\n", + "-1.066667 0.135912 0.143061 0.950025 1.383000\n", + "-0.800000 0.201272 0.211855 0.950045 1.505000\n", + "-0.533333 0.282061 0.296901 0.950017 1.641000\n", + "-0.266667 0.375148 0.394863 0.950071 1.793000\n", + "0.000000 0.475002 0.500000 0.950004 1.960000\n", + "0.266667 0.574901 0.605137 0.950034 2.144000\n", + "0.533333 0.668002 0.703099 0.950083 2.344000\n", + "0.800000 0.748771 0.788145 0.950042 2.558000\n", + "1.066667 0.814162 0.856939 0.950082 2.786000\n", + "1.333333 0.863434 0.908789 0.950093 3.025000\n", + "1.600000 0.898037 0.945201 0.950102 3.273000\n", + "1.866667 0.920602 0.969026 0.950028 3.527000\n", + "2.133333 0.934454 0.983551 0.950081 3.787000\n", + "2.400000 0.942229 0.991802 0.950016 4.049000\n", + "2.666667 0.946425 0.996170 0.950064 4.314000\n", + "2.933333 0.948510 0.998323 0.950103 4.580000\n", + "3.200000 0.949431 0.999313 0.950084 4.846000\n", + "3.466667 0.949786 0.999737 0.950036 5.112000\n", + "3.733333 0.949989 0.999906 0.950079 5.379000\n", + "4.000000 0.949983 0.999968 0.950014 5.645000\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjL0lEQVR4nO3dd3hVVdrG4d+iKCUW0Ih0HBCkqCBNxQLICBawoaOO2OWzYEcRGWDUcSyM2FCBEVSUMiCONAuKiSKCCEgTCDIU6cUCBKQkWd8fLxFQSsjZJ/ucnee+rnOlnezzboNPVtZe+13Oe4+IiERHkbALEBGRYCnYRUQiRsEuIhIxCnYRkYhRsIuIRIyCXUQkYooFcRDn3FJgM5ANZHnvGwVxXBEROXSBBPsuLbz3GwI8noiI5IOmYkREIsYFceepc24J8DPggX7e+/77eE5HoCNAiRIlGlapUiXm101UOTk5FCkS3d+ZUT6/KJ8b6PyS3cKFCzd471MP9ryggr2C936Vc+444BPgbu/9F/t7fq1atXxGRkbMr5uo0tPTad68edhlxE2Uzy/K5wY6v2TnnJuel2uYgfxq896v2vV2HfBfoEkQxxURkUMXc7A750o7547IfR84H5gb63FFRCR/glgVUw74r3Mu93hDvPcfBXBcERHJh5iD3Xu/GDg1gFpERCQA0b18LCJSSCnYRUQiRsEuIhIxCnYRkYhRsIuIRIyCXUQkYhTsIiIRo2AXEYkYBbuISMQo2EVEIkbBLiISMQp2EZGIUbCLiESMgl1EJGIU7CIiEaNgFxGJGAW7iEjEKNhFRCJGwS4iEjEKdhGRiFGwi4hEjIJdRCRiFOwiIhGjYBcRiRgFu4hIxCjYRUQiRsEuIhIxgQW7c66oc+5b59zYoI4pIiKHLsgR+73A/ACPJyIi+RBIsDvnKgEXAa8HcTwREcm/oEbsLwAPAzkBHU9ERPLJee9jO4BzFwMXeu/vdM41Bzp77y/ex/M6Ah0BUlNTGw4fPjym101kmZmZpKSkhF1G3ET5/KJ8bqDzS3YtWrSY7r1vdLDnBRHsTwEdgCygBHAk8J73/rr9fU+tWrV8RkZGTK+byNLT02nevHnYZcRNlM8vyucGOr9k55zLU7DHPBXjve/qva/kva8GXA18dqBQFxGR+NI6dhGRiCkW5MG89+lAepDHFBGRQ6MRu4hIxCjYRUQiRsEuIhIxCnYRkYhRsIuIRIyCXUQkGRzCzaQKdhGRRPfDD6w/o22en65gFxFJVNnZbP3nC2yrXodSX6fl+dsU7CIiiWjmTH6udTqlut1PWtY59LljXp6/VcEuIpJItm4ls1MXsk9rxI7//UCXqsMoN20cXV6tmudDBNpSQERE8s9/PJ7MDrdzxPolDCxyKxu7PcM/epSlePFDO46CXUQkbGvWsOX2Byk9agirqMnLddPo9G5zTjopf4dTsIuIhCUnh5zX+rHzoa4U+/VXni7enaOeeZSX7i1BkRgmyhXsIiJhmDmTX2+8nZKzvuZLWjKk2at0f6cW1arFfmhdPBURKUiZmWTf/yA5DRuxedZi/q/U2ywf+CmvTwwm1EEjdhGRguE9jBrFjtvv5rC1K+hHR6a0e4qn+pXl+OODfSmN2EVE4m3ZMrLbXgKXXcaCtWW4uMwkUkf2441RwYc6aMQuIhI/O3fC88+T3fMxtm+HHvRi4w338vbzxSlTJn4vq2AXEYmHSZPI7ng7RefNZSzteKbiy/x9YBXOPz/+L62pGBGRIP30E9x2G5x1FmsyNnIp7/PZPaMYv6BgQh00YhcRCUZODrz9NjmdH8L/+BO96cyw6j15aWAKzZoVbCkasYuIxGr2bDjnHLjxRqZvrEFjN52N3XoxaVbBhzpoxC4ikn8bN0LPnvg+fdhUtAz3M4DZ9W5k4MAi1K8fXlkasYuIHCrvYfBg/Ekn4V96iTeK3sZJZFDr6ZuZMjXcUAeN2EVEDs28eXDXXZCezsIjG3GdH03Jpo1J/zfUqhV2cUYjdhGRvMjMhIcfxp96KtumzuKew/rSJHsKN73SmPT0xAl10IhdROTAvId334X774eVKxl73C3cvO4pGrVJZU4/qFIl7AL/SCN2EZH9yciA1q3hqqtYk5PKOcW+4sas1+k9KJUPPkjMUIcAgt05V8I5N9U5N8s5951z7rEgChMRCc3WrdCtG5x8MlmTp/JUhZepuHoa5S8/g3nzoEMHcC7sIvcviKmY7UBL732mc6448KVz7kPv/ZQAji0iUnC8h/fegwcfhGXLmFG3AxfP7wVHlGPkf+HSS8MuMG9iHrF7k7nrw+K7Hj7W44qIFKg5c+C886B9ezKLHsm1lT6n4XeDuPCmcsyblzyhDuC8jz2DnXNFgelADeAV732XfTynI9ARIDU1teHw4cNjft1ElZmZSUpKSthlxE2Uzy/K5wY6v30ptnEjJ7zxBhXGjGFn6RQGVHuEu+c8xHHld/Lggxk0bPhLfIrNhxYtWkz33jc66BO994E9gKOBNKDegZ5Xs2ZNH2VpaWlhlxBXUT6/KJ+b9zq/vezc6X2fPt6XLet90aJ+SdtOvl6FH71z3t93n/eZmXErM9+AaT4PWRzoqhjv/S9AOtAmyOOKiAQqLQ1OOw06dWJH3fp0vWAmJ4x5mZyjy/LVV/D881C6dNhF5l8Qq2JSnXNH73q/JNAKWBDrcUVEArd0KbRvDy1b4jdvZtIDI6k0/1P+9VE9uneHGTPg9NPDLjJ2QayKKQ+8tWuevQgw3Hs/NoDjiogEY8sWeOYZ6NULihRh00NPcOv8BxnRuySNGsGnA+CUU8IuMjgxB7v3fjbQIIBaRESC5T385z/w0EOwYgX+mmsYWv9Z7vxnJbZvt5y/7z4oFrF78HXnqYhE07ffwrnnwjXXQGoqK/8zkVZrh/DXLpWoX99WN3buHL1QB/WKEZGoWbkSuneHN9+EY44hp29/Xsq8mUdvLEqxYtC3r+1cVyTCw1oFu4hEQ2Ym1d54wxp2ZWXBgw8y77Ju3HT/0UydChddZKFeqVLYhcZfhH9niUihkJUF//431KhBtUGDoF07dsxewGMpvajf/GgWL4YhQ2DMmMIR6qARu4gkK+/ho4/swuh330GzZkzv0YPsRndyS3uYO9em1198EVJTwy62YGnELiLJZ9YsOP98uPBC2L4dRo5k68cTefbztpxxBvz8M4webSP1whbqoGAXkWSyciXcfDM0aGB3E734Inz3HWllLufkUxzDh1fm1lttAN+2bdjFhkfBLiKJb/Nm6NEDTjwRBg+2trqLFrHxhnv4v7sPo2VL64/eu/dM+vWDo44Ku+BwKdhFJHFlZUH//hboTzwBl1wCCxZAr16M+bIMderA66/bevTZs6FBg1/CrjghKNhFJPHkbnhx6qnwf/8HNWrAlCkwdCjrU07g2muhXTsoWxYmT7Y7SEuVCrvoxKFgF5HE4T18/DE0aQJXXAHZ2bYufeJEfJOmDBkCtWvbpx57DKZPt6fK3hTsIpIYvvwSmjeHNm1g/Xp44w1bs3jFFaxY6WjXDv76V6he3a6b9ugBhx0WdtGJScEuIuGaMcOWLZ59NixcCH36QEYG3HgjOUWK0a8f1KkDEyZA797w1VdQr17YRSc23aAkIuGYP9+G3e++a5PlzzwDnTr9Nln+/ffW0+Xzz6FlS7uGWr16yDUnCY3YRaRgLVkCN95ow+6PPrJwX7wYHn4YSpUiKwv+9S/rj/7tt9Yt4NNPFeqHQiN2ESkYq1bBk09aUhctCg88AF26wLHH/vaU2bPhlltg2jRb9fLqq1CxYog1JykFu4jE148/2jTLyy/buvTbboNu3fZK7O3b4Z//tEeZMjB0KPzlL3bTkRw6BbuIxMe6dfDcczbs3rIFOnSAnj3hT3/a62lTptgofd48uO4620h6j0G85IPm2EUkWCtX2n5z1arZZHnbtrZs8a239gr1LVtsNubMM2HTJhg7Ft5+W6EeBI3YRSQYy5bB00/DwIF2Y1GHDtC1K9Ss+YenTphgMzJLlsAdd9i3HXlkCDVHlIJdRGKzaBE89RQMGmST4jffbBdFTzjhD0/95Rfr6zJggLV/SU+3bUklWAp2Ecmf+fNtlcvQoXYL6J132qYX+9mmaNQoG52vXWsrG//+dyhZsmBLLiwU7CJyaGbNgn/8A0aOtJuJHnjA2ugef/w+n752LdxzDwwfbmvTR4+GRo0KuOZCRsEuInnzzTcW6KNH24T4o4/aRdL9XO303lqn33svZGZa190uXaB48YItuzBSsIvI/nkP48fbssVPPrFF5o8/DnffDUcfvd9vW77cuu1++CGcfrrNqdepU3BlF3YKdhH5o+3bbbjdu7ftM1ehgi1dufNOOOKI/X5bTg707Wsj85wceOEFa/9StGjBlS4KdhHZ04YNlsx9+tjk+Cmn2Przq68+aI/chQvh1lth4kRo1cqadu1jYYwUgJhvUHLOVXbOpTnn5jvnvnPO3RtEYSJSgBYutCUrVapA9+5w2mnWeWvmTLj++gOGelaWdQw45RSYM8eWsY8fr1APUxAj9izgQe/9DOfcEcB059wn3vt5ARxbROLFe/jiC5s/HzvWwrtDB7j//jxPiM+aZcvWZ8yAyy6DV16B8uXjXLccVMwjdu/9au/9jF3vbwbmA+rHJpKodu60teeNG9uORZMn2yh92TLrvJiHUN+2Df72N1u2uGIFjBhhqx8V6onBee+DO5hz1YAvgHre+02/+1pHoCNAampqw+HDhwf2uokmMzOTlJSUsMuImyifX5TPrVhmJmXfe48/jRtHiXXr2Fq5Msvbt2dt69bkHH54no8zd+6R9OpVix9+KE3r1mu4445FHHVUVhwrz7so//wAWrRoMd17f/C7ALz3gTyAFGA6cPnBnluzZk0fZWlpaWGXEFdRPr9IntvUqd7ffLP3JUt6D963aOH9mDHeZ2cf0mE2b/b+7ru9d877KlW8/+ijONUbg0j+/PYATPN5yONAVsU454oDI4HB3vv3gjimiMRgyxabbunbF6ZPtztEr7uOaY0b0+i22w75cOPHQ8eONlvTqZP1TT/AqkcJWRCrYhwwAJjvve8de0kikm9z51ryVqhg7RO3b7crmqtWQf/+ZJ544iEd7qef4KaboHVrKFHCljK+/LJCPdEFMWJvBnQA5jjnZu763KPe+w8COLaIHMy2bbYhdN++MGkSHH44XHUV3H47nHFGvrcheu89uOsuWL/euu/26GHhLokv5mD33n8JaAMrkYL2/fd2F9Abb9j2czVq2MYWN9wQ024Va9bYoH/kSKhfHz74ABo0CK5siT/deSqSTHbutCZcffvaDUTFisGll9rovEULKJL/2VXv7SbTBx6ArVttHr1zZzXtSkYKdpFE5z18+62l7pAhdtt/lSrWafHmmwNZPL50qTXtGj8emjWD11+Hk06KvXQJh4JdJFGtXg3vvGM7E82da3eGXnKJTbW0aRNIZ62cHLu22rWrTcX36WOdBWIY+EsCULCLJJJff7Wtht56y4bPOTnW9/a11+Avf7G2uQFZsMCadk2aZKte+vWDqlUDO7yESMEuEjbvLV3fesu2Gdq0CSpXtmH09dfvczPoWOzcCb16wWOPQenS9rIdOuR78YwkIAW7SFiWLLFplkGDYPFiS9n27S3MmzePy3zIwoUp3H+/NW1s396mXsqVC/xlJGQKdpGCtHatLRAfNsw6Kzpnq1l69oTLL4c49Tn59Vfb+OjZZxty3HFWwmWXxeWlJAEo2EXibc0aS9IRIyzMc3JsycmTT8J119kKlzj68ku45RZruX7BBWsYPLh8kFP1koAU7CLxsHq13eEzYoTdh+891K5tvW6vvBLq1o37pPbmzTZN/8orUK2abVlarFgGZcqot27UKdhFgrJq1e4w//JLC/O6dW2apX17e7+AfPSRrUtfvhzuvdeWvKekQHp6gZUgIVKwi8RixYrdYf7VVxbm9erB3/9uI/PatQu0nB9/tDtHBw2yl540ydrFSOGiYBc5FN7DvHm2ldzo0RbmYBt+Pv64jcxDuGXTe+sD1qmTdWTs3h26dbN+YFL4KNhFDmb7dpvDGDvWHkuX2ucbNLA5jiuvDHyt+aFYvRruvBPefx8aNrT7mk49NbRyJAEo2EX2Zc0aGDfOgvyTT2zjipIloVUruyJ50UVQMdytfb23xo4PPGC/e5591vahLqb/qws9/RMQgd8abVV96y14+GH45hv7fOXKdsPQxRfbevOSJcOtc5clS2xHo08/hXPOsT2oQ/yjQRKMgl0Kr40bbYpl3Dh7rFpFNeesN8uTT1qYn3xyQt1rn51td4s++qj1AHvtNQt4Ne2SPSnYpfDYvh0mT7Zh7oQJNirPzrZ93lq3hrZt+erII2l26aVhV7pP8+ZZ067Jk+HCC60le+XKYVcliUjBLtGVnW1NUXKD/Msv7d76okWhcWObKz/vPDjzTGuJC+xMwIXeO3fCM8/AE0/Y76B33oFrr02oPyQkwSjYJTq8t+3iJkywME9Lg59/tq/VrWubO7dqZZPSRx0Vbq15NG2atQOYPRuuvhpefBGOOy7sqiTRKdgleXkPy5bZLfuffWZhvmKFfa1KFdsy7rzzoGXLQHYZKki//mo3rD73HBx/vLVob9cu7KokWSjYJXnk5NhOQhMn2rTKxImwcqV9rWxZC/DzzrNRefXqSTtX8fnnNpe+aJH9kfHss3D00WFXJclEwS6Ja9s2m4vIDfJJk2wlC0CFCnD22fY46yxbvZLkS0M2bYIuXeyi6J/+ZDNKLVuGXZUkIwW7JI5ffrFb9HODfOpU2LHDvla7tm0Nd9ZZ9qhWLWlH5PvywQfWtGvVKrvh6PHHbd8NkfxQsEs4du6EOXNsyeHUqfb47jubNy9WzO6Nv+ceC/FmzeDYY8OuOC42bID77oPBg+367rvvQtOmYVclyU7BLvGXk2MTxnuG+MyZNtUCcMwx0KSJ9Vw5+2xLtlKlQi053ry37U3vvtv+UOnZ02462rXqUiQmCnYJ3qpVe4f4tGmWXmCB3bChda1q0sQeEZtWOZiVK+30R4+25fQDBtglApGgKNgl/3KXG86cufsxbdrulSpFi1o726uusgBv3Bjq1Cm0Xaq8h9dfh86dbSbqX/+yaZiiRcOuTKImkP/DnHMDgYuBdd77ekEcUxLMjh12T/vMmdQYO9Y2kpg5c/cqFeegVi27+adpUwvx+vUjP6WSV//7n/V0+ewzaN7cmnbVqBF2VRJVQQ2d3gT6AIMCOp6E6eefYdasvUfi8+bZMBMoX6KEhfY119jb+vVt1yAt4/iD7Gy7W/Rvf4PixaF/f1ujXohmniQEgQS79/4L51y1II4lBWjjRgvsefNsRUru+8uX737O8cdbcF9wwW8hPnHFCpqfd15YVSeNuXOtHcDUqdC2rXViDLmFuxQSznsfzIEs2MfubyrGOdcR6AiQmpracPjw4YG8biLKzMwkJSUl7DJ+U2zzZkotXUrpZcsotWwZpZcupfTSpRy+YcNvz8k+/HC2VqnClmrV2FqtGptr1GBLjRrsKFv2D8dLtPMLUhDntnOnY8iQKrzzTlVKl87innsW0aLFuoQYpUf5ZwfRP78WLVpM9943OtjzCizY91SrVi2fkZERyOsmovT0dJo3b16wL+q9XbRcuBAyMmDBgt0j8dWrdz+vVCm72aduXbuQmfu2atU8X8UL5fwKSKznNnWqjdLnzrUOjC++mFhL8KP8s4Pon59zLk/BXjiXJySzTZssvHMDPPfx/fe2fVuu0qUtsFu3tre5IV6lStLfep+Itm6FHj3g+eet39iYMbZPh0gYFOyJaNs22zB50aK9A3zhwr1H30WK2BrwWrXg3HPtbe6jQgVdoSsgaWl2QXTxYrj9duudfuSRYVclhVlQyx2HAs2BY51zK4Ce3vsBQRw7kryH9estCRYvtrVwe76fuw481zHHWFi3aWMbW+aGd/XqcPjh4ZyDsHGjbY/av78tXUxPt9+vImELalXMNUEcJ1K2bIEffrAbeH4f3osXQ2bm3s+vUMGCulUra+1Xvbq9rVnTgl0SypgxNjpfswYeesiW9WvJviQKTcXkR+6Ie9my3eG9x9tmixbZXPieSpTYHdgtWtj7uR9XqwYlS4ZyKnJo1q+33mTDhlkbgFGjoNFBL2WJFCwF++95b39jr1xpu/Hkvt0zwH/4YXcDq1wpKbaypEoV1lWqRMXTT//tY6pXt/XgmvNOWt7D0KEW6ps2WVvdLl3UtEsSU+EK9uxsWLdud2D/Prxz327d+sfvPf54C+pTT7U9yqpU2R3cVavaFje7gvv79HQqRnjJVWGzfDnccQeMGwenn279XurWDbsqkf2LRrBnZcHatbZi5GCP7Oy9v7dYMbsdsGJFu7PyoougUqXdn6tUydav6SJloZOTYxdGH37Y/tm88AJ06qSmXZL4EjfYvbdWr2vX7v1Ys+aPYb1+vT3/94491kK5fHlbx71nWOe+TU3Vum75g++/t/1GP//ctlHt398uiYgkg1CCvciOHdbm7vehvedj3brd26LtqWhRmxYpX96mQZo23R3e5cvv/lq5cpoAlUOWne3o1ctuNjr8cOuVftNNujwiySWUYC+9dKkNg36rohgcd5yFcbly1imwXDkL6dzP5T6OOUYjbImL2bPhrrsakJEBl1wCr75qq1BFkk0owb6tfHnb5DE3rMuUUVhLaLZvhyefhKeegiOOKMHw4dC+vUbpkrxCCfadRxxha7lFQjZlijXtmjcPOnSAK66YyiWXnBV2WSIx0TBZCqUtW+D+++HMM2HzZvjgAxg0CI46Kivs0kRilrirYkTi5NNPbZu6JUvgrrtyp2DCrkokOBqxS6Hxyy827fLnP9v1+i++gD59FOoSPQp2KRTef99uZXjrLXjkEdvS9eyzw65KJD40FSORtnYt3H03jBhh3SDGjIGGDcOuSiS+NGKXSPIe3n7bRumjRtlyxm++UahL4aARu0TODz9Yr/QPP4QzzrC7R2vXDrsqkYKjEbtERk6O3S1at65dGH3pJZg4UaEuhY9G7BIJGRnWtGviRFv10r+/7V8iUhhpxC5JLSsLnn7aLozOmQNvvAEff6xQl8JNI3ZJWjNn2rr0GTPg8svhlVesb5xIYacRuySdbdvgb3+Dxo1t06t334WRIxXqIrk0Ypek8tVXNkpfsABuuAF694ayZcOuSiSxaMQuSSEz0zaSPuss25L2o4/gzTcV6iL7omCXhDd+vO290qeP7Tk6dy60bh12VSKJS8EuCeunn2xbutatoUQJW8r40ktq2iVyMAp2SUgjR1o7gLffhkcftRUwzZqFXZVIctDFU0koa9ZY065334UGDWwuvX79sKsSSS6BjNidc22ccxnOuUXOuUeCOKYULt7bxdA6dawD41NPwddfK9RF8iPmYHfOFQVeAS4A6gDXOOfqxHpcKTyWLrV59Jtusj4vs2ZZz/TixcOuTCQ5BTFibwIs8t4v9t7vAIYBlwRwXIm4nBx4+WVb8TJ5sq16+fxzqFUr7MpEklsQc+wVgeV7fLwCaBrAcSXC5s+HW2+1G47atIG+faFq1bCrEomGIILd7eNz/g9Pcq4j0BEgNTWV9PT0AF46MWVmZur89iMryzFsWGUGDapGyZLZPPLIIs4/fy1Lltjm0mHTzy65Rf388sx7H9MDOAP4eI+PuwJdD/Q9NWvW9FGWlpYWdglxld/zmz7d+/r1vQfvr7zS+zVrgq0rCPrZJbeonx8wzechl4OYY/8GONE5d4Jz7jDgamB0AMeViPj1V+jaFZo0seWM770Hw4dDuXJhVyYSTTFPxXjvs5xznYCPgaLAQO/9dzFXJpEwcaLNpS9caM27evWCMmXCrkok2gK5Qcl7/wHwQRDHkmjYvNmWLL76qm168ckn0KpV2FWJFA5qKSCB+/BDW4/+2mtw333WtEuhLlJwFOwSmB9/hOuvhwsvhJQUmDQJnn8eSpcOuzKRwkXBLjHzHkaMsHYAQ4dC9+7w7bdwxhlhVyZSOKkJmMRk9Wq48054/31o2NDm0k85JeyqRAo3jdglX7yHgQOhdm3rwPjsszBlikJdJBFoxC6HbNWqEvz5zzBhApxzDrz+Opx4YthViUguBbvkWXa2Ne3q2rUxxYvbqpeOHaGI/u4TSSgKdsmTefPsBqMpU6Bp018YMeIYKlcOuyoR2RcFuxzQjh3wzDPwj3/YXqODB0P58nOoXLl52KWJyH7oj2jZr2nToHFj6NEDLr/cWu1eey24ffXzFJGEoWCXP9i6FR5+GJo2hQ0bYNQoW5+emhp2ZSKSF5qKkb18/rk17Vq0CG67zZp2HXVU2FWJyKHQiF0A2LQJ7rgDmje3LesmTID+/RXqIslIwS6MG2dNu/r3hwcegDlzoGXLsKsSkfxSsBdiGzbAddfBxRfbyPyrr+C556BUqbArE5FYKNgLIe9h2DBrBzB8OPTsCTNm2MVSEUl+unhayKxcaU27Ro+2pYwDBsDJJ4ddlYgESSP2QsJ7+Pe/rbXuJ5/YlMvkyQp1kSjSiL0Q+N//bOliWhq0aGEBX7162FWJSLxoxB5h2dnQu7eNyqdPt1UvEyYo1EWiTiP2iJo715p2TZ0KbdtaJ8aKFcOuSkQKgkbsEbNjBzz2GJx2GixebK0ARo1SqIsUJhqxR8jUqTZKnzvXmnW9+CIce2zYVYlIQdOIPQK2boXOnW3z6J9/hjFjrL2uQl2kcNKIPcmlpVnTrsWL4fbbrXf6kUeGXZWIhEkj9iS1caNtS9eypW1Nl55uF0gV6iKiYE9CY8bYjUYDBsBDD8GsWXDuuWFXJSKJQsGeRNavh2uugXbt4Jhj4Ouv4dln1bRLRPYWU7A75650zn3nnMtxzjUKqijZm/cwZIg17Ro5Eh5/3Lata6T/4iKyD7GO2OcClwNfBFCL7MPy5XaD0V//CjVqwLffQvfucNhhYVcmIokqpmD33s/33mcEVYzslpMD/frZBhhpafD88zBpkn0sInIgznsf+0GcSwc6e++nHeA5HYGOuz6sh432o+pYYEPYRcRRlM8vyucGOr9kV8t7f8TBnnTQdezOuU+B4/fxpW7e+1F5rcZ73x/ov+uY07z3kZ0h1vklryifG+j8kp1zbr+D5z0dNNi9961iL0dERAqKljuKiERMrMsdL3POrQDOAMY55z7O47f2j+V1k4DOL3lF+dxA55fs8nR+gVw8FRGRxKGpGBGRiFGwi4hETOjB7pzr7JzzzrnIdA93zj3hnJvtnJvpnBvvnKsQdk1Bcs71cs4t2HWO/3XOHR12TUGKaqsM51wb51yGc26Rc+6RsOsJknNuoHNunXMucvfHOOcqO+fSnHPzd/27vPdg3xNqsDvnKgN/Bn4Is4446OW9P8V7Xx8YC/QIuZ6gfQLU896fAiwEuoZcT9Ai1yrDOVcUeAW4AKgDXOOcqxNuVYF6E2gTdhFxkgU86L2vDZwO3HWwn13YI/bngYeBSF3B9d5v2uPD0kTv/MZ777N2fTgFqBRmPUGLaKuMJsAi7/1i7/0OYBhwScg1BcZ7/wXwU9h1xIP3frX3fsau9zcD84ED7mIc2g5Kzrl2wErv/SznXFhlxI1z7kngemAj0CLkcuLpZuA/YRchB1URWL7HxyuApiHVIvnknKsGNAC+PtDz4hrsB2pHADwKnB/P14+ng7Va8N53A7o557oCnYCeBVpgjPLSSsI51w37M3FwQdYWhKBaZSSRfY2eIvWXZNQ551KAkcB9v5sV+IO4Bvv+2hE4504GTgByR+uVgBnOuSbe+zXxrCkoh9BqYQgwjiQL9oOdn3PuBuBi4DyfhDdDFMJWGSuAynt8XAlYFVItcoicc8WxUB/svX/vYM8PZSrGez8HOC73Y+fcUqCR9z4SXdmccyd677/f9WE7YEGY9QTNOdcG6AKc673fGnY9kiffACc6504AVgJXA9eGW5LkhbPR7wBgvve+d16+J+yLp1H1tHNurnNuNjbddNDlSUmmD3AE8MmuJZ19wy4oSDG0ykhYuy52dwI+xi6+DffefxduVcFxzg0FJgO1nHMrnHO3hF1TgJoBHYCWu/5/m+mcu/BA36CWAiIiEaMRu4hIxCjYRUQiRsEuIhIxCnYRkYhRsIuIRIyCXUQkYhTsIiIR8/9gXP3YUuGv/QAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import numpy as np\n", + "from scipy.stats import norm\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "\n", + "# Resolution\n", + "sigma_thetahat = 1.0\n", + "\n", + "confidenceinterval = 0.95\n", + "\n", + "def denominator(mu): \n", + " return 1-norm.cdf(0, mu, sigma_thetahat)\n", + "\n", + "def numerator(mu, up):\n", + " return norm.cdf(up, mu, sigma_thetahat) - norm.cdf(0, mu, sigma_thetahat) \n", + "\n", + "prob_left = 1 - confidenceinterval\n", + "\n", + "theta_min = -4.0\n", + "theta_max = 4.0\n", + "\n", + "theta_obsmin = -4.0\n", + "theta_obsmax = 4.0\n", + "\n", + "thetas = [] # array to collect the scanned theta\n", + "lbounds = [] # array to collect the left bounds\n", + "bbounds = [] # bayesian bounds\n", + "\n", + "nsteps = 30\n", + "step = 0.2\n", + "\n", + "print (\"theta_obs\", \"numerator\", \"denominator\", \"Ratio\", \"limit\")\n", + "for i in range(nsteps+1): \n", + " \n", + " theta = theta_min + i/nsteps*(theta_max-theta_min)\n", + " thetas.append(theta)\n", + " \n", + " # upper limit\n", + " left_bound = norm.ppf(prob_left, loc=theta, scale=sigma_thetahat)\n", + " # print (theta, left_bound, right_bound)\n", + " lbounds.append(left_bound)\n", + "\n", + " # bayesian upper limit: solve numerically to find the upper limit\n", + " up = 0.001 \n", + " scanstep = 0.001\n", + " while (numerator(theta,up)/denominator(theta) < confidenceinterval) and (up < theta+5) :\n", + " # print (numerator(theta,up), denominator(theta), numerator(theta,up)/denominator(theta), up)\n", + " up+=scanstep \n", + " stringa = \"{:.6f}\".format(theta) + \" {:.6f}\".format(numerator(theta,up)) + \" {:.6f}\".format(denominator(theta)) + \" {:.6f}\".format(numerator(theta,up)/denominator(theta)) + \" {:.6f}\".format(up) \n", + " print (stringa)\n", + " bbounds.append(up)\n", + " \n", + " \n", + "plt.plot(lbounds,thetas, 'b-')\n", + "plt.plot(thetas,bbounds, 'r-')\n", + "plt.axis([-4,2,-1,5])\n", + "plt.xticks(np.arange(-4,3, 1.0))\n", + "plt.yticks(np.arange(-1,6, 1.0))\n", + "\n", + "plt.grid()\n", + "plt.show()\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.9772498680518208\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtX0lEQVR4nO3deXiU1dnH8e+djSWEECHIEjYRBbQgISIiUGUTlK0VNYBiWS6gAlXrjhWx2lYrWq0iElARAQOIuKLQvgpSlCUsIghCZBcsgbCHkITc7x8z2AiBTCDJmeX+XFcu8szznJnfALlz5syZc0RVMcYYE7zCXAcwxhhTuqzQG2NMkLNCb4wxQc4KvTHGBDkr9MYYE+QiXAcoTLVq1bR+/fquYxhjTMBYuXLlPlWNL+ycXxb6+vXrk5aW5jqGMcYEDBHZfrZzNnRjjDFBzgq9McYEOSv0xhgT5KzQG2NMkPOp0ItIVxH5XkTSReSRc1x3tYicFJE+xW1rjDGmdBRZ6EUkHBgPdAOaAn1FpOlZrnsWmF/ctsYYY0qPLz36VkC6qm5R1RwgFehVyHWjgDnA3vNoa4wxppT4Mo++NrCzwPEu4JqCF4hIbeA3QAfg6uK0NSYQHD58mHHjxgEQFRVFt27dSExMREQcJzOmaL4U+sL+J5++iP2LwMOqevK0//i+tPVcKDIUGApQt25dH2IZU3oOHDjAyy+/TEJCAoMGDeLo0aM8/fTTAKgqjz/+OE2aNOGVV16hQ4cOjtMac26+DN3sAuoUOE4Adp92TRKQKiLbgD7AqyLS28e2AKhqiqomqWpSfHyhn+I1pkwsWrSIZs2aMXbsWJYuXQpArVq1yM/PJz8/n8zMTCZOnEjVqlWpWrUq4Cn+xvgrXwr9CqCRiDQQkSggGfiw4AWq2kBV66tqfeBd4G5Vfd+Xtsb4i9zcXP70pz9xww03UKFCBZYvX05KSsoZ18XFxTF06FAWL15M8+bNUVVuvfVW/vrXv3Ly5EkHyY05tyILvarmASPxzKbZAMxS1fUiMlxEhp9P2wuPbUzJ+/e//81f/vIXBg4cyKpVq0hKSvKp3YkTJ4iIiOCxxx5j0KBB5Ofnl3JSY4pH/PElZ1JSktqiZsaFDRs20KRJk2K3U1WeeuopnnjiCUaOHMk///lPe6PWlCkRWamqhfZO/HL1SmPKiqry6KOP0rlzZzp27HheRR5ARHj88cc5dOgQL7zwAvHx8YwZM6aE0xpzfqzQm5D297//nWeffRaAjh07XtB9iQjjxo1DROjWrVtJxDOmRNjQjQlZS5cupW3bttxyyy2kpqaWylBLVlYWFStWLPH7NeZ05xq6sUXNTEg6fPgw/fr1o06dOqSkpJRKkX/yySdp1aoVx48fL/H7NqY4rNCbkDRlyhS2b9/O9OnTiY2NLZXHuPbaa1m/fj33339/qdy/Mb6yQm9C0qhRo1ixYgVt2rQptcfo0qUL999/PxMmTODDD+3jI8YdG6M3IeXgwYMcPny4zJbZyMnJ4eqrr+bAgQNs2LCB6OjoMnlcE3psjN4YryeeeIIrrriC/fv3l8njRUVFMX78eA4cOMDKlSvL5DGNOZ1NrzQh49tvv2X8+PEMHTr05zVqykLbtm3ZuXMnVapUKbPHNKYg69GbkKCqjBo1itjYWJ566qkyf/wqVaqgqixcuLDMH9sYK/QmJMycOZNFixbx17/+tUx78wXNmDGDG264gXnz5jl5fBO6rNCbkLBu3TqSkpIYMmSIswy33norl112GQ899JAtfGbKlBV6ExKefvppvvrqK8LDw51liIqK4sknn2T9+vXMnj3bWQ4TeqzQm6B24sQJVq1aBUBkZKTjNJ5e/RVXXMHYsWNt7XpTZqzQm6D2xhtv0LJly5+LvWvh4eGMHTuWrKwstm3b5jqOCRH2gSkTtLKzs7n00kupX78+ixcv9pv14fPz88nLyyMqKsp1FBNELvgDUyLSVUS+F5F0EXmkkPO9RGStiKwRkTQRaVvg3DYR+fbUufN/GsYUz6RJk/jxxx/585//7DdFHiAsLIyoqCiys7PZsGGD6zgmBBTZoxeRcGAT0BnPZt8rgL6q+l2BayoBx1RVRaQZni0DG3vPbQOSVHWfr6GsR28uVG5uLpdccgmXXHIJCxcu9KtCf0r37t3ZvHkzGzZsICzMRlHNhbnQHn0rIF1Vt6hqDpAK9Cp4gaoe1f/9xogG/G88yISUb7/9lkOHDvHwww/7ZZEHuOOOO9i0aRMff/yx6ygmyPlS6GsDOwsc7/Le9gsi8hsR2Qh8AgwqcEqBBSKyUkSGnu1BRGSod9gnLSMjw7f0xpxFYmIiO3fupGvXrq6jnFWfPn2oW7cu48aNcx3FBDlfCn1h3aEzeuyqOtc7XNMbKPgZ8+tUNRHoBowQkfaFPYiqpqhqkqomxcfH+xDLmMIdPHgQVSU2Ntavh0QiIiK47777WLx4McuXL3cdxwQxX34KdgF1ChwnALvPdrGqfgk0FJFq3uPd3j/3AnPxDAUZU2r69evn1z35ggYPHkxsbCzvvvuu6ygmiPlS6FcAjUSkgYhEAcnAL3ZREJFLxTsQKiKJQBSwX0SiRSTGe3s00AVYV5JPwJiC1q1bx6effkr79oW+cPQ7MTExrF69+ucNyo0pDUUuU6yqeSIyEpgPhANvqOp6ERnuPf8acAswQERygePA7d4ZOBcDc72/AyKAGar6WSk9F2N44YUXqFChAsOHD3cdxWcNGjQAIC8vj4gIWznclDz7wJQJGvv27SMhIYFBgwbx6quvuo5TLFOnTuXRRx/l+++/p1KlSq7jmABkO0yZkPD2229z4sQJRowY4TpKsTVq1Ijdu3czffp011FMELIevQkaJ06c4IsvvgiYN2ILUlUSExPJz89nzZo1fjv33/gv69GbkFCuXLmALPIAIsLdd9/N2rVr+eqrr1zHMUHGCr0JCiNHjmTixImuY1yQfv36UblyZSZMmOA6igkyVuhNwNu+fTsTJkxg586dRV/sx6Kjo5k0aRIPPvig6ygmyNhcLhPwUlJSABg69KwrbASM2267zXUEE4SsR28CWm5uLpMnT+bmm2+mbt26ruOUiLVr13LvvffavrKmxFihNwHtk08+Ye/evUHRmz9l3bp1vPTSSyxcuNB1FBMkrNCbgFa9enX69+8fsLNtCvOb3/yGKlWq8Prrr7uOYoKEFXoT0Nq0acO0adOCaumAChUq0K9fP+bMmcOBAwdcxzFBwAq9CVhLlixh+/btrmOUisGDB3PixAlSU1NdRzFBIHi6QSakqCqDBg2iRo0aLFq0yHWcEteiRQs6depEXl6e6ygmCFihNwFpyZIlbNq0idGjR7uOUipEhH/961+uY5ggYUM3JiC9/vrrxMTE0KdPH9dRSlV+fj7p6emuY5gAZ4XeBJyjR48ya9Ysbr/9dqKjo13HKVV/+MMfuOaaa8jJyXEdxQQwK/Qm4CxbtowTJ07wu9/9znWUUnfzzTeTmZnJvHnzXEcxAcynQi8iXUXkexFJF5FHCjnfS0TWisgaEUkTkba+tjWmuDp27MiePXto06aN6yilrnPnzlx88cW8/fbbrqOYAFZkoReRcGA80A1oCvQVkaanXfZ/QHNVvQoYBEwuRltjfHZq/4T4+PiQWLM9IiKCfv368dFHH5GZmek6jglQvvToWwHpqrpFVXOAVKBXwQtU9aj+bweTaEB9bWtMcYwbN47rr7+e7Oxs11HKzIABA8jNzeXdd991HcUEKF8KfW2g4Pqvu7y3/YKI/EZENgKf4OnV+9zW236od9gnLSMjw5fsJsSoKm+99Ra5ubmUL1/edZwy07x5c/79738zcOBA11FMgPKl0Bf2+viM/QdVda6qNgZ6A08Vp623fYqqJqlqUnx8vA+xTKhZs2YN69evZ8CAAa6jlCkRoWPHjkRGRrqOYgKUL4V+F1CnwHECsPtsF6vql0BDEalW3LbGnMvUqVOJiooKyTXbVZUxY8bwyiuvuI5iApAvhX4F0EhEGohIFJAMfFjwAhG5VLzvjIlIIhAF7PelrTG+yMvLY8aMGfTo0YO4uDjXccqciLB06VJeeOEF/vd2mDG+KbLQq2oeMBKYD2wAZqnqehEZLiLDvZfdAqwTkTV4Ztncrh6Fti2F52GCXF5eHo899hijRo1yHcWZO+64g61bt/L111+7jmICjPhj7yApKUnT0tJcxzDGrxw5coTq1aszePBgG8IxZxCRlaqaVNg5+2Ss8XtZWVm8+eabHDp0yHUUp2JiYujZsyezZs0iNzfXdRwTQKzQG7/30UcfMWjQIFauXOk6inN33XUX119/vW1IYorFhm6M3+vVqxdpaWns2LGD8PBw13GM8Us2dGMCVmZmJp9++inJyclW5AvYtGkTx48fdx3DBAgr9MavzZkzh9zcXPr16+c6it9Yvnw5l19+OR9+aDOVjW+s0Bu/lpaWxuWXX05iYqLrKH6jZcuW1K5dmxkzZriOYgKEFXrj1yZOnMiyZctCYqVKX4WHh5OcnMynn35qK1oan1ihN37r1ESB2NhYx0n8T9++fcnNzWXu3Lmuo5gAYIXe+K0OHTrw+OOPu47hlxITE7n00kuZNWuW6ygmAFihN37phx9+YOHChdabPwsRYebMmaSmprqOYgJAhOsAxhRm5syZACG5UqWv7A1q4yvr0Ru/9M4773DddddRt25d11H82uzZs7n77rtdxzB+zgq98Tvr1q1j3bp19O3b13UUv5eens6ECRPYvn276yjGj1mhN36nSpUqPProo/Tp08d1FL93++23A9ibsuacbK0bYwLcNddcQ15eni36FuJsrRsTMNLT0/nss89sGd5iSE5OZtWqVWzatMl1FOOnfCr0ItJVRL4XkXQReaSQ8/1FZK336ysRaV7g3DYR+VZE1oiIddPNOU2cOJGePXty5MgR11ECxm233Ua7du1Cfr1+c3ZFTq8UkXA82wN2xrPZ9woR+VBVvytw2Vbg16p6QES6ASnANQXO36Cq+0owtwlCqsqsWbPo0qULF110kes4AaN27dp8+eWXrmMYP+ZLj74VkK6qW1Q1B0gFehW8QFW/UtVTOyEsBRJKNqYJBUuXLmXHjh0kJye7jhKQ9u/fz7591p8yZ/Kl0NcGdhY43uW97WwGA58WOFZggYisFJGhZ2skIkNFJE1E0jIyMnyIZYJNamoq5cqVo2fPnq6jBJzDhw+TkJDASy+95DqK8UO+FPrClg0sdKqOiNyAp9A/XODm61Q1EegGjBCR9oW1VdUUVU1S1aT4+HgfYplg89VXX3HTTTdRuXJl11ECTuXKlbnuuuuYOXMm/jiTzrjlS6HfBdQpcJwA7D79IhFpBkwGeqnq/lO3q+pu7597gbl4hoKMOcOyZcuYNGmS6xgBKzk5mc2bN7N69WrXUYyf8aXQrwAaiUgDEYkCkoFfbG0jInWB94A7VXVTgdujRSTm1PdAF2BdSYU3wSUsLIyqVau6jhGwfvvb3xIREfHzOkHGnFJkoVfVPGAkMB/YAMxS1fUiMlxEhnsvGwNUBV49bRrlxcB/ROQbYDnwiap+VuLPwgS0vLw8rrrqKqZOneo6SkC76KKL6Ny5sw3fmDP4tHqlqs4D5p1222sFvh8CDCmk3Rag+em3G1PQ559/zjfffENMTIzrKAHvmWeeITo62nbkMr9gyxQb52bOnEnlypXp1q2b6ygBr1mzZq4jGD9kSyAYp3Jycnjvvffo3bs35cuXdx0nKCxZsoRBgwZx8uRJ11GMn7BCb5xasGABBw8etA9JlaBdu3bx5ptv8p///Md1FOMnrNAbpxISEhg+fDidOnVyHSVodO/enYoVK9rsG/MzW6bYmCCUnJzM559/zu7du4mIsLfiQoEtU2z80tq1a1mzZo1NBSwFt99+OxkZGXzxxReuoxg/YIXeOPPUU09x4403kp+f7zpK0OnWrRstW7bk2LFjrqMYP2Cv6YwTR44c4eOPP2bw4MGEh4e7jhN0ypcvjw1/mlOsR2+c+Oijj8jOzv55z1NTOk6cOMHevXtdxzCOWaE3TrzzzjskJCRw3XXXuY4StFSVK6+8kvvuu891FOOYFXpT5rKzs1myZAnJycmEhdl/wdIiInTo0IEPPviArKws13GMQ/ZTZspc+fLl2blzJ488csb2w6aE9e3bl2PHjvHxxx+7jmIcskJvnIiOjrYlictAu3btqFWrFu+8847rKMYhK/SmTO3Zs4fExEQWL17sOkpICA8P57bbbmPevHkcPHjQdRzjiBV6U6Zmz57N6tWrqVatmusoIWPkyJF88cUXxMbGuo5iHLF59KZMvfPOOzRv3pwmTZq4jhIyGjZsSMOGDV3HMA751KMXka4i8r2IpIvIGe+giUh/EVnr/fpKRJr72taEjq1bt7J06VL69u3rOkrI2bJlCyNHjuS///2v6yjGgSILvYiEA+OBbkBToK+IND3tsq3Ar1W1GfAUkFKMtiZEpKamAtiHpBzIyspi/PjxzJ4923UU44AvPfpWQLqqblHVHCAV6FXwAlX9SlUPeA+XAgm+tjWho1mzZtx///3Ur1/fdZSQc+WVV/KrX/2KGTNmuI5iHPCl0NcGdhY43uW97WwGA58Wt62IDBWRNBFJy8jI8CGWCTQ333wz48aNcx0jZPXv35+vv/6aLVu2uI5iypgvhb6wXYYLXVdWRG7AU+gfLm5bVU1R1SRVTYqPj/chlgkkixcvZvfu3a5jhLRT741Yrz70+FLodwF1ChwnAGf8xIpIM2Ay0EtV9xenrQlu+fn59O/fn6FDh7qOEtLq1q1Lr169ECms/2WCmS/TK1cAjUSkAfAjkAz0K3iBiNQF3gPuVNVNxWlrgt/ixYvZuXMnzzzzjOsoIe/99993HcE4UGShV9U8ERkJzAfCgTdUdb2IDPeefw0YA1QFXvX2FvK8wzCFti2l52L81PTp04mOjqZXL3sf3h+oKrt376Z27XO91WaCie0Za0rViRMnqFGjBjfffDPTpk1zHccAw4YN4+OPP2bHjh226UsQsT1jjTMrVqzg0KFD9O/f33UU49WhQwd2797NokWLXEcxZcQKvSlVbdu2ZdeuXXTq1Ml1FOPVo0cPYmJimD59uusopoxYoTelrlatWkRGRrqOYbwqVqxInz59mD17tm1IEiKs0JtSM3nyZLp27crhw4ddRzGnGTBgAEeOHLFZOCHCVq80pebNN9/k0KFDxMTEuI5iTtO+fXvmzp1L165dXUcxZcB69KZU/PDDD3z11Vfceeed9gEdPxQWFkbv3r0pX7686yimDFihN6Xi7bffRkRsto0fy8/P5+mnn2bq1Kmuo5hSZkM3psSpKlOnTqVjx44kJCQU3cA4ERYWxkcffUR2djYDBgxwHceUIuvRmxKXk5PD8OHDueeee1xHMUUYMGAAa9eu5ZtvvnEdxZQiK/SmxJUrV46HHnqI7t27u45iipCcnExkZCRvvfWW6yimFFmhNyXq2LFjzJgxg+PHj7uOYnxQtWpVevTowbRp08jJyXEdx5QSK/SmRM2ZM4f+/fuzYsUK11GMj4YMGULr1q3Zv39/0RebgGSLmpkSdf311/Pjjz+yadMmm1ZpTBmyRc1Mmfjhhx9YtGgRgwYNsiIfgLZu3UpmZqbrGKYUWKE3JWbKlCmEhYXZVL0AtHPnTho2bMjkyZNdRzGlwKdCLyJdReR7EUkXkUcKOd9YRL4WkRMi8sBp57aJyLciskZEbDwmiKWlpXHjjTfahhYBqE6dOrRp04Y333wTfxzONRemyEIvIuHAeKAb0BToKyJNT7ssE/gDMO4sd3ODql51tvEjExzmzZtHamqq6xjmPA0aNIiNGzeydOlS11FMCfOlR98KSFfVLaqaA6QCv9gTTlX3quoKILcUMpoAkJOTg4hQuXJl11HMebr11luJjo624Zsg5Euhrw3sLHC8y3ubrxRYICIrRWRoccKZwPDf//6XmjVrMnv2bNdRzAWIiYkhOTmZuXPncuLEidJ7oHvv9XyZMuNLoS9s+kRxBvGuU9VEPEM/I0SkfaEPIjJURNJEJC0jI6MYd29cmzJlCpmZmVx55ZWuo5gLNHbsWDZu3Ei5cuVK70HWrPF8mTLjS6HfBdQpcJwA7Pb1AVR1t/fPvcBcPENBhV2XoqpJqpoUHx/v690bx/Lz80lJSaF9+/Y0adLEdRxzgRISEqhevbrrGKaE+VLoVwCNRKSBiEQBycCHvty5iESLSMyp74EuwLrzDWv8z//93/+xZcsWhg0b5jqKKSGbN2+mXbt2LFu2zHUUU0KKXKZYVfNEZCQwHwgH3lDV9SIy3Hv+NRGpAaQBlYF8EbkXzwydasBc74dnIoAZqvpZqTwT40RKSgpVq1bllltucR3FlJAaNWqwevVqJk6cyDXXXOM6jikBtgSCuSDffPMNP/zwA7/97W9dRzElaOjQoUybNo09e/YQGxtbsnd+/fWePxcuLNn7DXG2BIIpNc2bN7ciH4SGDRvG8ePHmTZtmusopgRYoTfn5eTJk4wcOZK1a9e6jmJKQcuWLWnZsiUTJkywT8oGASv05rx89NFHjB8/ns2bN7uOYkrJmDFjGD58OCdPnnQdxVwg2zPWnJeXX36ZOnXq0KtXr6IvNgGpZ8+eriOYEmI9elNs69ev5/PPP+fuu+8mIsL6CsEsKyuLiRMnsn37dtdRzAWwQm+K7ZVXXqFcuXIMGTLEdRRTyvbv38+IESMYP3686yjmAlihN8VWs2ZNRowYQbVq1VxHMaWsTp069O7dm8mTJ5OVleU6jjlPVuhNsY0ZM4bnn3/edQxTRkaNGsWBAweYMWOG6yjmPFmhNz47efIkCxYsID8/33UUU4bat29Ps2bNeOmll2yqZYCyQm98NnfuXG688UY++8xWsQglIsJ9991HxYoV2bdvn+s45jzYEgjGJ6pK69atyczMZOPGjYSHh7uOZMpQfn4+IlIym77bEgil4lxLINjcOOOTxYsXs3z5ciZMmGBFPgSFhXle/O/bt4/s7GwSEhIcJzLFYUM3xifjxo2jWrVq3HXXXa6jGEdyc3O58sorefjhh11HMcVkhd4U6ejRo6xdu5aRI0dSoUIF13GMI5GRkfTv35+ZM2eyY8cO13FMMVihN0WqVKkSmzdv5sEHH3QdxTh2zz33APDiiy+6DWKKxQq9OaeDBw+SnZ1NZGQkFStWdB3HOFa3bl2Sk5NJSUmxGTgBxKdCLyJdReR7EUkXkUcKOd9YRL4WkRMi8kBx2hr/NmbMGBo3bkx2drbrKMZPjB49muPHj9s02wBS5KwbEQkHxgOd8WwUvkJEPlTV7wpclgn8Aeh9Hm2Nn9qzZw8pKSnccccdlC9f3nUc4yeaNm3Ktm3bqFOnjusoxke+9OhbAemqukVVc4BU4Bdr06rqXlVdAeQWt63xX8899xx5eXk8+uijrqMYP3OqyB86dMhxEuMLXwp9bWBngeNd3tt84XNbERkqImkikpaRkeHj3ZvSsnfvXl577TX69+9Pw4YNXccxfujFF1/kkksusWIfAHwp9IV9FM7Xj9P63FZVU1Q1SVWT4uPjfbx7U1pmzJhBdnY2o0ePdh3F+Kn27duTmZnJK6+84jqKKYIvhX4XUHAwLgHY7eP9X0hb49A999zDqlWruPzyy11HMX4qMTGR7t278/zzz3Pw4EHXccw5+FLoVwCNRKSBiEQBycCHPt7/hbQ1jmRlZSEiXHXVVa6jGD/35z//mQMHDjBu3DjXUcw5FFnoVTUPGAnMBzYAs1R1vYgMF5HhACJSQ0R2AX8E/iQiu0Sk8tnaltaTMRduy5Yt1K5dmw8++MB1FBMAWrRoQd++fZk4caJNwfVjPi1qpqrzgHmn3fZage9/wjMs41Nb47+eeOIJsrOzSUoqdBE8Y87w3HPPISI2BdeP2eqV5mdr165l+vTpPPTQQ9Su7evEKhPqTv1fUVWOHz9un6D2Q7YEgvnZ6NGjiY2NtdUJTbGpKjfddBODBw92HcUUwgq9AWDjxo3MmzePhx9+mLi4ONdxTIARERITE0lNTWXZsmWu45jTWKE3ADRu3JiVK1dy7733uo5iAtQjjzxCzZo1GTVqlO0r7Ges0BsyMzMBzwwKe0PNnK+YmBieffZZVqxYwdSpU13HMQVYoQ9xGRkZXHrppbz00kuuo5gg0L9/f1q3bs2rr76KP+5HHaps1k2I+9Of/sSRI0fo0qWL6ygmCISFhZGamkq1atVKZiNxUyKsRx/Cli9fzqRJkxg5ciRNmjRxHccEiXr16hEdHU12dja7du1yHcdgPfqQlZOTw+DBg6lVqxZPPvmk6zgmyKgqnTp1Ii8vjyVLlhAeHu46UkizHn2IWrNmDVu3bmXChAlUrlzZdRwTZESE3//+9yxbtozx48e7jhPyrNCHqFatWrF161Z69OjhOooJUv369aNbt26MHj2a7du3u44T0qzQh5j8/Hw++eQTVBVb99+UJhFhwoQJAAwbNsxm4ThkhT7EvPDCC3Tv3p358+e7jmJCQL169fjb3/5GZmam7UTlkPjjb9mkpCRNS0tzHSPorFq1itatW9O9e3fmzJlj099MmcjPz+fkyZNERkZ6brj+es+fCxe6ihSURGSlqha67Kz16ENEVlYW/fr1Iz4+nkmTJlmRN2UmLCyMyMhI9u/fz8MPP0y+H3Yug50V+hDxxz/+kU2bNjF16lSqVq3qOo4JQV9//TV///vf2bp1q+soIcenefQi0hV4CQgHJqvqM6edF+/5m4As4Hequsp7bhtwBDgJ5J3tpYUpXT169KBu3bp07NjRdRQTorp3786wgcMo92Yq/835kYtdBwohRY7Ri0g4sAnojGez7xVAX1X9rsA1NwGj8BT6a4CXVPUa77ltQJKq7vM1lI3Rl5zs7GxbqMz4jYzt+6hYvy5R5LL0jQW0G3iD60hB40LH6FsB6aq6RVVzgFSg12nX9AKmqsdSoIqI1Lyg1OaCZWZm0rx5cyZOnOg6ijEAxNerRlateghKpUEjOJaR5TpSSPCl0NcGdhY43uW9zddrFFggIitFZOjZHkREhopImoikZWRk+BDLnMvJkyfp168fW7dupVmzZq7jGPOz6CrRZFKR5mxkWb1b0Xx7c7a0+VLoC5uecfq/zLmuuU5VE4FuwAgRaV/Yg6hqiqomqWqSfZDnwqgq9957L/Pnz+fll1/m2muvdR3JmF+4iCy+pD0djs9jWs2BruMEPV8K/S6gToHjBGC3r9eo6qk/9wJz8QwFmVL07LPP8sorr/DHP/6RYcOGuY5jTKF+zSI+5jru3PsWLzUe6TpOUPOl0K8AGolIAxGJApKBD0+75kNggHi0Bg6p6h4RiRaRGAARiQa6AOtKML8pRLly5ejfvz/PPfec6yjGnJUAHVjGFyQy4vvXeK2jraJaWoqcXqmqeSIyEpiPZ3rlG6q6XkSGe8+/BszDM+MmHc/0ylOvxS4G5no/nBMBzFDVz0r8WRgADh06RGxsLPfddx+qah+KMn6vInlcwfd8w2X87vO/8dbtsdw1817XsYKOT/PoVXUenmJe8LbXCnyvwIhC2m0Bml9gRuOD9957jyFDhrBgwQKSkpKsyJuAUZ1jHOYnfiCBPrMeY1HDBvz6r6dP7DMXwj4ZGwRmzZrFbbfdRuPGjWnUqJHrOMYU26UcIJJD7KYGSX/rx5JHPnIdKahYoQ9wU6ZMoW/fvlx77bXMnz+f2NhY15GMOS+XsY8qHOJHapP0bB9SbhjrOlLQsEIfwD755BMGDhxIhw4d+PTTT4mJiXEdyZgLEs9+qrCXNVzGkIV/5pWGto59SbBCH8BuvPFG/vnPfzJv3jwqVarkOo4xJaI6h2hEOgu4mpFbUni7Ug+yDtknaC+EFfoAs2XLFnr06MFPP/1EREQEo0aN+t8638YEiYvIpiMrmE57BmR9wsq4dvz0zR7XsQKWFfoAMmfOHFq0aMF//vMfNm3a5DqOMaUqEqU/XzKFNrTUDXBVC5b/xXZGOx9W6APAoUOHGDZsGH369KFx48asXr2a9u0LXUnCmKDzO75iF7U5RjRJf+rGlLhb2LfT58VwDVboA8Jjjz3G5MmTefDBB1m8eDH169d3HcmYMnUZ6dRkN+/Slt8dfI/9dVuROtxWZfWV7Rnrp9LT08nJyaFp06ZkZGSwbds2rr76atexjLlgWVckEfXdGiI4eV7t59CCFuzjEnYyp1wnGrz7NxK7235GtmdsANm/fz/33nsvTZs25f777wcgPj7eirwxXrewmipk8Ba/pvuJL7m8x6/5V6tHObTjkOtofssKvZ/YvXs3DzzwAPXr1+fll19m4MCBvPnmm65jGeOXLiKbu1jEeqrzDc3pvOIZTtarz+SEAWxeZhMVTmeF3iFVJT8/H4Bp06bxj3/8g549e7J27VomTpxIjRo1HCc0xr8lsos2fM0GLieNyxjy49vUat2CnvVfZ86cDa7j+Q0r9A7s3buXV199lZYtW/L2228D8Pvf/57Nmzczffp0rrjiCscJjQksTfieLiznCy5lBu35ZHtf+vRpQo0a23n9dTh40HVCt6zQlxFVZdKkSXTq1ImaNWsyYsQITp48+fPaNDExMVxyySWOUxoT2G4gnb4sphfTgf+wd291hgyBatXyqVZtKQMHfsmaNaH3wSsr9KUkIyOD999/n5SUFABEhIkTJ7Jjxw5Gjx7N2rVrWbNmDb1793Yb1JggFEsE0BbVCgCo7mD//jpMmdKeFi1qUr78Bq6+eiGzZh0Jid6+T+vRm7MruMHH+++/z6xZs0hLS2Pz5s0AxMXFMWjQICIiIliwYAFxcXG2VrwxZSw/vz6ebaw3I/ITJ05UJi2tDbffHgXARRftp3Ll72jePI927WK58cbaXHFF9aD5WfWp0ItIV+AlPDtMTVbVZ047L97zN+HZYep3qrrKl7b+7tTnDESE9PR0Fi5cyNatW9m6dSs//PADGzduZNu2bcTFxbFmzRqWLFlCixYtGDJkCG3btqVly5ZERHj+mi+66CKXT8WYECdAI1RP7dmQC0BEBBw4kEFm5uVs21adDz6ABx6AsLADtGgRR9OmcOTIWuLiDnHZZRVo2jSW5s3jqVMnlrCwwPhFUGShF5FwYDzQGc8m4CtE5ENV/a7AZd2ARt6va4AJwDU+ti0xJ0+eJCsri9zcXHJzc8nJySEnJ4fq1asTExNDZmYmq1at4vjx4xw7doxjx45x9OhRevfuTb169Vi6dCnPP/88mZmZ7N+/n4yMDDIyMli6dCmJiYl8/vnnDBs2jPDwcOrVq0eDBg248847ycnJAeDxxx9n7NixpfHUjDElzrMYYF4eQGPvbfuA7YSFZSFShdWr41i5EuBKTh/pFsnh4oujqF4d9u5dR0TEYWJicomJyScmBho0qEabNr8iOho2blxFxYoQExNJTEwElSpFUqNGFWrVqkZUFBw7dpAKFSKoUqUi5cuX/Ii6Lz36VkC6d1tARCQV6AUULNa9gKneLQWXikgVEakJ1PehbYlJS0ujdevWZ9yemprK7bffzurVq+ncufMZ5xs2bEi9evU4evQo69evJy4ujjp16tCyZUuqV6/+c0/81ltvpUuXLiQkJPzcSy8oPDy85J+UMUFIKL1P5IeRfwGtqwHVyD/jLk4CPwIHgGOEheUA1di79wp++gk8AxY1gUpAZaAcAJMnn2qfWMTjVgEgJiaPgwfDCCvhWu9Loa8N7CxwvAtPr72oa2r72BYAERkKDAWoW7euD7HOVK9ePZ577jkiIyOJjIwkKiqKcuXKcc01nods2bIlX375JRUqVCA6OpqKFStSqVIlqlSpAkCnTp347ruz/w6Ki4sjLi7uvLIZYzzyr2yOfLeSo5TOHgqVOMqB8rWoVKLvQEYCdbxfZ8rPb0J2NgV+QeTgGRqK9h5v8N6W5/06CVwMNMTz3sGXgHLkyNXs2xdB9eolmd23Ql/YINTpv47Pdo0vbT03qqYAKeBZ68aHXGeoUaMGDzzwwFnPV6lShXbt2p3PXRtjSkilma/DzNdLqcx7vF+K9+2bKO/XKU3Oca0Avy7VNL4U+l388tdYArDbx2uifGhrjDGmFPkyErQCaCQiDUQkCkgGPjztmg+BAeLRGjikqnt8bGuMMaYUFdmjV9U8ERkJzMfzjsMbqrpeRIZ7z78GzMMztTIdz/TKgedqWyrPxBhjTKFsPXpjjAkCth69McaEMCv0xhgT5KzQG2NMkLNCb4wxQc4v34wVkQxgu/ewGp4FKPyZZSwZlrFkWMaSEWgZ66lqfGEX+WWhL0hE0s72TrK/sIwlwzKWDMtYMoIpow3dGGNMkLNCb4wxQS4QCn2K6wA+sIwlwzKWDMtYMoImo9+P0RtjjLkwgdCjN8YYcwGs0BtjTJALmEIvIqNE5HsRWS8if3ed53QiMlZEfhSRNd6vm1xnOhsReUBEVESquc5yOhF5SkTWev8OF4hILdeZTiciz4nIRm/OuSJSxXWm04nIrd6flXwR8ZspgiLS1ftznC4ij7jOUxgReUNE9orIOtdZCiMidUTkCxHZ4P03vqeoNgFR6EXkBjx7zTZT1SuAcY4jnc0/VPUq79c812EKIyJ18GzWvsN1lrN4TlWbqepVwMfAGMd5CvMv4EpVbQZsAh51nKcw64Df4tmjzi+ISDgwHugGNAX6ikhTt6kKNQXo6jrEOeQB96tqE6A1MKKov8eAKPTA74FnVPUEgKrudZwnkP0DeIizbOnomqoeLnAYjR/mVNUFqprnPVyKZ+c0v6KqG1T1e9c5TtMKSFfVLaqaA6Ti6cD5FVX9Esh0neNsVHWPqq7yfn8Ez4a0tc/VJlAK/WVAOxFZJiKLRORq14HOYqT35fwbIuJ3u4iLSE/gR1X9xnWWcxGRv4jITqA//tmjL2gQ8KnrEAGiNrCzwPEuiihQ5txEpD7QAlh2rutKdJ/0CyEi/wZqFHLqMTw54/C8TLkamCUil2gZzw0tIuME4Ck8PdCngOfxFIEyVUTG0UCXsk10pnNlVNUPVPUx4DEReRQYCTxRpgEpOqP3msfwvIyeXpbZTvElo5+RQm7zu1dsgUJEKgFzgHtPeyV8Br8p9Kra6WznROT3wHvewr5cRPLxLOaTUVb54NwZCxKRSXjGl8vc2TKKyK+ABsA3IgKe4YZVItJKVX8qw4g+/z0CM4BPcFDoi8ooIncB3YGOZd3hOKUYf4/+YhdQp8BxArDbUZaAJiKReIr8dFV9r6jrA2Xo5n2gA4CIXAZE4WeryolIzQKHv8HzZpjfUNVvVbW6qtZX1fp4fugSy7rIF0VEGhU47AlsdJXlbESkK/Aw0FNVs1znCSArgEYi0kBEooBk4EPHmQKOeHpqrwMbVPUFn9oEwidjvf8p3gCuAnKAB1T1c6ehTiMib+PJp8A2YJiq7nGZ6VxEZBuQpKr+9gtzDnA5kI9nqerhqvqj21S/JCLpQDlgv/empao63GGkM4jIb4CXgXjgILBGVW90GgrwTjt+EQgH3lDVv7hNdCYReQe4Hs+owX+BJ1T1daehChCRtsBi4Fs8PycAo8810y8gCr0xxpjzFyhDN8YYY86TFXpjjAlyVuiNMSbIWaE3xpggZ4XeGGOCnBV6Y4wJclbojTEmyP0/CMe6BEbKxNcAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Example plot to show the ratio Num / Den\n", + "mu = -2\n", + "sigma = 1\n", + "x1 = np.linspace(norm.ppf(0.0001,mu, sigma), norm.ppf(0.9999,mu, sigma), 100)\n", + "x2 = np.linspace(norm.ppf(norm.cdf(0,mu, sigma),mu, sigma), norm.ppf(0.9999,mu, sigma),100)\n", + "x3 = np.linspace(norm.ppf(norm.cdf(0,mu, sigma),mu, sigma), norm.ppf(0.995,mu, sigma),100)\n", + "y1 = norm.pdf(x1,mu, sigma)\n", + "y2 = norm.pdf(x2,mu, sigma)\n", + "y3 = norm.pdf(x3,mu, sigma)\n", + "print(norm.cdf(0,mu, sigma))\n", + "plt.plot(x1, y1,'k--')\n", + "plt.plot(x2, y2,'b-')\n", + "plt.plot(x3, y3,'r-')\n", + "p1x = np.array([0,0])\n", + "p1y = np.array([0,0.15])\n", + "plt.plot(p1x, p1y, color='r' )\n", + "plt.fill_between(x2,y2,color='b')\n", + "plt.fill_between(x3,y3,color='r')\n" + ] + }, + { + "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.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/book/interactive-nbs/LimitNearBoundariesBayesianApproachPoisson.ipynb b/book/interactive-nbs/LimitNearBoundariesBayesianApproachPoisson.ipynb new file mode 100644 index 0000000..2a8a5a4 --- /dev/null +++ b/book/interactive-nbs/LimitNearBoundariesBayesianApproachPoisson.ipynb @@ -0,0 +1,209 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Compute the bayesian upper limit for a Poisson near the physical boundary" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "95% CL Upper limit on Nobs = 6 with an expected background vb = 3 \n", + "vs_up = 8.91\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "from scipy.stats import poisson\n", + "import matplotlib.pyplot as plt\n", + "from math import exp\n", + "%matplotlib inline\n", + "\n", + "confidenceinterval = 0.95\n", + "beta = 1 - confidenceinterval\n", + "\n", + "# parameters\n", + "nobs = 6\n", + "vb = 3\n", + "\n", + "max_obs = 20 # range to compute the Poisson distribution\n", + "k = np.arange(0, max_obs, 1)\n", + "\n", + "def numerator(nobs, vs_up, vb):\n", + " N = poisson.pmf(k, vs_up + vb)\n", + " sum = 0\n", + " for i in range(nobs+1): # range stops at nobs -1\n", + " sum+=N[i]\n", + " return sum\n", + "\n", + "def denominator(vb):\n", + " N = poisson.pmf(k, vb)\n", + " sum = 0\n", + " for i in range(nobs+1): # range stops at nobs -1\n", + " sum+=N[i]\n", + " return sum\n", + "\n", + "# Compute upper boundary\n", + "vs_up = 0\n", + "s = 0\n", + "scan_step = 0.01\n", + "ck = True\n", + "while (vs_up < 2000) and (ck):\n", + " vs_up = s * scan_step\n", + " # print (\"\\nSignal = \", Signal)\n", + "\n", + " ratio = numerator(nobs, vs_up, vb)/denominator(vb)\n", + " if (ratio < beta):\n", + " ck = False\n", + " # print (Signal, sum)\n", + " s+=1\n", + "\n", + "print (\"95% CL Upper limit on Nobs = \", nobs, \" with an expected background vb = \",vb, \"\\nvs_up = \", vs_up)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "95% CL Upper limit on Nobs = 6\n", + "vb = 0.00 vs_up = 11.90\n", + "vb = 1.00 vs_up = 10.90\n", + "vb = 2.00 vs_up = 9.90\n", + "vb = 3.00 vs_up = 9.00\n", + "vb = 4.00 vs_up = 8.10\n", + "vb = 5.00 vs_up = 7.40\n", + "vb = 6.00 vs_up = 6.80\n", + "vb = 7.00 vs_up = 6.30\n", + "vb = 8.00 vs_up = 5.90\n", + "vb = 9.00 vs_up = 5.60\n", + "vb = 10.00 vs_up = 5.30\n", + "vb = 11.00 vs_up = 5.10\n", + "vb = 12.00 vs_up = 4.90\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfJklEQVR4nO3deZgU1dXH8e+ZYdgEZVQG0AGRIIigMEiUiCCrgnHfiQtGTczixqtxjVlN4pLE+KiJEiWaGFFQcUEDIkhwRXYcQERZFAHRqCAqCHLeP26PM2wydFd199T8Ps9TT+91TgNzprh16x5zd0REJFkKcp2AiIhET8VdRCSBVNxFRBJIxV1EJIFU3EVEEkjFXUQkgXZY3M1suJmtMrPyKs/dYmZvmNkcMxttZk1izVJERHZKdY7c7wMGbvHceKCTux8EvAlcE3FeIiKSgR0Wd3efDHy0xXPPuvvG1MNXgdIYchMRkTTViWAf5wEPb+9FM/sh8EOA+vXrH9yqVasIQlbPpk2bKCjI3mkFxauZsRRP8fI93ptvvvmhuzfdqQ+5+w43oDVQvo3nrwNGA1ad/TRv3tE3bfKsef7557MXTPFqbCzFU7x8jwdM82rU2Kpb2r96zGwIcAxwZir4Dq1cWZ9Bg2Dp0nSjiohIdaRV3M1sIHAVcJy7f17dz5WUrOfFF6FTJ7jzTti0KZ3oIiKyI9WZCjkCeAVob2bLzOx84A6gMTDezGaZ2V3VCdakyZeUl8Nhh8FFF8ERR8CCBRnlLyIi21Cd2TKD3b2Fuxe5e6m73+vubd29pbt3SW0/qm7A1q1h7Fi47z6YOxc6d4Ybb4QNGzL5GiIiUlVOrlA1gyFDYN48OPZYuOYaOPRQmDkzF9mIiCRPTpcfaN4cRo2CRx+F5cvh29+Ga6+FdetymZWISM2XF2vLnHQSzJ8P55wDf/gDdOkCL72U66xERGquvCjuAMXFMHw4jBsXjtx79oSLL4ZPP811ZiIiNU/eFPcKRx4J5eWhsN95Z5g2OW5crrMSEalZ8q64AzRqBLfdBi+8AA0bwsCBcO658NFHO/yoiIiQp8W9Qo8eYQbNddfBAw/AAQeEk68iIvLN8rq4A9SvDzfcANOmwV57wSmnwMknw4oVuc5MRCR/5X1xr9ClC7z2Wrjg6emnw1H8ffdB9Va1ERGpXdLtxHSqmc01s01m1i3eFCvVqQNXXQVz5sCBB8L3vw9HHQVLlmQrAxGRmiHdTkzlwEnA5KgTqo527WDSpDCb5pVXwoya22+Hr77KRTYiIvkn3U5M8909p0t+FRTAT34S1qfp2RMuuQR69QoXQ4mI1HZWnaXYzaw1MMbdO23x/CTgCnef9g2f/boTU9OmTQ8eOXJkJvlukzuMH9+MO+9syxdfFHLOOUs444x3WbfuUxo1ahR5vO1Zu3at4tXAWIqnePker0+fPtPdfeeGwKvT0YPtd2KaBHSrbmeQdu3aRdibZGsrV7qfdpo7uHfu7H733VNjjbelpHeDUScmxVO83MQjm52Y8lGzZvDwwzB6NKxaBT/+8cFcfTV88UWuMxMRya5EFfcKJ5wQxuIHDlzJTTeFaZQvvJDrrEREsietTkxmdqKZLQO+AzxtZnm3+ktxMfzsZwsYPz40AunVC376U1izJteZiYjEr86O3uDug7fz0uiIc4lF//7w+uvw85+H9WqeegruvhsGDcp1ZiIi8UnksMyWdtkFbr0VXn4ZGjeGo48Oa8f/73+5zkxEJB61orhX6N4dZsyA66+HESOgQwcYOVJLGIhI8tSq4g5Qrx785jcwfTq0agWnnx46QS1fnuvMRESiU+uKe4WDDoJXX4VbboGxY8NCZPfeq6N4EUmGWlvcISxEdsUV4YRrly5wwQXhBOyiRbnOTEQkM7W6uFdo2xYmToS77oKpU8OKk3/5ixYiE5GaS8U9paAALrwQ5s2DPn1g6NDQCWru3FxnJiKy81Tct1BaGubC//vf8NZbUFYWTsB++WWuMxMRqT4V920wg+99LywffPLJ8MtfQrduYchGRKQmSLcT0+5mNt7MFqZui+NNMzeaNg3z4Z94Ilzw1L07/Oxn8Pnnuc5MROSbpduJ6WpggrvvB0xIPU6s444LY/EXXAB//CN07hw6QYmI5Ku0OjEBxwP3p+7fD5wQbVr5Z7fdwpo0EyeGufB9+sCPfgSrV+c6MxGRraU75t7M3VcApG5Lokspv/XpExp0X345/P3v0LEjjBmT66xERDaXVps9M/vE3ZtUef1jd9/muHs22uxtT9ytsObPb8wtt7Rn8eJG9O37Pt///mxKS+vGFm9LSW4tluTvpniKt7Oy1mYPWAC0SN1vASyozn7ibrO3pWy0wlq/3v1Xv3IvKnLfZZcNfuut7l9+GXtYd092a7EkfzfFU7ydRRbb7D0JDEndHwI8keZ+ary6dcNUyTlzoEOHNQwdGk64PvdcrjMTkdosrU5MwI3AADNbCAxIPa7V9t8fbr55Dk88AevXw4ABYbXJxYtznZmI1EbVmS0z2N1buHuRu5e6+73u/j937+fu+6Vut5xNUyuZhWmTc+fC738P48aFNeOvvx4++yzX2YlIbaIrVGNQvz5ccw0sWBCucL3hhnBk/9BDWlJYRLJDxT1GpaVhjZoXXghXuw4eDEccAbNm5TozEUk6FfcsOPzwsC7NsGFhvZqDD4Yf/xg+/DDXmYlIUqm4Z0lhIfzgB/Dmm3DRReECqHbt4I47YOPGXGcnIkmj4p5lxcVw220wezZ07QoXXxyWFX7++VxnJiJJouKeIx07wvjx8OijsHYt9O0Lp54KS5fmOjMRSQIV9xwyC3Ph580LDUGefjrMqvnVr7SssIhkRsU9DzRoEObCL1gAJ5wAv/51mB//yCOaOiki6VFxzyMtW4bmIP/9bxibP/VU6NcPXn8915mJSE2j4p6HevWC6dPhb38LJ17LysKJ1490HbCIVFNGxd3MLjWzcjOba2aXRZSTEKZO/uhHsHBhuP3rX8PUybvugq++ynV2IpLv0i7uZtYJ+AFwCNAZOMbM9osqMQl23z3MhZ85Ew48MFz81K1buOpVRGR7Mjly7wC86u6fu/tG4L/AidGkJVs66KDQ4m/UqDA806sX/Pa3HXj33VxnJiL5qFqdmLb5QbMOhHXcvwN8QWiUPc3dL97ifYntxJSreOvWFfDQQ60YMaIlBQVw5plLOe20ZdStuynWuOrEpHiKl5t4sXVi2t4GnA/MACYDdwG3ftP7k9iJKZfxRox42U85xR3c993XffRo902b4ounTkyKp3i5iUcWOzFV/GK41927unsv4CNgYSb7k53TvPl6Ro2CCROgYUM48UQ48shwUZSI1G6ZzpYpSd22Ak4CRkSRlOycvn3DMsK33w7TpoXx+aFD4ZNPcp2ZiORKpvPcHzWzecBTwE/d/eMIcpI01KkTVptcuBAuuCAsTtauHdxzj6ZOitRGmQ7L9HT3A9y9s7tPiCopSd+ee4a58NOnQ/v2YZnhQw+Fl1/OdWYikk26QjWhyspg8mR48EFYuRJ69ICzzoL33st1ZiKSDSruCWYWWvstWADXXRcWImvfHv7wB1i/PtfZiUicVNxrgV12CU26582DAQPg2mvDevJPPaVVJ0WSSsW9FmnTBkaPhmefhbp14bjjYNAgeOONXGcmIlFTca+FBgwIq03eeiu88kpYs+aKK2D16lxnJiJRUXGvpYqK4LLLwtTJc8+FP/85TJ38xz9gU7yrGIhIFqi413IlJfD3v8PUqfCtb8F550H37jBlSq4zE5FMqLgLAAcfDC+9BP/6FyxbFgr8uefCihW5zkxE0qHiLl8zC3PhFyyAq68OLf/atYNbboEvv8x1diKyMzJdW2ZoqgtTuZmNMLP6USUmudO4cZgLX14OvXvDlVdCp04wblwzFXmRGiKTTkx7A5cA3dy9E1AInBFVYpJ7++0X5sI/8wzUqwc33tiBNm3Ckbxm1ojkt0yHZeoADcysDtAQWJ55SpJvBg2COXPgxhvn0L59OJJv2TJMn1QnKJH8lHYnJggNsoHfEToxPevuZ27jPerElLB4b77ZiIcfbsmkSSWYOX36fMDpp79L27ZrI4+VLYqnePkcL6udmIBiYCLQFCgCHgfO+qbPqBNTsuItWeI+dKh7o0ahG1T//u5jx0bTDSrX303xFC+f4pHlTkz9gcXu/oG7bwAeAw7LYH9Sw+yzT7j46Z134MYbYe5cGDgQOneGf/5TM2xEcimT4v4O0N3MGpqZAf2A+dGkJTVJcTFcdRUsXlx5heuQIejkq0gOpV3c3X0K8AihQfbrqX0NiygvqYHq1QsXPr3+ephho5OvIrmTaSemX7r7/u7eyd3PdnetEi6YhRk2EyaEjlDHHAN/+Us4kj/77NDvVUTipStUJVZdu4ZuUG+/DRdfDI8/HrpEDRgQlh7WevIi8VBxl6zY1snXo46CLl3CejY6+SoSLRV3yaotT75+9RWcc45OvopETcVdcmJHJ19XraqX6xRFajQVd8mp7Z18HTy4OyecAGPGwMaNuc5SpOZRcZe8UfXk62mnvcurr8Kxx4bx+p//HBYtynWGIjWHirvknX32gQsvXMS774aG3mVlYQnib30L+vUL68yvW5frLEXym4q75K2iIr4emlm6FH7723Ai9nvfg732gksuCatVisjWVNylRigtDUMzb70Fzz0XplHefXdYx+aQQ2DYMFizJtdZiuSPTJp1tDezWVW2NWZ2WYS5iWyloKByaGb5crjtNvjiC7jwQmjRIjT4fuklXRwlksnaMgvcvYu7dwEOBj4HRkeVmMiO7LFH5dDMlClw5pkwahQcfjgccAD86U+walWusxTJjaiGZfoBb7v70oj2J1JtZpVDMytWwPDhsPvuYb58aSmccgqMHRsumBKpLTLqxPT1TsyGAzPc/Y5tvKZOTIqXk1hLlzbkmWdaMG5cM1avrktJyToGDVrJwIEraN588zXukvxnqXg1P15WOzFVbEBd4EOg2Y7eq05MipeLWOvXu48a5X7UUe5mYTvySPeRI93XrYs+XnUonuLtDLLcianCIMJR+/sR7EskcnXrVg7NLFkCv/wlvPEGnHZaGLa5/HJYvLhhrtMUiVQUxX0wMCKC/YjErlWrUNwXLQrFvk8fuP12OO+8Q+jYEX7xC5g9W7NtpObLqLibWUNgAKF/qkiNUVgY5sqPHAnvvQeXXLKQkhL43e/CMsT77RdWr3ztNRV6qZky7cT0ubvv4e5aqFVqrKZN4cQT3+P558Nsm2HDoG3bsP78oYdC69YwdCi8+GLoDytSE+gKVZEqSkrgBz8IQzarVsH994erYP/2N+jZE/beG37yk7CKpVarlHym4i6yHcXFoZHIk0/CBx+Eq2IPPzwU/P79oXlzOP/8sB69OklJvlFxF6mGxo3hjDPCFbAffACPPQYDB8Ijj8B3vxuGds46K/SI/eKLXGcrouIustMaNoQTT4QHHghDN2PGwMknw3/+E55v2jRMs3z4Yfj001xnK7WVirtIBurVC0fuw4fDypUwfjycfTZMnhyO9Js2heOPh3/+Ez7+ONfZSm1SJ9cJiCRFUVEYi+/fH+64A15+GR59NGxPPgl16oQVLU8+GfbcsyjX6UrCqbiLxKCwMMyu6dkzTKmcOrWy0P/whwA96NIF+vYNBb9nzzCuLxIVDcuIxKygIMyXv/nm0Gxk5kw477zFFBfDnXeGYZ3iYjjsMLj+epg0SW0EJXMq7iJZZBaugD377KVMnBjG4Z97Dq68Mlwg9fvfhyURiothwIDQO/a11zSnXnZeRsMyZtYEuAfoBDhwnru/EkFeIrVCgwZhWKZfv/B49epwMnbCBJg4Ea69Njy/667Qu3flME7HjuEXhcj2ZDrmfhsw1t1PMbO6gJbWE8nAbrvBsceGDeD998MwTUWxf/LJ8HxJSSj0FcW+TZucpSx5Ku3ibma7Ar2AcwHc/UtA1+mJRKhZMzj99LABLF0aivyECWF76KHwfOvWlcW+b9/QT1Zqt7Q7MZlZF2AYMA/oDEwHLnX3z7Z4nzoxKV5ex6qp8dzhnXcaMmNGE2bOLGbmzCasXRumWO6zz2d07foxZWWf0KXLJ5h9UuO+n+JVymonJqAbsBE4NPX4NuC33/QZdWJSvHyMlZR4Gze6T5vmftNNoetUw4buEDpPtW27xs8/3/2OO9xfftn9s88iD7+ZJPx55lM80ujElMmY+zJgmbtPST1+BLg6g/2JSAYKC+Hgg8N25ZVhMbMpU8LwzdNPb+Dxx+Hee8N7Cwpg//2ha9fKrUuXMOYvyZB2cXf3lWb2rpm1d/cFQD/CEI2I5IG6dSsvpOrdew5HHNGbZctgxozKbeLEsEZOhbZtQ6EvK6ss+nvumbvvIOnLdLbMxcC/UzNlFgHfzzwlEYmDGbRsGbbjj698/v33w4VVFQV/6tTQoapCy5abH+F37RpO2GoqZn7LqLi7+yzC2LuI1FDNmoXliwcOrHzu449h1qzNj/KffLKy5WCzZpsf3XftGmbsqODnD60tIyJbKS4OV8r26VP53Nq1oXn4jBmVR/o331x59WyTJpWFvn79Elq0CL1oC3QdfE6ouItItTRqBD16hK3CunVQXr75Ef7tt8P69Qdwww3hM507b36E36FDWEFT4qXiLiJpq18funULW4UNG+Bf/5pKYeG3vy74w4eHog9hDfyDDtp8WOfAA8O+JDoq7iISqaIiaNPmM3r3hiFDwnObNsHChZsP6YwaBcOGhdcLC8N6OVWP8Dt3Dkf+kh4VdxGJXUEBtG8ftsGDw3PuYTmFqkM6zzwD990XXjeDdu02L/hlZeF8gOyYiruI5IRZmGHTujWcdFJ4zh1WrNj8CP+ll2DEiMrP7bvv1nPxZWsq7iKSN8xgr73Cdswxlc9/+OHmc/FnzAhdrSrsumsP2rULhb9ia9Mm3O6zT7igq7ZRcReRvLfnnqF5yYABlc+tXh2mZk6fDpMmrWL9+r2ZNQueeCIsvVDBDEpLt1349903XJCVxOmaKu4iUiPtthv06hW2srKF9O69NxBO3i5fDosWweLFlduiRaHr1fLllRdjQZi907r19ot/TR3jz7QT0xLgU+ArYKPv7JKUIiIRKygIR+qlpaHwb2n9+nAit2rxr7g/ZUq4Oreq3XbbvNjvuy+sXbs7zZuHXwr5OoUziiP3Pu7+YQT7ERGJXb16YRZOu3bbfv2TTzY/4q8o/vPmhdk8oXn5QVydWgO3RYuti3/F4733DtM8c0HDMiIiVTRpEmbilJVt/dqmTbByJTz66AyKi7tuVvwnT4YHHwzvqVBUBK1abb/477FHfOvxpN2JCcDMFgMfE5pj3+3uw7bxHnViUry8jqV4ihdVvA0bjFWr6rFiRQNWrKjPypX1v76/YkV9Vq/efNpOgwYbadFiXWr7gubN17HXXuto3jzcb9Ag/KbIaiem1C+FvVK3JcBsoNc3vV+dmBQvH2MpnuJlK96aNe5z5rg//rj7rbe6X3KJ+7HHunfsWNk5q+pWUuLevXv2OzHh7stTt6vMbDRwCDA5k32KiCRV48ZhHZ0DD9z6NXf44IPNT/BWbOlIu7ib2S5Agbt/mrp/JPCbdPcnIlKbmUFJSdgOPXTr13ZWJkfuzYDRFqLWAR5097EZ7E9ERCKSSQ/VRUDnCHMREZGIJPCiWxERUXEXEUkgFXcRkQRScRcRSSAVdxGRBFJxFxFJIBV3EZEEUnEXEUkgFXcRkQTKuLibWaGZzTSzMVEkJCIimYviyP1SYH4E+xERkYhkVNzNrBT4LnBPNOmIiEgUMu3E9AjwB6AxcIW7H7ON96gTk+LldSzFU7x8j5fVTkzAMcBfU/d7A2N29Bl1YlK8fIyleIqX7/FIoxNTJsMyPYDjzGwJ8BDQ18weyGB/IiISkbSLu7tf4+6l7t4aOAOY6O5nRZaZiIikTfPcRUQSKKMG2RXcfRIwKYp9iYhI5nTkLiKSQCruIiIJpOIuIpJAKu4iIgmk4i4ikkAq7iIiCaTiLiKSQCruIiIJpOIuIpJAaRd3M6tvZq+Z2Wwzm2tmv44yMRERSV8myw+sB/q6+1ozKwJeNLP/uPurEeUmIiJpSru4p9YYXpt6WJTa0u/8ISIikcm0E1MhMB1oC9zp7ldt4z3qxKR4eR1L8RQv3+NltRNT1Q1oAjwPdPqm96kTk+LlYyzFU7x8j0eWOzFV/QXxCWHJ34FR7E9ERDKTyWyZpmbWJHW/AdAfeCOivEREJAOZzJZpAdyfGncvAEa6+5ho0hIRkUxkMltmDlAWYS4iIhIRXaEqIpJAKu4iIgmk4i4ikkAq7iIiCaTiLiKSQCruIiIJpOIuIpJAKu4iIgmk4i4ikkCZrC3T0syeN7P5qU5Ml0aZmIiIpC+TtWU2Ape7+wwzawxMN7Px7j4votxERCRNaR+5u/sKd5+Ruv8pMB/YO6rEREQkfZGMuZtZa8IiYlOi2J+IiGQmozZ7AGbWCPgv8Dt3f2wbr6vNnuLldSzFU7x8j5f1NnuEptjjgP+rzvvVZk/x8jGW4ilevscjm232zMyAe4H57v7ndPcjIiLRy2TMvQdwNtDXzGaltqMjyktERDKQSSemFwGLMBcREYmIrlAVEUkgFXcRkQRScRcRSSAVdxGRBFJxFxFJIBV3EZEEUnEXEUkgFXcRkQRScRcRSSAVdxGRBMqouJvZcDNbZWblUSUkIiKZy/TI/T5gYAR5iIhIhDIq7u4+GfgoolxERCQiUXRiag2McfdO23ldnZgUL69jKZ7i5Xu8rHdiSv1iaA2UV+e96sSkePkYS/EUL9/jkc1OTCIikr9U3EVEEijTqZAjgFeA9ma2zMzOjyYtERHJRNpt9gDcfXBUiYiISHQ0LCMikkAq7iIiCaTiLiKSQCruIiIJpOIuIpJAKu4iIgmk4i4ikkAq7iIiCaTiLiKSQJkuPzDQzBaY2VtmdnVUSYmISGbSLu5mVgjcCQwCDgAGm9kBUSUmIiLpy+TI/RDgLXdf5O5fAg8Bx0eTloiIZCKThcP2Bt6t8ngZcOiWb6raiQlYn+Vm2nsCHypejYyX5O+meIq3s9rv7AcyKe62jee26tnn7sOAYQBmNs13tlVUBhSv5sZL8ndTPMVLJ97OfiaTYZllQMsqj0uB5RnsT0REIpJJcZ8K7Gdm+5pZXeAM4Mlo0hIRkUykPSzj7hvN7CJgHFAIDHf3uTv42LB046VJ8WpuvCR/N8VTvNjjWWisLSIiSaIrVEVEEkjFXUQkgbJS3LO9TIGZDTezVdmYU29mLc3seTObb2ZzzezSmOPVN7PXzGx2Kt6v44xXJW6hmc00szFZiLXEzF43s1npTAFLI14TM3vEzN5I/T1+J8ZY7VPfq2JbY2aXxRUvFXNo6t9KuZmNMLP6Mca6NBVnblzfa1s/32a2u5mNN7OFqdvimOOdmvqOm8wssimR24l1S+rf5hwzG21mTaq1M3ePdSOcbH0baAPUBWYDB8QcsxfQFSjPwvdrAXRN3W8MvBnn9yNcX9Aodb8ImAJ0z8L3/D/gQWBMFmItAfaMO06VePcDF6Tu1wWaZCluIbAS2CfGGHsDi4EGqccjgXNjitUJKAcaEiZrPAfsF0OcrX6+gZuBq1P3rwZuijleB8KFRZOAbjHHOhKok7p/U3W/WzaO3LO+TIG7TwY+ijNGlVgr3H1G6v6nwHzCD1Rc8dzd16YeFqW2WM+Km1kp8F3gnjjj5IKZ7Ur4gboXwN2/dPdPshS+H/C2uy+NOU4doIGZ1SEU3riuR+kAvOrun7v7RuC/wIlRB9nOz/fxhF/SpG5PiDOeu8939wVRxdhBrGdTf54ArxKuKdqhbBT3bS1TEFvxyyUzaw2UEY6m44xTaGazgFXAeHePNR7wF+BKYFPMcSo48KyZTU8tXxGnNsAHwD9Sw073mNkuMcescAYwIs4A7v4e8EfgHWAFsNrdn40pXDnQy8z2MLOGwNFsfqFjnJq5+woIB1xASZbiZtt5wH+q88ZsFPdqLVNQ05lZI+BR4DJ3XxNnLHf/yt27EH6DH2JmneKKZWbHAKvcfXpcMbahh7t3Jaw4+lMz6xVjrDqE/wb/zd3LgM8I/62PVerCv+OAUTHHKSYc1e4L7AXsYmZnxRHL3ecThg3GA2MJQ7Abv/FDUm1mdh3hz/Pf1Xl/Nop74pcpMLMiQmH/t7s/lq24qeGDScDAGMP0AI4zsyWEIbW+ZvZAjPFw9+Wp21XAaMLQXlyWAcuq/O/nEUKxj9sgYIa7vx9znP7AYnf/wN03AI8Bh8UVzN3vdfeu7t6LMLywMK5YW3jfzFoApG5XZSluVpjZEOAY4ExPDb7vSDaKe6KXKTAzI4zXznf3P2chXtOKs+Vm1oDww/tGXPHc/Rp3L3X31oS/u4nuHsuRH4CZ7WJmjSvuE04mxTbryd1XAu+aWcWqe/2AeXHFq2IwMQ/JpLwDdDezhql/q/0I54ViYWYlqdtWwElk5ztCqClDUveHAE9kKW7szGwgcBVwnLt/Xu0PRnkW+xvOAB9NmEXyNnBdFuKNIIwvbiAcmZ0fY6zDCcNMc4BZqe3oGOMdBMxMxSsHfpGNv8NU7N7EPFuGMAY+O7XNzdK/ly7AtNSf6eNAcczxGgL/A3bL0t/brwkHAOXAv4B6McZ6gfDLcTbQL6YYW/18A3sAEwj/U5gA7B5zvBNT99cD7wPjYoz1FuG8ZUV9uas6+9LyAyIiCaQrVEVEEkjFXUQkgVTcRUQSSMVdRCSBVNxFRBJIxV1EJIFU3EVEEuj/AfPMxHVxQhbvAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import numpy as np\n", + "from scipy.stats import poisson\n", + "import matplotlib.pyplot as plt\n", + "from math import exp\n", + "%matplotlib inline\n", + "\n", + "confidenceinterval = 0.95\n", + "beta = 1 - confidenceinterval\n", + "\n", + "# parameters\n", + "nobs = 6\n", + "vb_min = 0\n", + "vb_max = 12\n", + "\n", + "max_obs = 20 # range to compute the Poisson distribution\n", + "k = np.arange(0, max_obs, 1)\n", + "\n", + "def numerator(nobs, vs_up, vb):\n", + " N = poisson.pmf(k, vs_up + vb)\n", + " sum = 0\n", + " for i in range(nobs+1): # range stops at nobs -1\n", + " sum+=N[i]\n", + " return sum\n", + "\n", + "def denominator(vb):\n", + " N = poisson.pmf(k, vb)\n", + " sum = 0\n", + " for i in range(nobs+1): # range stops at nobs -1\n", + " sum+=N[i]\n", + " return sum\n", + "\n", + "# Compute upper boundaries\n", + "\n", + "v_vs_up = []\n", + "v_vb = []\n", + "\n", + "print (\"95% CL Upper limit on Nobs = \", nobs)\n", + "\n", + "nsteps = 12\n", + "scan_step = 0.1\n", + "for i in range(nsteps+1): \n", + " vb = vb_min + i/nsteps*(vb_max-vb_min)\n", + " \n", + " vs_up = 0\n", + " s=0\n", + " ck = True\n", + " while (vs_up < 2000) and (ck):\n", + " vs_up = s * scan_step\n", + " ratio = numerator(nobs, vs_up, vb)/denominator(vb)\n", + " if (ratio < beta):\n", + " ck = False\n", + " s+=1\n", + " print (\"vb = {:.2f}\".format(vb), \" vs_up = {:.2f}\".format(vs_up))\n", + " v_vs_up.append(vs_up)\n", + " v_vb.append(vb)\n", + "\n", + "plt.plot(v_vb, v_vs_up, 'b-')\n", + "plt.axis([0,12,0,12])\n", + "plt.xticks(np.arange(0, 13, 1.0))\n", + "plt.yticks(np.arange(0, 13, 1.0))\n", + "\n", + "plt.grid()\n", + "plt.show()\n", + "\n" + ] + }, + { + "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.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}