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
x <- 0:24
plot(x, dbinom(x, size=24, prob=0.9), xlab="Number of Bits Received", ylab="P(X = x)", type="h", main="Bits, n = 24, p = 0.9,", font.main=1)

x <- 0:24
plot(x, pbinom(x, size=24, prob=0.9), xlab="Number of Bits Received", ylab="P(X <=x)", ylim=c(0,1), type="s", main="Bits, n= 24, p=0.9", font.main=1)

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)

prob <- 0.9
x <- 0:24
mean_value <- round(dbinom(x, size=24, prob), 3)

simple_mean = 24 * 0.9
cat("Mean number of correctly received bits ", simple_mean, "\n")

Mean number of correctly received bits  21.6 

variance <- 1 - prob
standard_dev = sqrt(24 * variance * prob)
cat("Standard Deviation of correctly received bits ", standard_dev)

Standard Deviation of correctly received bits  1.469694

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

#Answer 1c

size = 24
prob = 0.9
final_prob = 1 - pbinom(3, size, prob)

cat("The probability of more than 3 bit errors occuing is ", final_prob, "\n")

The probability of more than 3 bit errors occuing is  1 

print("Note that the probability is so low that it is auto rounded to 1, true value is 1-1.498*10^(-18)")

[1] "Note that the probability is so low that it is auto rounded to 1, true value is 1-1.498*10^(-18)"

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

#Answer 1d

size = 24
prob = 0.9
median <- qbinom(0.5, size, prob)

cat("Interpretation: The median is", median, "this value tells us that half of the expected values will be less than that and half will be more. It also means this is middle or central value of the distribution")

Interpretation: The median is 22 this value tells us that half of the expected values will be less than that and half will be more. It also means this is middle or central value of the distribution

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

#Answer 1e

size = 24
prob = 0.9
quantile <- qbinom(0.6, size, prob)

#Interpretation
cat("Interpretation: The 60th quantile is", quantile,"this value represents the minimum number of successes required for the cumulative probability of success to be greater than or equal to 60%. In other words, it is the point in the distribution where there is a 60% probability or more of achieving at least that many successes.")

Interpretation: The 60th quantile is 22 this value represents the minimum number of successes required for the cumulative probability of success to be greater than or equal to 60%. In other words, it is the point in the distribution where there is a 60% probability or more of achieving at least that many successes.

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
given_mean = 20
x1 = 20
x2 = 30

call_a <- ((given_mean^x1)* (exp(-1 * given_mean))) / (factorial(x1))
call_b <- ((given_mean^x2)* (exp(-1 * given_mean))) / (factorial(x2))

cat("The probability that the operator can process 20 calls an hour is ", call_a,"\n")

The probability that the operator can process 20 calls an hour is  0.08883532 

cat("The probability that the operator can process 30 calls an hour is ", call_b)

The probability that the operator can process 30 calls an hour is  0.008343536

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
given_mean = 20
operator = 10
x1 = 24

#The value is risen to the number of operators to show the percentage being multiplied by the number of instances it has to occur in.
call_a <- ((given_mean^x1)* (exp(-1 * given_mean))) / (factorial(x1))
prob <- call_a^operator

cat("The probability that 10 operators can process 240 calls an hour is", prob)

The probability that 10 operators can process 240 calls an hour is 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

given_mean = 85
x1 = 100

call_a <-((given_mean^x1)* (exp(-1 * given_mean))) / (factorial(x1))
cat("The probability that the operator can process 20 calls an hour is ", call_a)

The probability that the operator can process 20 calls an hour is  0.01139877

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 RGN

#Answer 3a

lcg <- function(a,c,m,rlength,seed) {
    x <- rep(0,rlength)
    x[1] <- seed
    for (i in 1:(rlength-1))
    {
       x[i+1] <- (a * x[i] + c) %% m
    }
    U <- x/m 
    return(list(x=x,U=U))
}

z1 <- lcg(9,1,16,20,5)
z2 <- lcg(7,3,32,20,10)
z1

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

$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 0.8750 0.9375
[20] 0.5000

z2

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

$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 0.56250 0.03125 0.31250
[18] 0.28125 0.06250 0.53125

print("This shows that for LCG A the system has a period of 17 and LCG B has a period of 8. This means that form a security and randomness measurment LCG A is stronger as it takes longer before it repeats itself.")

[1] "This shows that for LCG A the system has a period of 17 and LCG B has a period of 8. This means that form a security and randomness measurment LCG A is stronger as it takes longer before it repeats itself."

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

#Answer 
cat("Answer:", "The parameters of 'a' and 'seed or Z0' effect the period of a LCG")

Answer: The parameters of 'a' and 'seed or Z0' effect the period of a LCG

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) {
    x <- rep(0,rlength)
    x[1] <- seed
    for (i in 1:(rlength-1))
    {
       x[i+1] <- (a * x[i] + c) %% m
    }
    U <- x/m 
    return(list(x=x,U=U))
}

z1 <- lcg(9,1,16,20,5)
z2 <- lcg(7,3,32,20,10)

ZT = z1$z
length(ZT)

[1] 0

z1 = ZT[1:19]
z2 = ZT[2:20]

plot(z1,z2, main="Scatter Plot of Zi Values", xlim = c(0, 30), ylim = c(0, 30))

Quesion 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

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

#Answer

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

#Answer

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)
}


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

Quesion 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

Quesion 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

n_sims <- 10000
counts <- numeric(n_sims)

# Perform the Monte Carlo simulation
for (i in 1:n_sims) {
  sum_numbers <- 0
  count <- 0
  while (sum_numbers <= 1) {
    sum_numbers <- sum_numbers + runif(1)
    count <- count + 1
  }
  counts[i] <- count
}

mean_count <- mean(counts)
cat("Mean count:", mean_count, "\n")

Mean count: 2.7091 

print("I'm unclear on why this number would be familiar. However, as for why every studetn would get a diffrent answer that is becuase the answer changes with every run of the simulation based on the numbers generated, creating diffrent results every time.")

[1] "I'm unclear on why this number would be familiar. However, as for why every studetn would get a diffrent answer that is becuase the answer changes with every run of the simulation based on the numbers generated, creating diffrent results every time."

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

Quesion 1 (b)

#Answer