Inversion Method

Similar to generating continuous random numbers, the inversion method is also the simplest and most straightforward approach for generating discrete random numbers.

Let \(X\) be a discrete random variable assuming values \(x_1,x_2, \cdots\) with respective probabilities \(p_i\). Also let us denote the respective distribution function as \(F_i=P(x\leq x_i)=p_1+p_2+\cdots+p_i\) Then we can have \[X=\begin{cases}x_1 & \text{if } u \leq p_1 \\x_2& \text{if } u \leq p_1+p_2 \\ x_3 & \text{if } u \leq p_1 +p_2+p_3\\ \vdots &\\x_n & \text{if } u \leq p_1+p_2+\cdots+p_n \\ \vdots & \end{cases}\] In in terms of sitribution function it can be summarised as follows \[X=\begin{cases}x_1 & \text{if } u \leq F_1 \\x_2& \text{if } u \leq F_2 \\ x_3 & \text{if } u \leq F_3\\ \vdots &\\x_n & \text{if } u \leq F_n \\ \vdots & \end{cases}\]

Now to show that the generated random number has the same distribution as the original density function, consider

\[P(X=x_i)=P(F_{i-1}<u \leq F_i)=F_i-F_{i-1}=p_i\]

Therefore, the inversion method can be algorithmically expressed as follows.

  1. Generate uniform random numbers

  2. Construct the cumulative distribution function (CDF) for the given probability distribution.

  3. Now if the random number corresponding to the given uniform random number is \(x_j\) if \(F(x_{j-1})<u\leq F(x_j)\) where \(x_i\)’s are ordered realizations of the given random variable \(X\).

Example 1

Generate five random numbers of a random variable with density function given below

X 1.00 5.00 8.00 12.00 15.00 20.00
pmf 0.12 0.19 0.29 0.21 0.11 0.08 
import random

samples=6

values=[1,5,8,12,15,20]
pmf=[0.12,0.19,0.29,0.21,0.11,0.08]
df=[0]*len(pmf)
for i in range(0,len(pmf)):
  df[i]=sum(pmf[:i+1])
print(df)
## [0.12, 0.31, 0.6, 0.8099999999999999, 0.9199999999999999, 0.9999999999999999]
# Generate uniform random numbers
u=[random.random() for i in range(samples)]
rnos=[]
for s in u:
  for i,t in enumerate(df):
    if s<t:
      rnos.append(values[i])
      break
print(u)
## [0.5033873789936537, 0.42478731283826365, 0.7785532009351623, 0.38086448654769023, 0.3395025352197536, 0.5878342578161238]
print(rnos)
## [8, 8, 12, 8, 8, 8]
Exercise 1

Generate 10 random numbers from a Bernoulli distribution with probability of success as \(0.4\)

Example 2: Generating Uniform Random Numbers

Suppose \(X\) is a uniform discrete random variable assuming values \(1,2,\cdots,n\). Then the distribution function is given by \(F(x_j)=P(X\leq x_j)=\dfrac{j}{n}\). Hence therefore using the inversion method \[X=x_j \text{ if } \begin{cases}\dfrac{j-1}{n} <u \leq \dfrac{j}{n} \\ j-1 <nu\leq j\end{cases}\].Therefore \(X=[nu]+1\) Generate 10 uniform random numbers with the possible values \(1,2,3,4,5\)

import random
u=[random.random() for i in range(10)]
rnos=[int(5*i)+1 for i in u]
print(rnos)
## [4, 4, 4, 4, 4, 2, 4, 2, 1, 2]

Example 3: Generating Geometric Random Numbers

Consider \(G(x) =1-e^{-\lambda x}\) \[\begin{eqnarray*} G(x+1)-G(x)&=& \left[ 1- e^{- \lambda (x+1)}\right]- \left[1- e^{- \lambda x}\right] = e^{-\lambda x} - e^{-\lambda (x+1)}\\ &=& e^{-\lambda x} \left[1- e^{-\lambda}\right] \end{eqnarray*}\] The RHS can be written in the form \((1-p)^x\;p\) where \(p = 1- e^{-\lambda}\). Hence the RHS is the pdf of geometric distribution. This implies that \(G(x)\), for integer values of \(x\), is the distribution function of geometric random variable with parameter \(1- e^{-\lambda}\) But as the distribution function lies between \(0\) and \(1\), we can generate its values with \(u\), standard uniform random number. Now \(u = 1-e^{-\lambda x} \implies e^{-\lambda x} = 1-u \implies x=\frac{-\log u}{\lambda}\)

Since geometric random variable can be integers only and we are dealing with distribution function, we round the value to the integer just above. Hence \(\lceil\frac{ -\log u }{\lambda}\rceil\) follows geometric distribution with parameter \(p=1-e^{-\lambda}\).

Now \(p=1-e^{-\lambda} \implies \lambda = -\log p\). Hence \(\lceil\frac{ \log u }{\log (1-p)}\rceil\) follows geometric distribution with parameter \(p\)

import random
import math
import numpy as np

no_sam=10
p=0.4
u=[random.random() for i in range(no_sam)]
x=[math.ceil(np.log(i)/np.log(1-p)) for i in u]
x
## [4, 2, 3, 1, 5, 1, 1, 1, 7, 4]
import random
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
no_sam=1000
p=0.3
u=[random.random() for i in range(no_sam)]
x=[math.ceil(np.log(i)/np.log(1-p)) for i in u]
#Now to generate the frequency table we have to convert the given data into dataframe. 
data=pd.DataFrame(x)
#change the column label as Random Numbers
data.columns=['Random Numbers']
#value_counts generates frequency table with values as index and frequencies and the value and then it is sorted by values. 
freq=data['Random Numbers'].value_counts().sort_index()/no_sam
plt.bar(freq.index,freq.values,width=0.05)
## <BarContainer object of 18 artists>
plt.xlabel('Values')
plt.ylabel('Empirical Probability ')
plt.title('Empirical probability Plot of Random Numbers')

vals=[i for i in range(1,data['Random Numbers'].max()+3)]
pmf=[(1-p)**(i-1)*p for i in vals]
plt.plot(vals,pmf,marker='*',linestyle='none', markersize=8, label='Pmf',color='r')
# Show the plot
plt.show()

Generating Poisson Random Numbers

A method for generating Poisson random numbers relies on the relationship that connects \(f(x+1)\) to \(f(x)\), which can be expressed as \(f(x+1) = \frac{\lambda}{x+1} f(x)\). We will use the above recurrence relationa and the inverse transformation to obtain the random numbers.

The algorithm for generating Poisson random numbers with mean \(\lambda\) is given as follows (Ross,2022). The quantity \(i\) refers to the value presently under consideration; \(p = p_i\) is the probability that \(X\) equals \(i\), and \(F = F (i)\) is the probability that \(X\) is less than or equal to \(i\).

  1. Generate uniform random number \(u\)

  2. For \(i=0, p=e^{-\lambda}, F=p\)

  3. If \(u<F\), set the random number as \(i\) and stop

  4. \(p=\dfrac{\lambda\times p}{i+1}, F=F+p, i=i+1\)

  5. Goto step 3

import random
import numpy as np

n_samples=3
mean=1
u=[random.random() for i in range(n_samples)]
print(u)
## [0.9447168833283385, 0.5705148395420875, 0.23332530932637763]
i=0
j=0
p=np.exp((-1)*mean)
F=p
print(F)
## 0.36787944117144233
for s in u:
  while (s>F):
    i=j
    j=j+1
    p=mean*p/j
    F=F+p
    print("Value of the distribution function for the value {} is {:.5f}".format(j,F))
  print("Poisson random number corresponidng to uniform random number {:.5f} is {}".format(s,i))
## Value of the distribution function for the value 1 is 0.73576
## Value of the distribution function for the value 2 is 0.91970
## Value of the distribution function for the value 3 is 0.98101
## Poisson random number corresponidng to uniform random number 0.94472 is 2
## Poisson random number corresponidng to uniform random number 0.57051 is 2
## Poisson random number corresponidng to uniform random number 0.23333 is 2

Generating Poisson Random Numbers using Exponential random numbers

The generation of Poisson random numbers using exponential random numbers is rooted in a fundamental property of the Poisson process. This property states that the count of exponential random numbers, each with a mean of \(\frac{1}{\lambda}\), whose sum exceeds one, follows a Poisson distribution with a mean of \(\lambda\).

To generate Poisson random numbers using exponential random numbers in Python, you can use the following approach:

  1. Generate exponential random numbers with mean \(\frac{1}{\lambda}\).

  2. Sum the exponential random numbers until the sum exceeds 1.

  3. The count of exponential random numbers generated until the sum exceeds 1 follows a Poisson distribution.

The following

import random
import numpy as np
nsamples=3
P_mean=2
P_rnos=[]
def expo_rn(mean):
  u=random.random()
  exp=(-1)*mean*np.log(u)
  return(exp)
for i in range(nsamples):
  i=0
  x=expo_rn(1/P_mean)
  while (x<1):
    i=i+1
    x=x+expo_rn(1/P_mean)
  P_rnos.append(i)
P_rnos
## [0, 1, 3]

Exercise