Skip to main content

Essentials of Mathematics and Statistics

G Solution: Computational Testing

G.1 Solution: Exercise 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 samples

Construct the z-statistic:

z <- (x_bar - mu) / (sigma / sqrt(n))
print(z)
#> [1] 7.643287

Check 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.326348

Thus, the z-statistic is much greater than the threshold and there is evidence to suggest the cartons are overfilled.

G.2 Solution: Exercise 2

Parameters given by the problem:

x_bar <- 103.11
s <- 53.5
mu <- 100
n <- 45

Compute the z-statistic assuming large sample assumptions apply:

z <- ( x_bar - mu )/(s/sqrt(n))
print(z)
#> [1] 0.3899535

Now, work out the thresholds of the critical regions:

z_upper <- qnorm(1 - 0.025, mean = 0, sd = 1)
print(z_upper)
#> [1] 1.959964
z_lower <- qnorm(0.025, mean = 0, sd = 1)
print(z_lower)
#> [1] -1.959964

The z-statistic is outside the critical regions and therefore we do not reject the null hypothesis.

G.3 Solution: Exercise 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.055 
#>  one-tailed probability = 0.478 
#>  two-tailed probability = 0.956
print(out)
#> $z
#> [1] 0.055
#> 
#> $one_p
#> [1] 0.478
#> 
#> $two_p
#> [1] 0.956
x <- 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.584 
#>  one-tailed probability = 0.005 
#>  two-tailed probability = 0.01
print(out)
#> $z
#> [1] 2.584
#> 
#> $one_p
#> [1] 0.005
#> 
#> $two_p
#> [1] 0.01

G.4 Solution: Exercise 4

Define some parameters

mu <- 5.4
n <- 5
x_bar <- 5.64
s2 <- 0.05

Compute the t-statistic:

t <- (x_bar - mu) / sqrt(s2 / n)
print(t)
#> [1] 2.4

Work out the thresholds of the critical regions:

t_upper <- qt(1 - 0.025, df = n - 1)
print(t_upper)
#> [1] 2.776445
t_lower <- qt(0.025, df = n - 1)
print(t_lower)
#> [1] -2.776445

The t-statistic is outside of the critical regions so we do not reject the null hypothesis.

G.5 Solution: Exercise 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 <- 0

Compute 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.2379

Work out for the 5% significance level, the critical values:

z_upper <- qnorm(1 - 0.05, mean = 0, sd = 1)
print(z_upper)
#> [1] 1.644854

There is evidence to support the claim that process \(A\) yields higher pressurisation.

G.6 Solution: Exercise 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.294654

Compute the t-statistic:

t = ( x_bar_A - x_bar_B ) / sqrt( s2_p*(1/n_A + 1/n_B) )
print(t)
#> [1] -0.3535909

Work out the critical values:

t_upper <- qt(1 - 0.025, df = n_A + n_B - 2)
print(t_upper)
#> [1] 2.144787
t_lower <- qt(0.025, df = n_A + n_B - 2)
print(t_lower)
#> [1] -2.144787

Since \(|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.7325

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

G.7 Solution: Exercise 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*n

Compute the chi-squared statistic:

chi2 <- sum((x - Ex)^2 / Ex)
print(chi2)
#> [1] 2.866667

Compute the critical value form the chi-squared distribution:

chi_upper <- qchisq(1 - 0.05, df = K-1)
print(chi_upper)
#> [1] 5.991465

Thus 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.2385

G.8 Solution: Exercise 8 

y <- c( 0, 1, 2 )
x <- c( 32, 12, 6 )

You will need the vcdExtra package to use the expand.dft command. If you are using BearPortal you don’t need to install it, but if you are running the practical on your own computer you will need to install the package install.packages("vcdExtra").

Now before we can use the expand.dft() function we need to load the package containing this function. The expand.dft command allows one to convert the frequency table into a vector of samples:

library(vcdExtra)
#> Loading required package: vcd
#> Warning: package 'vcd' was built under R version 4.3.3
#> Loading required package: grid
#> Loading required package: gnm
samples <- expand.dft(data.frame(y,Frequency = x), freq = "Frequency")

Now 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 incorrect
print(out)
#> 
#>  Chi-squared test for given probabilities
#> 
#> data:  x
#> X-squared = 1.3447, df = 2, p-value = 0.5105

Actually, 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.34466
chi2_lower <- qchisq(1 - 0.01, df = 1)
print(chi2_lower)
#> [1] 6.634897

Hence, 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.