DDC-6: A Poisson counter

A data challenge a day helps you master machine learning

About these daily data challenges

Each post is an exercise that helps you learn about data in Python.

Try to solve the exercise before checking my solution at the bottom of the post 🤓

You can share your solution or visualization in the comments!

Today’s challenge

A Poisson random variable is the number of iterations it takes for a randomly decreasing number to be smaller than a target value. It is used in physics and computational biology to estimate the duration of different events.

Part 1. Write a Python function that calculates a Poisson random number. Translate the following algorithm into code.

1. Initialize a counter to zero.

2. Initialize a “current value” numerical variable to 1.

3. Initialize a target variable as exp(-λ) (λ is a parameter).

4. In a while-loop, generate a random number from a standard uniform distribution (that is, equal probability of a number between 0 and 1) and scale that by the current-value variable. Increment the counter by 1.

5. Remain in the while-loop until the current-value is less than the target.

6. Return the counter value. This is the random Poisson variable.

Check your code using λ=10. The result will be an integer that will differ each time you call the function, but the counter should be close to 10.

Part 2. Call the function 500 times using λ=10, and show a bar plot of the distribution of random Poisson variables. An interesting feature of a Poisson distribution is that the mean and variance equal each other, which also equals λ. In finite sample sizes, the mean and variance won’t equal exactly λ, but it should be fairly close.

.

.

.

.

Scroll down for the solution…

.

.

.

.

a fish doing mathematics

.

.

.

.

keep scrolling!

.

.

.

.

import numpy as np
import matplotlib as mpl

def poissonCounter(lam):

  # initialize
  counter,currval = 0,1
  target = np.exp(-lam)

  # run algorithm
  while currval>target:
    currval *= np.random.rand()
    counter += 1
  
  # return result
  return counter-1 # -1 to find value<target

nIterations = 500
lambda_param = 10
poissonRand = np.zeros(nIterations,dtype=int)

for i in range(nIterations):
  poissonRand[i] = poissonCounter(lambda_param)

# descriptive characteristics
mean = poissonRand.mean()
var = poissonRand.var()

# get histogram values
y,x = np.histogram(poissonRand,np.arange(poissonRand.min(),poissonRand.max()+1))

# show as bar plot
plt.figure(figsize=(8,3))
plt.bar(x[:-1],y,width=x[1]-x[0],facecolor=[.7,.9,.7],edgecolor='k')
plt.gca().set(xlabel='Value',ylabel='Count',xticks=x,
              title=f'Distribution of Poisson random variables ($\\lambda={lambda_param}$)\nmean = {mean:.2f}, var. = {var:.2f}')

plt.show()