The following assignment is to generate the random numbers from various discrete distributions using the algorithm mentioned in the book Simulation by Sheldon M. Ross.

Binomial distribution

Generating random numbers from \(X \sim {\sf Binom}(n=7,p_1=0.4)\)using the algorithm given in simulation book

n=7 # Binomial parameter
p=0.4
set.seed(100)
n1<- 1000 # sample size and not binomial parameter
U<-runif(n1)
# Finding distribution function values for values of X
c=p/(1-p)
X<-c(0:7)
pr<-c()
pr[1]=(1-p)^n
F<-c()
F[1]=pr[1]
for(i in 2:(n+1)){
  pr[i]=((c*(n-X[i-1])/X[i]))*pr[i-1]
  F[i]<-F[i-1]+pr[i]
}
# distribution function of X
F
## [1] 0.0279936 0.1586304 0.4199040 0.7102080 0.9037440 0.9811584 0.9983616
## [8] 1.0000000

Generating samples of size 1000 from \({\sf Binomial}(n=7,p_1=0.4)\) distribution

x<-c()
for(i in 1:n1){
  x[i]<-min(which(U[i]<F))-1
}
# Frequency distribution of the generated distribution 
table(x)
## x
##   0   1   2   3   4   5   6   7 
##  29 113 262 284 202  86  22   2
# Plot
barplot(table(x)/n1)

## Poisson distribution
Generating random numbers from Poisson distribution with parameters \(\lambda = 3\) using the algorithm given in simulation book.

lambda= 3 # parameter of the distribution 
n=1000 # sample size
set.seed(1001)
U<-runif(n) # sample from uniform distribution

# Finding distribution function values for values of X
X<-c(0:20)
p<-c()
p[1]=exp(-lambda)
F<-c()
F[1]<-p[1]
for(i in 2:length(X)){
  p[i]=(lambda*p[i-1])/X[i]
  F[i]<-F[i-1]+p[i]
}
# distribution function of X
F
##  [1] 0.04978707 0.19914827 0.42319008 0.64723189 0.81526324 0.91608206
##  [7] 0.96649146 0.98809550 0.99619701 0.99889751 0.99970766 0.99992861
## [13] 0.99998385 0.99999660 0.99999933 0.99999988 0.99999998 1.00000000
## [19] 1.00000000 1.00000000 1.00000000

Generating samples of size 1000 from Poisson(\(\lambda=3\)) distribution.

x<-c()
for(i in 1:n){
  x[i]<-min(which(U[i]<F))-1
}
# Frequency distribution of the generated distribution 
table(x)
## x
##   0   1   2   3   4   5   6   7   8   9  10  11 
##  54 156 225 227 155 103  46  25   5   2   1   1
# Plot
barplot(table(x)/n)

Hypergeometric distribution

Generating random numbers from Hypergeometric distribution with parameters m = 15(the number of white balls in the urn), n = 10(the number of black balls in the urn), k= 7(the number of balls drawn from the urn)

n1<-1000 # sample size
set.seed(1002)
U<-runif(n1)
#defining the parameters as defined above
m=15;n=10;k=7
# range of X is 0 to 7 as max(0,k-n)<=x<=min(k,m)
max(0,k-n)
## [1] 0
min(k,m)
## [1] 7
# distribution function of X
F<- phyper(0:k,m,n,k)
F
## [1] 0.0002496359 0.0068025796 0.0618473060 0.2606199293 0.6013729977
## [6] 0.8824942792 0.9866132723 1.0000000000

Generating samples of size 1000 from Hypergeometric distribution with the given parameters

x<-c()
for(i in 1:n1){
  x[i]<-min(which(U[i]<F))-1
}
table(x)
## x
##   1   2   3   4   5   6   7 
##   6  53 186 344 294 109   8
barplot(table(x)/n1)

Geometric Distribution

Generating random numbers from Geometric distributions Consider p=0.3, we are using the following formula to get the random numbers \[X={\frac{\log(U)}{\log(q)}}\]
#### Note: In the above expression integer part of the ratio is missing, while in code it is taken care, kindly consider.

p=0.3
q=1-p
# Generating random numbers from Uniform distribution
n=100
set.seed(100)
U<-runif(n)
# Generating numbers from Geometric distribution 
X=ceiling(log(U)/log(q))
hist(X,freq = FALSE)

Negative Binomial distribution

Generating random numbers from Negative Binomial distribution with parameters n=5(target number of successful trials) prob=0.7 (probability of success in each trail)

n1<-1000 # sample size
set.seed(1002)
U<-runif(n1)
#defining the parameters as defined above
n=5;prob=0.7

# distribution function of X
F<- pnbinom(0:20,size=n,prob)
F
##  [1] 0.1680700 0.4201750 0.6470695 0.8058956 0.9011913 0.9526510 0.9783808
##  [8] 0.9905106 0.9959690 0.9983343 0.9993278 0.9997342 0.9998967 0.9999605
## [15] 0.9999851 0.9999944 0.9999980 0.9999993 0.9999997 0.9999999 1.0000000

Generating samples of size 1000 from Negative Binomial distribution with the given parameters

x<-c()
for(i in 1:n1){
  x[i]<-min(which(U[i]<F))-1
}
table(x)
## x
##   0   1   2   3   4   5   6   7   8   9  10 
## 150 259 232 155 107  59  22  10   2   3   1
barplot(table(x)/n1)

Multinomial distribution

Generating random vector of size 1000 from \((X_1,X_2,X_3) \sim {\sf Multinom}(n=10,p_1=0.25,p_2=0.45,p_3=1-p_1-p_2=0.3)\)

# defining the parameters of the distribution
n=10 # number of trials
p1= 0.25 # probability of getting outcome of class 1 
p2 = 0.45 # probability of getting outcome of class 2
p3 = 1- p1-p2 # probability of getting outcome of class 3

As the output is having three classess we need to generate two random variables i.e \(X_1\), \(X_2\).
\(X_3\) can be obtained by the following equation \(X_3 = n-X_1-X_2\).
So first we will generate observations from \(X_1 \sim {\sf Binom}(n,p_1)\) using the algorithm above for Binomial distribution.

set.seed(103)
n1<- 1000 # sample size and not binomial parameter
U<-runif(n1)
# Finding distribution function values for values of X1
c=p1/(1-p1)
X_1<-c(0:10)
pr_1<-c()
pr_1[1]<-(1-p1)^n
F_1<-c()
F_1[1]=pr_1[1]
for(i in 2:(n+1)){
  pr_1[i]=((c*(n-X_1[i-1])/X_1[i]))*pr_1[i-1]
  F_1[i]<-F_1[i-1]+pr_1[i]
}
# distribution function of X1
F_1
##  [1] 0.05631351 0.24402523 0.52559280 0.77587509 0.92187309 0.98027229
##  [7] 0.99649429 0.99958420 0.99997044 0.99999905 1.00000000

Generating samples of size 1000 from \(X_1 \sim {\sf Binom}(n=10,p_1=0.25)\) distribution

x1<-c()
for(i in 1:n1){
  x1[i]<-min(which(U[i]<F_1))-1
}
# Frequency distribution of the generated distribution 
table(x1)
## x1
##   0   1   2   3   4   5   6   7 
##  38 207 275 233 166  70   9   2
# Plot
barplot(table(x1)/n1)

Once the observations from \(X_1\) are generated, we will now generate observations from \(X_2|X_1=x1 \sim {\sf Binom}(n-x_1,\frac{p_2}{1-p_1})\) using the above algorithm. As we note in the above equation that value of the size parameter depends on the value of the generated sample \(x_1\), hence we are using inbuilt distribution function for getting the distribution function of \(X_2|X_1=x_1\)

x2<-c()
for(i in 1:n1){
  F_2<-pbinom(0:(n-x1[i]),size=n-x1[i],prob = (p2/(1-p1))) # distribution function of X2 given X1 =x1
  x2[i]<-min(which(U[i]<F_2))-1 # sample which is getting generated 
}
# Frequency distribution of the generated distribution 
table(x2)
## x2
##   1   2   3   4   5 
##   1   7  75 385 532

Once the observations from \(X_1 \& X_2\) are generated, we will now generate observations from \(X_3|X_1=x1,X_2=x2 \sim {\sf Binom}(n-x_1-x_2,\frac{p_3}{1-p_1-p_2})\) using the above algorithm.
But if we notice \(p_3 = 1-p_1-p_2\) therefore \(X_3|X_1=x1,X_2=x2 \sim {\sf Binom}(n-x_1-x_2,1)\) So we can use the fact that \(\sum_{i = 1}^{r} x_i = n\) and find \(x_3 = n-x_1-x_2\).
Therefore, the sample from \((X_1,X_2,X_3) \sim {\sf Multinom}(n,p_1,p_2,p_3)\) is generated and is given below.

xmult_sample<-cbind(x1,x2,x3=n-x1-x2)
dim(xmult_sample)
## [1] 1000    3
head(xmult_sample)
##      x1 x2 x3
## [1,]  1  4  5
## [2,]  1  3  6
## [3,]  2  5  3
## [4,]  2  5  3
## [5,]  1  4  5
## [6,]  1  3  6