Directions:

Submit the project as a group or as an individual work.

Note: If you Rmd file submission knits you will receive total of (5 points)

1. PROBLEM SOLVING PART

Problem 1 (5 pts) - similar to Pr. 5 / page 238 in text A box contains 30 consecutive balls numbered 1 to 30. If four numbers are drawn at random, how many ways are there for the largest number to be 19 and the smallest number to be 9?

Use the Fundamental Principle of Counting!

Solution:

YOUR CODE HERE:

# Define the range of middle numbers
middle_numbers <- 10:18

# Calculate the number of ways to choose 2 numbers from the middle numbers
number_of_ways <- choose(length(middle_numbers), 2)

# Print the result
print(number_of_ways)
## [1] 36

Problem 2 (5 pts) - similar to Pr. 14 / page 239 in text On a multiple-choice exam with four possible answers for each of the five questions, what is the probability that a student would get four or more correct answers just by guessing?

Hint: Use the fact that \(P(E ∩ F) = P(E) · P(F)\) for two independent event (generalized for more than events) Getting one answer correct is independent of another.

Also

\[P(\text{at least } 4) = P(\text{exactly } 4) + P(\text{exactly } 5)\]

the last events are mutually exclusive so \(P(A \cup B) = P(A) + P(B)\)

YOUR CODE HERE:

# Define the probability of guessing a question correctly
p_correct <- 1/4

# Calculate the probability of getting exactly 4 correct answers
p_exactly_4 <- dbinom(4, size = 5, prob = p_correct)

# Calculate the probability of getting all 5 correct answers
p_exactly_5 <- dbinom(5, size = 5, prob = p_correct)

# Calculate the probability of getting at least 4 correct answers
p_at_least_4 <- p_exactly_4 + p_exactly_5

# Print the result
print(p_at_least_4)
## [1] 0.015625

Out of 200 students who adopt this test taking approach how many are expected to get at least 4 correct? Hint: Use Binomial experiment settings to answer this questions.

Solution:

3.125

Problem 3 (5 pts) (IND-OPT) - similar to Pr. 44 page 243 in text Consider tossing four fair coins. There are 16 possible outcomes, e.g HHHH,HHHT,HHTH,HTHH,THHH, ... possible outcomes. Define X as the random variable “number of heads showing when four coins are tossed.” Obtain the mean and the variance of X. Simulate tossing four fair coins 10,000 times. Compute the simulated mean and variance of X. Are the simulated values within 2% of the theoretical answers?

Hint: To find the theoretical values use dbinom (x= , size = , prob = )

Solution:

YOUR CODE HERE:

library(MASS)
# Simulate tossing four fair coins 10,000 times
set.seed(123) # Setting a seed for reproducibility
tosses <- replicate(10000, sum(sample(c(0, 1), size = 4, replace = TRUE)))

# Compute the simulated mean and variance
sim_mean <- mean(tosses)
sim_variance <- var(tosses)

# Theoretical mean and variance
theo_mean <- 2
theo_variance <- 1

# Check if the simulated values are within 2% of the theoretical values
within_2_percent_mean <- abs(sim_mean - theo_mean) / theo_mean <= 0.02
within_2_percent_variance <- abs(sim_variance - theo_variance) / theo_variance <= 0.02

# Results
list(simulated_mean = sim_mean, 
     simulated_variance = sim_variance, 
     within_2_percent_mean = within_2_percent_mean, 
     within_2_percent_variance = within_2_percent_variance)
## $simulated_mean
## [1] 1.9938
## 
## $simulated_variance
## [1] 1.013063
## 
## $within_2_percent_mean
## [1] TRUE
## 
## $within_2_percent_variance
## [1] TRUE

Problem 4 (10 pts) - similar to Pr. 11 / page 307 in text

Traffic volume is an important factor for determining the most cost-effective method to surface a road. Suppose that the average number of vehicles passing a certain point on a road is 5 every 60 seconds.

  1. Find the probability that more than 12 cars will pass the point in a period of 3 minutes.
  2. What is the probability that more than 100 cars pass the point in half an hour?

Hint: Adjust the Poisson parameter \(\lambda = 5\) every 60 seconds to # per every 3 minutes

YOUR CODE HERE:

1 - ppois(11, 15)
## [1] 0.8152482
1 - ppois(99, 150)
## [1] 0.9999941

Problem 5 (10 pts) - similar to Pr. 18 / page 308 in text Suppose the percentage of drinks sold from a vending machine are 70% and 30% for soft drinks and bottled water, respectively.

  1. What is the probability that on a randomly selected day, the first soft drink is the fourth drink sold?
  2. Find the probability that exactly 4 out of 10 drinks sold is a soft drink.

Hint: Let X = number of waters (failures) purchased before the first soft drink is purchased. Then, \(X \sim Geo(0.70)\).

Solution:

YOUR CODE HERE::

# Part (a)
p_success <- 0.70
prob_first_soft_drink_fourth <- (1 - p_success)^3 * p_success
# Result
prob_first_soft_drink_fourth
## [1] 0.0189
  1. Let X = number of soft drinks sold. Then, \(X \sim Bin(10, 0.70)\) and \(P(X = 1) = 0\) since one has:
# Part (b)
n_drinks <- 10
k_soft_drinks <- 4
prob_exactly_4_soft_drinks <- choose(n_drinks, k_soft_drinks) * p_success^k_soft_drinks * (1 - p_success)^(n_drinks - k_soft_drinks)

# Result
prob_exactly_4_soft_drinks
## [1] 0.03675691

Problem 6 (10 pts) - Exponential Distribution: Light Bulbs If the life of a certain type of light bulb has an exponential distribution with a mean of 10 months, find

  1. The probability that a randomly selected light bulb lasts between 6 and 15 months.
  2. The 95th percentile of the distribution. How the probability of the light bulb to last more than 36 months compare to 0.05?
  3. The probability that a light bulb that has lasted for 10 months will last more than 25 months.

Solution:

YOUR CODE HERE:

# Part (a): Probability that a light bulb lasts between 6 and 15 months
lambda <- 1/10 # Rate parameter for the exponential distribution
prob_6_to_15 <- pexp(15, rate = lambda) - pexp(6, rate = lambda)

# Part (b): 95th percentile of the distribution
percentile_95 <- qexp(0.95, rate = lambda)
# Probability of lasting more than 36 months
prob_over_36 <- pexp(36, rate = lambda, lower.tail = FALSE)
# Compare to 0.05
comparison <- prob_over_36 < 0.05

# Part (c): Probability that a bulb lasts more than 25 months given it has lasted 10 months
# Using the memoryless property of the exponential distribution
prob_over_25_given_10 <- pexp(15, rate = lambda, lower.tail = FALSE)

# Output the results
list(
  prob_6_to_15 = prob_6_to_15,
  percentile_95 = percentile_95,
  prob_over_36 = prob_over_36,
  comparison_to_05 = comparison,
  prob_over_25_given_10 = prob_over_25_given_10
)
## $prob_6_to_15
## [1] 0.3256815
## 
## $percentile_95
## [1] 29.95732
## 
## $prob_over_36
## [1] 0.02732372
## 
## $comparison_to_05
## [1] TRUE
## 
## $prob_over_25_given_10
## [1] 0.2231302

#Problem 7 (10 pts) (IND-OPT) - Pr. 8 page 400 A population has the following elements: 2, 5, 8, 12, 13 (finite population).

  1. Enumerate all the samples of size 2 that can be drawn with and without replacement. (Hint: Use the srs() function from the PASWR2 package)

  2. Calculate the mean of the population.

  3. Calculate the variance of the population.

  4. Calculate the standard deviation of the population.

  5. Calculate the mean of the sample mean, \(E[X]\).

  6. Calculate the variance of the sampled mean, \(Var(\bar X)\)

  7. Calculate the standard deviation of the sample mean.

  8. Calculate the mean of the sample variance, \(E[S^2]\)

  9. Is the variance of \(\bar X\) larger when sampling with or without replacement? Explain your answer

Solution:

YOUR CODE HERE:

Problem 8 (10 pts) - similar to Pr. 10 / page 400 in text

Use the data frame WHEATUSA2004 from the PASWR2 package; draw all samples of sizes 4, 5, and 6; and calculate the mean of the means for the variable acres (wheat surface area measured in thousands of acres). . What size provides the best approximation to the population mean? What is the variance of these means?

Hint: Use the srs() function from the PASWR2 package to draw simple random samples.

YOUR CODE HERE:

SRS2 <- srs(WHEATUSA2004$acres, 4)
xbarSRS2 <- apply(SRS2, 1, mean)
c(mean(xbarSRS2), var(xbarSRS2))
## [1]   1148.733 645755.872
SRS3 <- srs(WHEATUSA2004$acres, 5)
xbarSRS3 <- apply(SRS3, 1, mean)
c(mean(xbarSRS3), var(xbarSRS3))
## [1]   1148.733 496720.646
SRS4 <- srs(WHEATUSA2004$acres, 6)
xbarSRS4 <- apply(SRS4, 1, mean)
c(mean(xbarSRS4), var(xbarSRS4))
## [1]   1148.733 397374.397

Problem 9 (10 pts) - Pr. 16 / page 513 in text

A group of engineers working with physicians in a research hospital is developing a new device to measure blood glucose levels. Based on measurements taken from patients in a previous study, the physicians assert that the new device provides blood glucose levels slightly higher than those provided by the old device. To corroborate their suspicion, 15 diabetic patients were randomly selected, and their blood glucose levels were measured with both the new and the old devices. The measurements, in mg/100 ml, appear in data frame GLUCOSE from the PASWR2 package:

  1. Are the samples independent? Why or why not?
  2. If the blood glucose level is a normally distributed random variable, compute a 95% confidence interval for the mean differences of the population.
  3. Use the results in (b) to decide whether or not the two devices give the same results.

Solution:

  1. The samples aren’t independent because the measurments are dependent on both the old and new devices and are related to each other.

  2. Hint: look at the QQ-plot to see if normality is reasonable.

YOUR CODE HERE:

GLUCOSE <- GLUCOSE %>% mutate(DIFF = old-new) # create variable DIFF for the difference between old and new level

head(GLUCOSE)
ggplot(data = GLUCOSE, aes(sample = DIFF)) +
  stat_qq() +
  theme_bw()

CI <- t.test(GLUCOSE$DIFF)$conf # use t-test since the sample size is small < 30
CI
## [1] -16.36145 -12.75855
## attr(,"conf.level")
## [1] 0.95

The 95% confidence interval for the mean differences of the population is [−16.3614,−12.7586].

  1. The confidence interval suggests that on average, the new device reports higher blood glucose levels than the old device.

#Problem 10 (10 pts) (IND-OPT) - similar to Pr. 9 / page 511 in text

A large company wants to estimate the proportion of its accounts that are paid on time.

  1. How large a sample is needed to estimate the true proportion within 2% with a 95% confidence level?

  2. Suppose 700 out of 800 accounts are paid on time. Construct 95% confidence intervals for the true proportion of accounts that are paid on time using an asymptotic confidence interval, an Agresti Coull confidence interval.

Hint: Use the nsize() and binom.confint() from the library(binom)

Solution:

YOUR CODE HERE:

library(binom)
## Warning: package 'binom' was built under R version 4.3.3
library(binom)
#(a) 

#(b)

Problem 11 (10 pts) - Pr. 10 / page 511 in text

In a study conducted at Appalachian State University, students used digital oral thermometers to record their temperatures each day they came to class. A randomly selected day of student temperatures is provided in the following table and in the data frame STATTEMPS. Information is also provided with regard to subject gender and the hour of the day when the students’ temperatures were measured.

Direction: load the data with data(STATTEMPS) having the library(PASWR2) loaded.

Hint: Use t.test(temperature ~ gender, data = STATTEMPS, mu = 0, paired = FALSE) the x ad y values are given with formula expression formula = temperature ~ gender. Then use the $ access to obtain the confidence interval: CI <- t.test(temperature ~ gender, data = STATTEMPS, mu = 0, paired = FALSE)$conf

  1. Construct a 95% confidence interval for the true average temperature difference between females and males. Does the interval contain the value zero? What does this suggest about gender temperature differences?

  2. Construct a 95% confidence interval for the true average temperature difference between students taking their temperatures at 8 a.m. and students taking their temperatures at 9 a.m. Give a reason why one group appears to have a higher temperature reading.

Solution:

YOUR CODE HERE:

library(PASWR2)
data(STATTEMPS)
#(a)
t.test(temperature ~ gender, data = STATTEMPS, mu = 0, paired = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  temperature by gender
## t = 0.467, df = 19.026, p-value = 0.6458
## alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
## 95 percent confidence interval:
##  -1.054064  1.659598
## sample estimates:
## mean in group Female   mean in group Male 
##             95.43913             95.13636
CI <- t.test(temperature ~ gender, data = STATTEMPS, mu = 0, paired = FALSE)$conf


#(b)
t.test(temperature ~ class, data = STATTEMPS, mu = 0, paired = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  temperature by class
## t = -2.514, df = 28.858, p-value = 0.01777
## alternative hypothesis: true difference in means between group 8 a.m. and group 9 a.m. is not equal to 0
## 95 percent confidence interval:
##  -2.496501 -0.256440
## sample estimates:
## mean in group 8 a.m. mean in group 9 a.m. 
##             94.65294             96.02941
CI <- t.test(temperature ~ class, data = STATTEMPS, mu = 0, paired = FALSE)$conf

A 95% confidence interval for the true average difference between students taking their temperatures at 8 a.m. and students taking their temperatures at 9 a.m. is [−2.4965,−0.2564]. Note that this interval does not contain 0, indicating that the there is evidence to suggest students in the 8 a.m. class have temperatures that are not as warm as the 9 a.m. class. One possible explanation is that students roll straight out of bed and into the 8 a.m. class. Consequently, their temperatures are closer to their sleeping temperatures which are lower than their waking temperatures.

2. BRAINSTORMING PART:

How do we form confidence intervals when normality assumption is not reasonable and we have small sample size so the Central Limit Theorem does not apply?

Look for the answer with the following :

Empirical bootstrap confidence interval for the mean.

For the data contained in a vector x = c(36,30,37,43,42,43,43,46,40,42) Estimate the mean μ of the underlying distribution and give an 90% bootstrap confidence interval.

x = c(36,30,37,43,42,43,43,46,40,42)
qqnorm(x) # see if the data appear to follow normal distribution

Data is not quite normal as it appears on the Q-Q PLOT!

n = length(x)
set.seed(25)  

# sample mean
xbar = mean(x)

cat("data mean = ",xbar,'\n')
## data mean =  40.2
nboot = 10000 # number of bootstrap samples 

# Generate 10000 bootstrap samples, i.e. an n x 10000 array of random resamples from x.
tmpdata = sample(x,n*nboot, replace=TRUE) # sample with replacement
bootstrapsample = matrix(tmpdata, nrow=n, ncol=nboot)

# Compute the sample mean xbar for each bootstrap sample
xbarstar = colMeans(bootstrapsample)

# Compute delta* for each bootstrap sample (difference between the original mean and the bootstrap sample mean)
deltastar = xbarstar - xbar

# Find the 0.1 and 0.9 quantiles for deltastar
d = quantile(deltastar,c(0.05,0.95)) # for the 905 CI we need 90% area between the quantiles


# Calculate the 90\% confidence interval for the mean.
ci = xbar - c(d[2],d[1])
ci
##  95%   5% 
## 38.1 42.6

We can achieve the same CI with the functions boot() and boot.ci() from the package boot

MEAN <- function(data, i){
  d <- data[i]
   M <- mean(d)
  }
boot.obj <- boot(x, statistic = MEAN, 10000)
CI <- boot.ci(boot.obj, conf = 0.90,type = "all") # all type of bootstrap intervals will be returned
## Warning in boot.ci(boot.obj, conf = 0.9, type = "all"): bootstrap variances
## needed for studentized intervals
CI
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.obj, conf = 0.9, type = "all")
## 
## Intervals : 
## Level      Normal              Basic         
## 90%   (37.91, 42.47 )   (38.00, 42.60 )  
## 
## Level     Percentile            BCa          
## 90%   (37.8, 42.4 )   (37.3, 42.0 )  
## Calculations and Intervals on Original Scale

3. BOOTSTRAP ESTIMATION PART

Review the problem below.

Pr. 30 / page 695 in text

The “Wisconsin Card Sorting Test” is widely used by psychiatrists, neurologists, and neuropsychologists with patients who have a brain injury, neurodegenerative disease, or a mental illness such as schizophrenia. Patients with any sort of frontal lobe lesion generally do poorly on the test. The data frame WCST and the following table contain the test scores from a group of 50 patients from the Virgen del Camino Hospital (Pamplona, Spain).

  1. Use the function eda() from the PASWR2 package to explore the data and decide if normality can be assumed. For details type in console ?eda

  2. What assumption(s) must be made to compute a 95% confidence interval for the population mean?

  3. Compute the confidence interval from (b).

  4. Compute a 95% BCa bootstrap confidence interval for the mean test score.

  5. Should you use the confidence interval reported in (c) or the confidence interval reported in (d)?

Solution:

  1. Assuming the variable score has a normal distribution is not reasonable.

See the results of EDA analysis below:

data(WCST)

WCST %>% pull(score) %>% eda()

## Size (n)  Missing  Minimum   1st Qu     Mean   Median   TrMean   3rd Qu 
##   50.000    0.000    4.000    8.500   21.480   17.000   19.413   26.000 
##      Max    Stdev      Var  SE Mean   I.Q.R.    Range Kurtosis Skewness 
##   94.000   18.406  338.785    2.603   17.500   90.000    4.511    2.033 
## SW p-val 
##    0.000

  1. In order to construct a 95% confidence interval for the population mean, one assumes that the values in the variable score are taken from a normal distribution. Although this is not a reasonable assumption, the sample size might be sufficiently large to overcome the skewness in the parent population. Consequently, one might appeal to the Central Limit Theorem and claim that the sampling distribution of X is approximately normal due to the sample size (50). In this problem, the skewness is quite severe, and one should not be overly confident in the final interval.

  2. If we assume normality one has:

# CI <- with(data = WCST, t.test(score)$conf) # 

# or using piping
CI <- WCST %>% pull(score) %>% t.test() %>% .$conf
CI
## [1] 16.24904 26.71096
## attr(,"conf.level")
## [1] 0.95
  1. Use the boot.ci() nonparametric bootstrap CIs functions to obtain the bootstrap CI.
library(boot) # use boot package
MEAN <- function(data, i){
  d <- data[i]
   M <- mean(d)
  }

set.seed(13)

B <- 10^4 - 1 # number of bootstrap samples

b.obj <- boot(data = WCST$score, statistic = MEAN, R = B) # bootstrap object of class "boot" containing the output of a bootstrap calculation
CIB <- boot.ci(b.obj, conf = 0.95, type = "bca")
CIB
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 9999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = b.obj, conf = 0.95, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   (17.26, 27.86 )  
## Calculations and Intervals on Original Scale

Problem 12 (10 pts) Set the seed to 23 (group 1), 45 (group 2), 67 (group 3), and 89 (group 4). Draw random sample of size 200 from exponential distribution with mean \(\lambda = 4\) (\(rate = \frac{1}{4}\)). Produce all bootstrap CI with the function boot.ci (use type = “all”).

How many of them contain the true mean for the sampled distribution?

YOUR CODE HERE:

set.seed(45) 
my.mean = function(x, indx){
  mean( x[indx] )
}
exps200 <- rexp(200,1/4)

b.obj <- boot(exps200, statistic = my.mean, 10000)
CIB <- boot.ci(b.obj, type = "all")
## Warning in boot.ci(b.obj, type = "all"): bootstrap variances needed for
## studentized intervals
CIB
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = b.obj, type = "all")
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 3.729,  4.935 )   ( 3.716,  4.919 )  
## 
## Level     Percentile            BCa          
## 95%   ( 3.737,  4.940 )   ( 3.769,  4.978 )  
## Calculations and Intervals on Original Scale

EXTRA CREDIT (10 pts) Set the seed to 23 (group 1), 45 (group 2), 67 (group 3), and 89 (group 4). Draw random sample of size 25 from \(Bin(20,\frac{1}{5})\). Produce all bootstrap CI with the function boot.ci (use type = “all”).

How many of them contain the true mean for the \(Bin(20,\frac{1}{5})\)?

YOUR CODE HERE:

set.seed(45) 
my.mean = function(x, indx){
  mean( x[indx] )
}
exps25 <- rexp(25,1/5)

b.obj <- boot(exps25, statistic = my.mean, 10000)
CIB <- boot.ci(b.obj, type = "all")
## Warning in boot.ci(b.obj, type = "all"): bootstrap variances needed for
## studentized intervals
CIB
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = b.obj, type = "all")
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 2.668,  6.688 )   ( 2.369,  6.369 )  
## 
## Level     Percentile            BCa          
## 95%   ( 2.969,  6.969 )   ( 3.290,  7.954 )  
## Calculations and Intervals on Original Scale