Is there an inverse probability mass function for Poisson statistics in Python?

Question

I am looking for an inverse pmf for Poisson statistics. By inverse I mean a function inv_pmf(p, k) that returns the distribution parameter lambda. For clarity the parameters are used as follows: p = lambda^k / k! * exp(-lambda). Thanks

Answer 1

So you have equation for probability p(k,λ) = λ^k e^-λ/k!. And you know p and k but want to know λ. Well, take log from lhs and rhs and get simple equation.

log(p) = k*log(λ) - λ - log(k!)

or

λ = k*log(λ) - log(p) - log(G(k+1)), where G() is Gamma-function, which is available in Python lib. You could plot difference between RHS and LHS and see that it could have multiple solutions. Them using Python fsolve function, you could get root out of that non-linear equation.

Code (Python 3.7, Anaconda, Windows 10 x64)

#%%    
import math
import numpy as np
import matplotlib.pyplot as plt

def logrhs(p, k, λ):
    return k*np.log(λ) - math.log(p) - math.lgamma(k+1)

def poissonPMF(k, λ):
    lp = k*np.log(λ) - λ - math.lgamma(k+1)
    return np.exp(lp)

p = 0.2
k = 3.0

λλ = np.linspace(0.001, 10.0, 101)

#%%    
rhs = logrhs(p, k, λλ)
lhs = np.copy(λλ)

pmf = poissonPMF(k, λλ)

plt.plot(λλ, lhs - rhs, 'r')
plt.plot(λλ, pmf, 'g')
plt.show()

# %%
from scipy.optimize import fsolve

def f(x):
    return x - logrhs(p, k, x)

starting_guess = 4.0
λ = fsolve(f, starting_guess)
print((λ, poissonPMF(k, λ)))

starting_guess = 1.9
λ = fsolve(f, starting_guess)
print((λ, poissonPMF(k, λ)))

As an example, I use k=3 and p=0.2 for a test. For me, it prints two roots

(array([3.90263215]), array([0.2]))
(array([2.24859448]), array([0.2]))

and verifies that indeed probability is equal to 0.2. GRaph shows red line is crossing 0 in two places