Model Answers: Computational Testing
10.8 Problem 1
The data is:
x <- c(263.9, 266.2, 266.3, 266.8, 265.0)Now construct some summary statistics and define some given parameters:
x_bar <- mean(x) # compute sample mean
sigma <- 1.65 # population standard deviation is given
mu <- 260 # population mean to be tested against
n <- length(x) # number of samplesConstruct the z-statistic:
z <- (x_bar - mu) / (sigma / sqrt(n))
print(z)## [1] 7.643287Check if the z-statistic is in the critical range. First, work out what the z-value at the edge of the critical region is:
z_threshold <- qnorm(1 - 0.01, mean = 0, sd = 1)
print(z_threshold)## [1] 2.326348Thus, the z-statistic is much greater than the threshold and there is evidence to suggest the cartons are overfilled.
10.9 Problem 2
Parameters given by the problem:
x_bar <- 103.11
s <- 53.5
mu <- 100
n <- 45Compute the z-statistic assuming large sample assumptions apply:
z <- ( x_bar - mu )/(s/sqrt(n))
print(z)## [1] 0.3899535Now, work out the thresholds of the critical regions:
z_upper <- qnorm(1 - 0.025, mean = 0, sd = 1)
print(z_upper)## [1] 1.959964z_lower <- qnorm(0.025, mean = 0, sd = 1)
print(z_lower)## [1] -1.959964The z-statistic is outside the critical regions and therefore we do not reject the null hypothesis.
10.10 Problem 3
z_test <- function(x, mu, popvar){
one_tail_p <- NULL
z_score <- round((mean(x) - mu) / (popvar / sqrt(length(x))), 3)
one_tail_p <- round(pnorm(abs(z_score),lower.tail = FALSE), 3)
cat(" z =", z_score, "\n",
"one-tailed probability =", one_tail_p, "\n",
"two-tailed probability =", 2 * one_tail_p)
return(list(z = z_score, one_p = one_tail_p, two_p = 2 * one_tail_p))
}
x <- rnorm(10, mean = 0, sd = 1) # generate some artificial data from a N(0, 1)
out <- z_test(x, 0, 1) # null should not be rejected!## z = 0.387
## one-tailed probability = 0.349
## two-tailed probability = 0.698print(out)## $z
## [1] 0.387
##
## $one_p
## [1] 0.349
##
## $two_p
## [1] 0.698x <- rnorm(10, mean = 1, sd = 1) # generate some artificial data from a N(1, 1)
out <- z_test(x, 0, 1) # null should be rejected!## z = 2.836
## one-tailed probability = 0.002
## two-tailed probability = 0.004print(out)## $z
## [1] 2.836
##
## $one_p
## [1] 0.002
##
## $two_p
## [1] 0.00410.11 Problem 4
Define some parameters
mu <- 5.4
n <- 5
x_bar <- 5.64
s2 <- 0.05Compute the t-statistic:
t <- (x_bar - mu) / sqrt(s2 / n)
print(t)## [1] 2.4Work out the thresholds of the critical regions:
t_upper <- qt(1 - 0.025, df = n - 1)
print(t_upper)## [1] 2.776445t_lower <- qt(0.025, df = n - 1)
print(t_lower)## [1] -2.776445The t-statistic is outside of the critical regions so we do not reject the null hypothesis.
10.12 Problem 5
Define the parameters:
x_bar_a <- 88
s2_a <- 4.5
n_a <- 72
x_bar_b <- 79
s2_b <- 4.2
n_b <- 48
mu_a <- 0
mu_b <- 0Compute the z-statistic:
z <- ((x_bar_a - x_bar_b) - (mu_a - mu_b)) / sqrt(s2_a / n_a + s2_b / n_b)
print(z)## [1] 23.2379Work out for the 5% significance level, the critical values:
z_upper <- qnorm(1 - 0.05, mean = 0, sd = 1)
print(z_upper)## [1] 1.644854There is evidence to support the claim that process \(A\) yields higher pressurisation.
10.13 Problem 6
# Data vectors
x_A <- c(91.50, 94.18, 92.18, 95.39, 91.79, 89.07, 94.72, 89.21)
x_B <- c(89.19, 90.95, 90.46, 93.21, 97.19, 97.04, 91.07, 92.75)
# parameters based on data
x_bar_A <- mean(x_A)
s2_A <- var(x_A)
n_A <- length(x_A)
x_bar_B <- mean(x_B)
s2_B <- var(x_B)
n_B <- length(x_B)Compute the pooled variance estimator:
s2_p <- ((n_A - 1) * s2_A + (n_B - 1) * s2_B) / (n_A + n_B - 2)
print(s2_p)## [1] 7.294654Compute the t-statistic:
t = ( x_bar_A - x_bar_B ) / sqrt( s2_p*(1/n_A + 1/n_B) )
print(t)## [1] -0.3535909Work out the critical values:
t_upper <- qt(1 - 0.025, df = n_A + n_B - 2)
print(t_upper)## [1] 2.144787t_lower <- qt(0.025, df = n_A + n_B - 2)
print(t_lower)## [1] -2.144787Since \(|t|<2.14\) we have no evidence to reject the null hypothesis that the mean yields are equal.
Now, let us use the built-in t.test command:
out <- t.test(x = x_A, y = x_B, paired = FALSE, var.equal = TRUE,
conf.level = 0.95, mu = 0, alternative = "two.sided")
print(out)##
## Two Sample t-test
##
## data: x_A and x_B
## t = -0.35359, df = 14, p-value = 0.7289
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.373886 2.418886
## sample estimates:
## mean of x mean of y
## 92.2550 92.7325The options paired=FALSE means this is an unpaired t-test, var.equal=TRUE forces the estimated variances to be the same (i.e. we are using a pooled variance estimator) and we are testing at 95% confidence level with an alternative hypothesis that the true difference in means is non-zero.
The p-value for the t-test is between 0 and 1. In this case, the value is around 0.72 which means the hypothesis should not be reject.
10.14 Problem 7
Define parameters:
x <- c(23, 36, 31)
p <- c(1 / 3, 1 / 3, 1 / 3)
n <- sum(x)
K <- length(x)Compute the expected counts:
Ex = p*nCompute the chi-squared statistic:
chi2 <- sum((x - Ex)^2 / Ex)
print(chi2)## [1] 2.866667Compute the critical value form the chi-squared distribution:
chi_upper <- qchisq(1 - 0.05, df = K-1)
print(chi_upper)## [1] 5.991465Thus there is no evidence to reject the null hypothesis. The data provides no reason to suggest a preference for a particular door.
Now, we could have done this in R:
out <- chisq.test(x, p = c(1 / 3, 1 / 3, 1 / 3))
print(out)##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 2.8667, df = 2, p-value = 0.238510.15 Problem 8
y <- c( 0, 1, 2 )
x <- c( 32, 12, 6 )You will need the vcdExtra package to use the expand.dft command:
install.packages("vcdExtra")The expand.dft command allows one to convert the frequency table into a vector of samples:
library(vcdExtra)
samples <- expand.dft(data.frame(y,Frequency = x), freq = "Frequency")## Warning in type.convert.default(as.character(DF[[i]]), ...): 'as.is' should be
## specified by the caller; using TRUENow we can use the fitdistr function in the MASS package to estimate the MLE of the Poisson distribution
# loading the MASS package
library(MASS)
# fitting a Poisson distribution using maximum-likelihood
lambda_hat <- fitdistr(samples$y, densfun = 'Poisson')Let just solve this directly using R built in function. First compute the expected probabilities under the Poisson distribution using dpois to compute the Poisson pdf:
pr <- c(0, 0, 0)
pr[1] <- dpois(0, lambda = lambda_hat$estimate)
pr[2] <- dpois(1, lambda = lambda_hat$estimate)
pr[3] <- 1 - sum(pr[1:2])Then apply chisq.test:
out <- chisq.test(x, p = pr)## Warning in chisq.test(x, p = pr): Chi-squared approximation may be incorrectprint(out)##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 1.3447, df = 2, p-value = 0.5105Actually, in this case the answer is wrong(!), we need to apply an additional loss of degree of freedom to account for the use of the MLE. However, we can re-use values already computed by chisq.test:
chi2 <- out$statistic
print(chi2)## X-squared
## 1.34466chi2_lower <- qchisq(1 - 0.01, df = 1)
print(chi2_lower)## [1] 6.634897Hence, there is no evidence to reject the null hypothesis. There is no reason to suppose that the Poisson distribution is not a plausible model for the number of accidents per week at this junction.