Homework 3

Question 1. Consider the Binomial distribution with \(n = 24\) and \(p = .9\) that is used to model the number of correctly received bits on a satellite link that transmits data in \(24\)-bit blocks.

Quesion 1 (a) Plot the probability density function and cumulative distribution function for the number of correctly received bits

#Answer 1a
n <- 24
x <- 0:n
p <- 0.9
#PDF
plot (x,dbinom(x,n,p), type = 'h')

xlab = "Number of correctly-received bits"
ylab= "probability"

#CDF
cdf<-pbinom (0:n, 24, 0.9)
plot(x,cdf, type ='s')

xlab= "Number of correctly-received bits" 
ylab= "CDF"

Quesion 1 (b) What is the mean number of correctly received bits? What is the standard deviation?

#Answer 1b
#(make sure that what you print make sense to me) - for example to print mean you can use cat command once you done with the calculation and have value in the variable
#variable_name=10 
#cat("Mean number of correctly received bits ",  variable_name)
n <- 24
p <- 0.9
bitMean <- (n*p)
cat ("Mean number of correctly received bits is", bitMean) 

Mean number of correctly received bits is 21.6

cat ("\n")

stDev <- sqrt((n*p)*(1-p))
cat ("Standard Deviation of correctly received bits is", stDev)

Standard Deviation of correctly received bits is 1.469694

Quesion 1 (c) What is the probability of more than \(3\)-bit errors in the block of \(24\)?

#Answer 1c
probMoreThanThree <- 1-pbinom(3, 24, 0.1)
cat("The probability of more than 3 bit errors in the block is", probMoreThanThree)

The probability of more than 3 bit errors in the block is 0.2142622

Quesion 1 (d) What is the median of this distribution? What can you interpret from the median value?

#Answer 1d
medianDist<-qbinom(0.5, 24, 0.9)
medianDist

[1] 22

cat ("The median is", medianDist, "and it can tell you the amount of times (22 or fewer) the correct bit will be received 50% of the time.")

The median is 22 and it can tell you the amount of times (22 or fewer) the correct bit will be received 50% of the time.

Quesion 1 (e) What is the 60th quantile of this distribution? What can you interpret from this value?

#Answer 1e
quant60 <- qbinom (0.6, 24, 0.9)
quant60

[1] 22

#Interpretation
cat("The 60th quantile is 22, meaning that 60% of the time, there will be 22 or less bits being received correctly.")

The 60th quantile is 22, meaning that 60% of the time, there will be 22 or less bits being received correctly.

Question2. The Public Service Answering Point (PSAP) in San Francisco employs \(19\) operators in \(8\)-hour shifts to process \(911\) calls. There are at least \(5\) operators always answering calls. The number of calls processed per operator can be modeled with a Poisson random variable with rate \(\lambda_0 =20\) calls per hour.

Quesion 2 (a) What is the probability an operator can process 20 calls in an hour? Repeat for 30 calls in an hour?

#Answer 2a
lambda <- 20
prob20 <- dpois (20, lambda)
prob30 <-dpois (30, lambda)
cat (prob20,  "is the probability an operator can process 20 calls in an hour")

0.08883532 is the probability an operator can process 20 calls in an hour

cat("\n")

cat (prob30, "is the probability an operator can process 30 calls in an hour")

0.008343536 is the probability an operator can process 30 calls in an hour

Quesion 2 (b) Given that 240 calls occurred in an hour, what is the probability that the ten operators can process them all assuming they are split equally among the operators?

#Answer 2b
##each op takes 24 calls wtihin an hour
prob24 <- dpois (24, lambda)
probAll <- (prob24)^10
probAll

[1] 2.892317e-13

Quesion 2 (c) Now, If the aggregate call rate per hour \(\lambda_c\) , is measured by \(\lambda_c = 85\) calls per hour, what is the probability of more than \(100\) calls in an hour

#Answer 2c
newLambda <- 85
Prob100<-1-ppois(100, 85)
cat("The probability of more than 100 calls in an hour is", Prob100)

The probability of more than 100 calls in an hour is 0.04934533

Question 3. For the two random number generator below A and B(don’t forget to add your R code)

- [A] \(Z_i = (9Z_{i-1} + 1) \mod 16\) with \(Z_0 = 5\).

- [B] \(Z_i = (7Z_{i-1} + 3) \mod 32\) with \(Z_0 = 10\),

Quesion 3 (a) Compute \(Z_i\) and \(U_i\) for values of \(i\) until a number is repeated, what is the period of both the generator? Provide your comments about the period of both RNG

#Answer 3a
seed <- 5
a <- 9
c<- 1
m<- 16
rlength <- 50

lcg <- function(a, c, m, rlength, seed) {
z <-rep (0, rlength)
z[1] <- seed

# loop to gen random numbers
for (i in 1:(rlength-1)) {
  #compute next lcg value
  z[i+1] <- (a*z[i] +c) %% m
}
U <- z / m
return (list(z=z,U=U))
}

lcg (9,1,16, 25, 5)

$z
 [1]  5 14 15  8  9  2  3 12 13  6  7  0  1 10 11  4  5 14 15  8  9  2  3 12 13

$U
 [1] 0.3125 0.8750 0.9375 0.5000 0.5625 0.1250 0.1875 0.7500 0.8125 0.3750 0.4375 0.0000 0.0625 0.6250 0.6875 0.2500 0.3125
[18] 0.8750 0.9375 0.5000 0.5625 0.1250 0.1875 0.7500 0.8125

cat ("The period for equation A is 16. This is equal to the modulus, meaning this is its full period, the same value as the mod.")

The period for equation A is 16. This is equal to the modulus, meaning this is its full period, the same value as the mod.

seed <- 10
a <- 7
c <- 3
m <- 32
rlength <- 50

lcg <- function(a, c, m, rlength, seed) {
z <- rep(0,rlength)
z[1] <- Z0_forB
for (i in 1:(rlength-1)) {
  z[i+1] <- (a*z[i] +c) %% m
}

#Calculate uniform rand value U in [0,1]
U <- z/m

#Return
return (list(z=z,U=U))
}
lcg(7,3,32,25,10)

$z
 [1] 10  9  2 17 26 25 18  1 10  9  2 17 26 25 18  1 10  9  2 17 26 25 18  1 10

$U
 [1] 0.31250 0.28125 0.06250 0.53125 0.81250 0.78125 0.56250 0.03125 0.31250 0.28125 0.06250 0.53125 0.81250 0.78125
[15] 0.56250 0.03125 0.31250 0.28125 0.06250 0.53125 0.81250 0.78125 0.56250 0.03125 0.31250

Question 3 (b) Which of these parameters effect the period of LCG – \(a\), \(b\), \(Z_0\)

#Answer 
cat("Parameters a and c affect the period of the LCG. The seed is multiplied by a and added to by c, and this number is operated upon by the modulo value.")

Parameters a and c affect the period of the LCG. The seed is multiplied by a and added to by c, and this number is operated upon by the modulo value.

Quesion 3 (c) For both generators plot a scatter diagram of the Zi values 1 apart. What are your observations from these plots? What do you think about the property of randomness of these generators?

#Answer
lcg <- function(a, c, m, rlength, seed) {
z <- rep(0,rlength)
z[1] <- Z0_forB
for (i in 1:(rlength-1)) {
  z[i+1] <- (a*z[i] +c) %% m
}

#Calculate uniform rand value U in [0,1]
U <- z/m

#Return
return (list(z=z,U=U)) }
example <- lcg(9,1,16, 25, 5)
ZT = example$z
Z1 =ZT[1:16]
Z2 = ZT[2:17]
plot (Z1, Z2, main = "Scatterplot of Zi Values for Equation A")

example <- lcg(7,3,32,25,10)
ZT = example$z
Z1 =ZT[1:8]
Z2 = ZT[2:9]
plot (Z1, Z2, main = "Scatterplot of Zi Values for Equation B")

cat("In scatterplot A, a visible and predictable pattern is shown. This means that LCG equation A is not very random. In scatterplot B, there is a less predictable pattern and the plot points are more scattered about. This shows that LCG B is more random, and can generate more unpreditable sequences than A.")

In scatterplot A, a visible and predictable pattern is shown. This means that LCG equation A is not very random. In scatterplot B, there is a less predictable pattern and the plot points are more scattered about. This shows that LCG B is more random, and can generate more unpreditable sequences than A.

Question 3 (d) Randomness of default random generator of R – Run runif command to generate 100 random numbers and plot the scatter diagram of these numbers ( values 1 apart) and discuss your observations about randomness of this random generator.

#Answer
rand <-runif(100)
rand

  [1] 0.57608160 0.90633661 0.86341068 0.06518454 0.94756522 0.31624360 0.39735815 0.29912225 0.61354231 0.62707194
 [11] 0.58216907 0.82012066 0.19860669 0.08905146 0.21158870 0.30366478 0.31175807 0.96213756 0.22735787 0.81010627
 [21] 0.13580325 0.55613351 0.19740285 0.78554348 0.88587381 0.28749707 0.84680518 0.56461377 0.35905576 0.26812346
 [31] 0.35646099 0.36033073 0.16581994 0.32082877 0.66402564 0.95992280 0.95425132 0.47297589 0.48864110 0.52958244
 [41] 0.99667129 0.69244686 0.95172444 0.63557902 0.66240638 0.12003894 0.71171426 0.52332448 0.49005332 0.14609978
 [51] 0.77039162 0.94400060 0.30546333 0.74528827 0.83473715 0.84479895 0.79791058 0.79453675 0.82380320 0.95854745
 [61] 0.03976804 0.65061744 0.12354842 0.01488115 0.42201990 0.31200090 0.98822245 0.60617974 0.38253127 0.97231882
 [71] 0.87166879 0.61465042 0.58090161 0.29415328 0.33103456 0.93238347 0.38675807 0.90356363 0.46713331 0.92869827
 [81] 0.23277721 0.25709528 0.86156281 0.69523689 0.14707560 0.35884933 0.07953704 0.20422794 0.15875490 0.77472802
 [91] 0.92148488 0.08272197 0.04469671 0.67531416 0.24090663 0.57905143 0.59784965 0.84457965 0.89290587 0.90498868

##Scatterplot

plot(rand[-length(rand)],rand[-1])

cat("This random generator is very unpredictable, as shown by the pattern on the scatterplot. There is no discernable pattern, and it appears that it is much stronger at generating unpredictable numbers than an LCG.")

This random generator is very unpredictable, as shown by the pattern on the scatterplot. There is no discernable pattern, and it appears that it is much stronger at generating unpredictable numbers than an LCG.

Question 3 (e) Compute the mean value of \(U_i\) across the period

#Answer
meanUi <- mean(rand)
cat("The mean value of Ui across the period is", meanUi, "\n")

The mean value of Ui across the period is 0.5441508 

Question 3 (f) By providing a plot of density (histogram) discuss the uniformity of both of the generators

#Answer
exampleA <- lcg(9,1,16,16, 5)
hdataA<- exampleA$U
hist(hdataA)

exampleB <- lcg(7,3,32,32,10)
hdataB<- exampleB$U
hist(hdataB)

cat("A shows more uniformity, supporting that it fully utilizes the whole period with its generated values. B on the other hand is more dispersed, showing it does not use its full period.")

A shows more uniformity, supporting that it fully utilizes the whole period with its generated values. B on the other hand is more dispersed, showing it does not use its full period.

Question 4. Using the inverse transform method

Question 4.(a) Develop an algorithm for the random variable with cumulative distribution function F(x) below.

\[ F(x) = 1 - e^{-(x/\lambda)^k} \]

where \(x \geq 0\), \(\lambda \geq 0\), and \(k \geq 0\).

# Answer 4a
# This might be hard to print using R - you can write the step in your notebook and then include image here / or you can include as seperate file when upload

# Check if the IRdisplay package is already installed
if (!require(IRdisplay, quietly = TRUE)) {
  # If not installed, install it
  install.packages("IRdisplay")
  
  # Load the IRdisplay library
  library(IRdisplay)
} else {
  # If already installed, just load the library
  library(IRdisplay)
}

inversetransform <- function(amount, lambda, k) {
  U <- runif(amount, 0,1)
  values <- rep(0, amount)
  i <- 1
for (i in 1: length(U)) {
  values[i] <- lambda * (nthroot(log(1-U[i]),k))
}  
  return (values=values)
  values
}


# Embed an image
#display_png(file = "example.png")

Question 4 (b) Use the first three Ui values from part (a) and generator [A] of previous problem to create 3 values from the random variable when, \(\lambda = 1\), \(𝑘 = 5\)

#Answer
Ui1 <- -1.0307
Ui2 <- -0.6685
Ui3 <- -0.5780

Question 4 (c) sing R generate 10,000 values from your algorithm when = 1 , 𝑘 = 5 and plot a histogram of the density. Discuss what insights you obtain when you look at the plot generated here vs the plot you generated in question 3 (e)

#Answer

Question 5. Develop a Monte Carlo simulation in R that counts the number of uniform [0,1] random numbers that must be summed to get a sum greater than 1. Run a single simulation with n = 10,000 times and find the mean of the number of counts. Do you think this number looks somewhat familiar. Every student might get slightly different value any ideas why?

#Answer
simulation_counts = function() {
  values <- 0
  count <- 0
  while (values <=1) {
    count = count +1
    values = values +runif(1)
  }
  return (count)
}


n <- 10000
rep = replicate(n, simulation_counts())
meanCounts = mean(rep)
cat ("The mean number of counts is", meanCounts, "\n")

The mean number of counts is 2.7131 

cat ("Students get different values because the numbers are randomly generated.")

Students get different values because the numbers are randomly generated.

Question 6. For our specific random variable, we know that its PDF is defined by the function \(f(x) = 10x(1-x)\). Utilize the accept-reject algorithm to draw samples from this distribution. Take a look at Lecture 10 slide 8 and 9

#Answer 

accept_reject <- function(f, c, g, rg, n) {
n.accepts <- 0
result.sample <- rep(NA, n)
while (n.accepts < n) {
x.star <- rg(1) # step 1
u <- runif(1,0,1) # step 2
if (u < f(x.star)/(c*g(x.star))) { # step 3 (accept)
n.accepts <- n.accepts+1
result.sample[n.accepts] = x.star
}
}
result.sample
}

f <- function(x) 10*x*(1-x) 
g <- function(x) x/x 
rg <- function(n) runif(n,0,1) 
c <- 2 
vals <- accept_reject(f, c, g, rg, 1000)
hist(vals,freq=FALSE, main="Histogram")

Quesion 1 (b)

#Answer