Thirlwall_probability_calculation.py

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import scipy as sp
import numpy as np


# In[2]:


# https://thirlwall.public-inquiry.uk/wp-content/uploads/thirlwall-evidence/INQ0108786.pdf
# (See last page)


# In[3]:


# From paragraph 11.3
obs = np.array([1,3,3,2,3])
mu = obs.mean()

# Be sure to use standard error equation for Poisson distribution:
# https://onbiostatistics.blogspot.com/2014/03/computing-confidence-interval-for.html?m=1
se = np.sqrt(obs.mean()/len(obs))


# In[4]:


# Matches values in paragraph 11.5
print(mu, se)


# se - mu < 0, so this case is underdispersed... the Poisson-Gamma/negative binomial is mainly used for overdispersion.  
# 
# Atypical use but not incorrect I think, as long as the rate in the Poisson takes values >0 I don't see why we can't choose the mean and se of the Gamma distribution it is drawn from however we like.

# In[5]:


# Determine Gamma distribution parameters to match observed rate and se

# https://math.stackexchange.com/a/1810274
# Shape
a = (mu/se)**2
# Rate
b = mu/se**2

print(a,b)


# In[6]:


mu_gam, var_gam = sp.stats.gamma.stats(a,scale=1/b, moments='mv')


# In[7]:


# Verify that Gamma parameters were set correctly by ensuring the Gamma distribution has first two moments which match the ovserved ones. 
print(mu_gam, np.sqrt(var_gam))


# In[8]:


assert mu == mu_gam
assert se == np.sqrt(var_gam)


# In[9]:


# The negative binomial model is equivalent to the gamma-Poisson for correct choices of the their parameters.
# This answer shows us how to set the parameters of nbinom so they correspond to particular shape and rate values of the gamma-Poisson model:
# https://stackoverflow.com/a/77746724


# In[10]:


# Use the survival function to compute the probability of 8 or more deaths


# In[11]:


sp.stats.nbinom.sf(7, a, b/(1+b))


# In[12]:


# Rounded up this indeed matches the 0.008 quoted in paragraph 11.5! :)


# In[13]:


# Just verifying that the survival function is doing what we expect it to so we are not misusing it, it is:


# In[14]:


1-sp.stats.nbinom.cdf(7, a, b/(1+b))