INFSCI 0310 HW3

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
p <- .9
x <- 0:n
dbinom(x, size = n, prob = p)

 [1] 1.000000e-24 2.160000e-22 2.235600e-20 1.475496e-18 6.971719e-17 2.509819e-15
 [7] 7.152983e-14 1.655405e-12 3.165961e-11 5.065538e-10 6.838477e-09 7.833164e-08
[13] 7.637335e-07 6.344863e-06 4.486725e-05 2.692035e-04 1.362843e-03 5.772039e-03
[19] 2.020214e-02 5.741660e-02 1.291874e-01 2.214640e-01 2.717968e-01 2.127105e-01
[25] 7.976644e-02

plot(x, dbinom(x, size = n, prob = p), main = "Probability Density Function(PDF)",
  xlab = "Correctly Received Bits",
  ylab = "Probability")

#Answer 1a
pbinom(x, size = n, prob = p)

 [1] 1.000000e-24 2.170000e-22 2.257300e-20 1.498069e-18 7.121525e-17 2.581034e-15
 [7] 7.411087e-14 1.729516e-12 3.338913e-11 5.399430e-10 7.378420e-09 8.571006e-08
[13] 8.494436e-07 7.194307e-06 5.206155e-05 3.212650e-04 1.684108e-03 7.456147e-03
[19] 2.765828e-02 8.507489e-02 2.142622e-01 4.357263e-01 7.075230e-01 9.202336e-01
[25] 1.000000e+00

plot(x, pbinom(x, size = n, prob = p), main = "Cumulative Distribution Function(CDF)",
     xlab = "Correctly Received Bits",
     ylab = "Probability")

NA

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 <- .9

mean <- n * p

cat("Mean number of correctly received bits: ", mean)

Mean number of correctly received bits:  21.6

n <- 24
p <- .9

stdev <- sqrt(n * p * (1-p))

cat("Standard Deviation of correctly received bits: ", stdev)

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
n <- 24
p <- 0.9

prob20 <- pbinom(20, size = n, prob = p)
proberror <- 1 - prob20

cat("Probability of more than 3-bit errors in the block of 24: ", proberror)

Probability of more than 3-bit errors in the block of 24:  0.7857378

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

#Answer 1d

n <- 24
p <- 0.9

median <- 0
for (i in 0:n) {
  if(pbinom(i, size = n, prob = p) >= 0.5) {
    median <- i
    break
  }
  
}

cat("The median value for this distribution: ", median,"->", " The median is the middle value is a set of numbers. In this case the median represents the middle value in the distribution of 24 bit blocks.")

The median value for this distribution:  22 ->  The median is the middle value is a set of numbers. In this case the median represents the middle value in the distribution of 24 bit blocks.

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

#Answer 1e

n <- 24
p <- 0.9

quantile60 <- qbinom(0.6, size = n, prob = p)

cat("The 60th Quartile of this distribution: ", quantile60, "->", "The 60th quartile represents that 40 pe" )

The 60th Quartile of this distribution:  22 -> The 60th quartile represents that 40 pe

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

lambda0 <- 20

prob20 <- dpois(20, lambda0)
prob30 <- dpois(30, lambda0)

cat("Probability an operator can process 20 calls an hour: ", prob20, "\n")

Probability an operator can process 20 calls an hour:  0.08883532 

cat("Probabilty an operator can process 30 calls an hour: ", prob30)

Probabilty an operator can process 30 calls an hour:  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

lambda0 <- 20
operators <- 10
callsperoperator <- 24

prob24peroperator <- dpois(callsperoperator, lambda0)
prob10operators <- prob24peroperator^operators

cat("The Probability that 10 operators can process 24 calls each: ", prob10operators)

The Probability that 10 operators can process 24 calls each:  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

lambdaC <- 85

prob100 <- ppois(100, lambdaC)
prob100more <- 1 - prob100

cat("The probability of more than 100 calls (when Lambda = 85) is: ", prob100more)

The probability of more than 100 calls (when Lambda = 85) 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 RGN

#Answer 3a
#Generator A

Z_0_A <- 5
Z_values_A <- c(Z_0_A)
U_values_A <- c(Z_0_A / 16)

while(TRUE) {
  Z_next_A <- (9 * Z_values_A[length(Z_values_A)] + 1) %% 16
  if(Z_next_A %in% Z_values_A) break
  Z_values_A <- append(Z_values_A, Z_next_A)
  U_values_A <- append(U_values_A, Z_next_A / 16)
}


cat("Period of A: ", length(Z_values_A))

Period of A:  16

#Answer 3a
#Generator B

Z_0_B <- 10
Z_values_B <- c(Z_0_B)
U_values_B <- c(Z_0_B / 32)

while(TRUE) {
  Z_next_B <- (7 * Z_values_B[length(Z_values_B)] + 3) %% 32
  if(Z_next_B %in% Z_values_B) break
  Z_values_B <- append(Z_values_B, Z_next_B)
  U_values_B <- append(U_values_B, Z_next_B / 32)
}


cat("Period of B: ", length(Z_values_B))

Period of B:  8

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

#Answer 
cat("Answer", "The parameter a, b, and z")

Answer The parameter a, b, and z

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
plot(Z_values_A[-length(Z_values_A)], Z_values_A[-1], main="Scatter Diagram for Generator A", xlab="Zi", ylab="Zi+1")

plot(Z_values_B[-length(Z_values_B)], Z_values_B[-1], main="Scatter Diagram for Generator B", xlab="Zi", ylab="Zi+1")

cat("The first plot shows a slight pattern while the second plot is mostly random. These plots show that the generators are producing random number. This is because if the generators were not producing random numbers, there would be stronger patterns in either of the scatterplots.")

The first plot shows a slight pattern while the second plot is mostly random. These plots show that the generators are producing random number. This is because if the generators were not producing random numbers, there would be stronger patterns in either of the scatterplots.

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

randomnumbers <- runif(100)

plot(randomnumbers[-100], randomnumbers[-1], xlab="Value at i", ylab="Value at i+1", main="Scatter Diagram of Random Numbers (1 Apart)")

cat("There are not many observations we can make about this scatterplot other than it very random. This is a great example of randomness as there is no observable pattern or strength (postive or negative). ")

There are not many observations we can make about this scatterplot other than it very random. This is a great example of randomness as there is no observable pattern or strength (postive or negative). 

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

#Answer

Z <- 5
Z_values_A <- c(Z)
m_A <- 16

while (TRUE) {
  Z <- (9 * Z + 1) %% m_A
  if (Z %in% Z_values_A) {
    break
  }
  Z_values_A <- c(Z_values_A, Z)
}

U_values_A <- Z_values_A / m_A
mean <- mean(U_values_A)

cat("Mean value across the period: ", mean)

Mean value across the period:  0.46875

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

#Answer

hist(U_values_A, breaks=20, main="Histogram of U_i Values for Generator A", xlab="U_i", probability=TRUE)

hist(U_values_B, breaks=20, main="Histogram of U_i Values for Generator B", xlab="U_i", probability=TRUE)

cat("Both generators provide uniform histograms. Generator B looks more uniform than Generator A as it is grouped better.")

Both generators provide uniform histograms. Generator B looks more uniform than Generator A as it is grouped better.

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

simulation <- function() {
  sum <- 0
  count <- 0
  while (sum <= 1) {
    sum <- sum + runif(1)
    count <- count + 1
  }
  return(count)
}

n <- 10000
results <- replicate(n, simulation())


mean_count <- mean(results)

cat("Mean number of counts: ", mean_count, "\n")

Mean number of counts:  2.7119 

cat("Answers will be different, because a Monte Carlo simulation has randomly generated numbers, therefore everytime we run this code chunk, we will get a different answer.")

Answers will be different, because a Monte Carlo simulation has randomly generated numbers, therefore everytime we run this code chunk, we will get a different answer.

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

f <- function(x) {
  10 * x * (1 - x)
}

g <- function(x) {
  1
}

c <- 2.5

acceptrejectsample <- function() {
  while(TRUE) {
    X <- runif(1)
    U <- runif(1)
    if(U <= f(X) / (c * g(X))) {
      return(X)
    }
  }
}

sample <- 10000
samples <- replicate(sample, acceptrejectsample())
hist(samples, main="Histogram of Samples", xlab="x", breaks=50)

Quesion 1 (b)

#Answer