---
title: "SQA AH Statistics Exam Papers, using R"
output: html_document
date: "2026-05-26"
---

INSTRUCTIONS

Expand (or Collapse) the code chunks by clicking on the grey triangle to the right of the line number, or use the menu item: Edit > Folding > Expand All or Collapse All

Locate the code chunk you want from the Year and Question number

Run any code chunk by clicking on the green triangle 'play button' in the top right of each section of code

```{r 2016 Question 1}

heights = c(26, 15, 28, 23, 25, 15, 16, 20, 22, 27, 17, 30, 48)

# boxplot display
boxplot(x = heights,
        horizontal = TRUE,
        las = TRUE)

# stem-and-leaf display
stem(heights)

```

```{r 2016 Question 3}

n = 12
mean = 147.8
st_dev = 2.379

# create a simulated data set with correct statistics
simulated_weights = scale(1:n) * st_dev + mean

t.test(simulated_weights,
       mu = 150,
       alternative = "less")

```

```{r 2016 Question 4b}

sample_size = 20
sample_mean = 142 / sample_size
sample_st_dev = sqrt( (1120.16 - (142^2) / sample_size) / (sample_size - 1 ) )

# calculate P(X-bar > 8.1)
pnorm(q = 8.1,
      mean = sample_mean,
      sd = sample_st_dev / sqrt(sample_size),
      lower.tail = FALSE)
```

```{r 2016 Question 5(a)(i)}

yields = c(4.00, 4.50, 4.75, 4.88, 4.99, 5.10, 5.00, 5.25, 5.63)

# the t.test command delivers both hypothesis test result and confidence interval
t.test(yields,
       mu = 4.75,
       alternative = "two.sided",
       conf.level = 0.90)

```

```{r 2016 Question 8(a)}

recoveries = c(75, 65)
patients = c(100, 100)

prop.test(recoveries,
          patients,
          alternative = "greater",
          correct = FALSE) # to prevent Yate's Continuity Correction

```

```{r 2016 Question 9(b)(c)}

mu = 4

# calculate P(X > 4 + 2 * sqrt(4) ) where X ~ Po(4)
ppois(q = mu + 2 * sqrt(mu),
      lambda = mu,
      lower.tail = FALSE)

# calculate approximate P(T < 140) where T ~ Po(152)
pnorm(q = 139.5,
      mean = 152,
      sd = sqrt(152),
      lower.tail = TRUE)
```

```{r 2016 Question 10(b)}

y_values = c(1, 2)
y_probabilities = c(2/5, 3/5)

# mean of y
sum(y_values * y_probabilities)

# variance of y
sum(y_values^2 * y_probabilities) - (sum(y_values * y_probabilities))^2

```

```{r 2016 Question 11(a)(b)}

cross = c(23.5, 12.0, 21.0, 20.1, 22.0, 21.5, 22.1, 20.4, 18.3, 21.1, 21.0, 12.0)
self = c(17.4, 20.4, 20.0, 20.1, 20.0, 18.6, 18.6, 15.3, 16.5, 18.0, 18.0, 18.0)
differences = cross - self

t.test(differences,
       mu = 0,
       alternative = "greater")

# remove zeros from differences before performing the Wilxocon test
updated_differences = differences[differences != 0]

wilcox.test(x = updated_differences,
            mu = 0, # the command uses mu rather than median
            alternative = "greater")

# note that the value of V stated in the Wilcoxon test output is *not* the minimum rank sum
```

```{r 2017 Question 2}

x_values = c(20, 30, 40, 50, 60)
x_probabilities = c(1/5, 1/5, 1/5, 1/5, 1/5)

# mean of y
sum(x_values * x_probabilities)

# variance of y
sum(x_values^2 * x_probabilities) - (sum(x_values * x_probabilities))^2

```

```{r 2017 Question 3(b)}

pianists_n = 16
pianists_mean = 77
pianists_st_dev = 10

violinists_n = 14
violinists_mean = 82
violinists_st_dev = 8

# create simulated data sets with correct statistics
simulated_pianist_marks = scale(1:pianists_n) * pianists_st_dev + pianists_mean
simulated_violinist_marks = scale(1:violinists_n) * violinists_st_dev + violinists_mean

t.test(simulated_pianist_marks,
       simulated_violinist_marks,
       paired = FALSE,
       var.equal = TRUE, #this will pool the samples
       alternative = "less")

```

```{r 2017 Question 7(a)}

trees_n = 18
trees_mean = 7.46
trees_st_dev = 1.46

# create simulated data sets with correct statistics
simulated_trees_per_hectare = scale(1:trees_n) * trees_st_dev + trees_mean

# the t.test command delivers both hypothesis test result and confidence interval
t.test(simulated_trees_per_hectare,
       mu = 5.87,
       alternative = "two.sided",
       conf.level = 0.95)

```

```{r 2017 Question 8}

# (a) W ~ Po(4) ... P(2 <= W <= 6 ) = P(W <= 6) - P(W <= 1)
ppois(q = 6, lambda = 4) - ppois(q = 1, lambda = 4)
# equivalent to using the diff(erence) command...
diff(ppois(q = c(1, 6), lambda = 4))

# (b) X ~ N(4, 4) .. P(2 < X < 6 )
pnorm(q = 6, mean = 4, sd = 2) - pnorm(q = 2, mean = 4, sd = 2)
# equivalent to using the diff(erence) command...
diff(pnorm(q = c(2, 6), mean = 4, sd = 2))

# (c) Y ~ U(6, 10) .. P(8 - sqrt(4/3) < Y < 8 + sqrt(4/3))
punif(q = 8 + sqrt(4/3), min = 6, max = 10) - punif(q = 8 - sqrt(4/3), min = 6, max = 10)
# equivalent to using the diff(erence) command...
diff(punif(q = 8 + c(-1, 1) * sqrt(4/3), min = 6, max = 10))
```

```{r 2017 Question 9(c)}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

lead_n = 25
lead_mean = 174.5
lead_st_dev = 23.1

# create simulated data sets with correct statistics
simulated_lead_concentration = scale(1:lead_n) * lead_st_dev + lead_mean

z.test(x = simulated_lead_concentration,
       sigma.x = lead_st_dev,
       mu = 165.6,
       alternative = "greater")

```

```{r 2017 Question 10(c)}

table = t(data.frame(
  recaptured = c(58, 51),
  not_recaptured = c(255, 182)
))
output = chisq.test(x = table,
                    correct = FALSE) # to prevent Yate's Continuity Correction

output # displays test statistic, degrees of freedom and p-value
print("Expected frequencies:")
output$expected # displays table of expected frequencies
```

```{r 2017 Question 11(a)}

in_favour = 61
sample_size = 100

# this uses the 2-sample proportion test command syntax, but sets the second sample to have zero 'successes'. This ensures that the generated confidence interval agrees with hand-calculated answers.
prop.test(c(in_favour, 0),
          c(sample_size, sample_size),
          alternative = "two.sided",
          conf.level = 0.99,
          correct = FALSE) # to prevent Yate's Continuity Correction
```

```{r 2018 Question 1(a)}
mu = 2

# calculate P(X = 3 ) where X ~ Po(2)
dpois(x = 3,
      lambda = mu)
```

```{r 2018 Question 3(a)}

hab_1 = c(313, 342, 366, 350, 376, 438, 400)
hab_2 = c(83, 94, 91, 86, 102, 113, 98)

# here's the scatterplot, for interest
plot(x = hab_1,
     y = hab_2)

# the question requires a test on beta, but the numerical output from a test on rho is equivalent in terms of test statistic, degrees of freedom and p-value
cor.test(x = hab_1,
         y = hab_2, 
         alternative = "two.sided",
         method = "pearson")
```

```{r 2018 Question 4}

# (a) P(X > 10100) where X ~ N(10000, 250^2)
pnorm(q = 10100,
      mean = 10000,
      sd = 250,
      lower.tail = FALSE)

# (b) we want x such that P(X > x) = 0.90 X ~ N(10000, 250^2)
qnorm(p = 0.90,
      mean = 10000,
      sd = 250,
      lower.tail = FALSE)

# (c) P(T < 3000) where T ~ N(2975, 157)
pnorm(q = 3000,
      mean = 2975,
      sd = sqrt(157),
      lower.tail = TRUE)

```

```{r 2018 Question 5(c)}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

pebble_n = 100
pebble_mean = 119.4
pebble_st_dev = 21.6

# create simulated data set with correct statistics
simulated_pebble_sizes = scale(1:pebble_n) * pebble_st_dev + pebble_mean

# the output from a z.test gives the confidence interval
z.test(x = simulated_pebble_sizes,
       sigma.x = pebble_st_dev ,
       alternative = "two.sided",
       conf.level = 0.90)
```

```{r 2018 Question 6}

m = 4
n = 7

# generate all 330 combinations
all_combinations = combn(m + n, m)

# calculate the rank sums all 330 combinations
rank_sums = colSums(all_combinations)

# locate combinations whose rank sum is <= 14
combinations_needed = rank_sums <= 14

# display the combinations that meet the criteria
all_combinations[, combinations_needed]
```

```{r 2018 Question 7(a)}

number_of_failed_bulbs = c(0, 1, 2, 3, 4, 5, 6)
observed_frequencies = c(59, 38, 19, 3, 1, 0, 0)

# generate binomial model probabilities
probabilities = dbinom(x = number_of_failed_bulbs,
                       size = 6,
                       prob = 0.1)

#conduct the goodness of fit hypothesis test
output = chisq.test(x = observed_frequencies,
                    p = probabilities)

# displays expected frequencies, test statistic, degrees of freedom, p-value
print("Expected frequencies:")
output$expected
output 

# Inspection of expected frequencies shows that too many are less than 5
# We therefore need to combine categories by merging the last 5 categories

# create a function called 'combine' that combines the final n categories in a list
combine <- function(list, col) {
  temp_col = length(list) - col # the number of columns lefts unchanged
  temp_start = list[1:temp_col] # capture the first columns unchanged
  temp_end = sum(list[(temp_col + 1):length(list)]) # sum the remaining columns
  c(temp_start, temp_end) # splice the two lists together
}

# combine categories before applying goodness of fit test
categories_to_combine = 5
new_observed_frequencies = combine(observed_frequencies, categories_to_combine)
new_probabilities = combine(probabilities, categories_to_combine)

#conduct the goodness of fit hypothesis test with the combined category data
output = chisq.test(x = new_observed_frequencies,
                    p = new_probabilities)

# displays expected frequencies, test statistic, degrees of freedom, p-value
print("Expected frequencies:")
output$expected
output 

```

```{r 2018 Question 8(a)}

accidents = 18
sample_size = 100

# the output of this test gives the correct p-value. All other output can be ignored.
prop.test(x = accidents,
          n = sample_size,
          p = 0.119,
          alternative = "greater",
          correct = FALSE) # to prevent Yate's Continuity Correction

```

```{r 2018 Question 11(b)}

welfare = c(19, 22, 39, 42, 59, 70, 79, 83, 117, 119, 140, 142, 144, 193, 204, 210, 236, 290)
employment = c(64.7, 51.1, 47.0, 54.0, 42.6, 48.3, 36.2, 47.9, 64.7, 54.9, 58.6, 69.9, 54.9, 61.9, 69.3, 76.0, 74.1, 78.7)

# fit the linear model and store it in a variable called 'model'
model = lm(employment ~ welfare)

# display linear model
model

# display fitted value and residual for 5th piece of data (for Czech Republic)
model$fitted.values[5]
model$residuals[5]

```

```{r 2018 Question 12}

location_A = c(13.9, 11.0, 11.3, 12.5, 13.4, 11.5, 14.5, 13.4, 12.2, 11.3)
location_B = c(14.4, 12.2, 17.4, 12.0, 16.1, 10.9, 19.4, 10.3)

t.test(location_A,
       location_B,
       paired = FALSE,
       var.equal = TRUE, #this will pool the samples
       alternative = "two.sided")

```

```{r 2019 Question 2(c)}

table = t(data.frame(
  colour_blind_no = c(450, 506),
  colour_blind_yes = c(40, 4)
))
output = chisq.test(x = table,
                    correct = FALSE) # to prevent Yate's Continuity Correction

output # displays test statistic, degrees of freedom and p-value
print("Expected frequencies:")
output$expected # displays table of expected frequencies

```

```{r 2019 Question 6}

energies = c(2725, 2650, 2421, 2793, 2239, 3225, 2156, 2692, 2369, 2725)

# the t.test command delivers both hypothesis test result and confidence interval
t.test(energies,
       alternative = "two.sided",
       conf.level = 0.95)

t.test(energies,
       alternative = "two.sided",
       conf.level = 0.99)

```

```{r 2019 Question 8(b)(c)}

Sgg = 531.5676
Sgc = 555.0811
Scc = 1731.2973

coeff_determination = Sgc^2 / (Sgg * Scc)
coeff_determination

r = sqrt(coeff_determination)
t = r * sqrt(37 - 2) / sqrt(1 - r^2)

p_value = 2 * pt(q = t,
                 df = 37 - 2,
                 lower.tail = FALSE)
p_value

```

```{r 2019 Question 9}

marks = c(86, 80, 78, 73, 69, 65, 62, 61, 59, 58, 54, 51, 49, 47, 43, 40, 38, 37, 35, 32, 29, 29)
differences = marks - 65

# remove zeros from differences before performing the Wilxocon test
updated_differences = differences[differences != 0]

wilcox.test(x = updated_differences,
            mu = 0, # the command uses mu rather than median
            alternative = "less")

```

```{r 2019 Question 10(a)}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

wingspan_n = 25
wingspan_mean = 48.3
wingspan_st_dev = 4

# create simulated data set with correct statistics
simulated_wingspans = scale(1:wingspan_n) * wingspan_st_dev + wingspan_mean

z.test(x = simulated_wingspans,
       sigma.x = wingspan_st_dev,
       mu = 50,
       alternative = "less")

```

```{r 2021 Paper 2 Question 1}

s_values = c(2, 4, 6, 8, 10)
s_probabilities = s_values / sum(s_values)

# mean of y
sum(s_values * s_probabilities)

# variance of y
sum(s_values^2 * s_probabilities) - (sum(s_values * s_probabilities))^2
```

```{r 2021 Paper 2 Question 3}

mites = c(6, 8, 11, 13, 6, 14, 11, 9, 6, 7, 11, 8, 6, 14)
differences = mites - 7

# remove zeros from differences before performing the Wilxocon test
updated_differences = differences[differences != 0]

wilcox.test(x = updated_differences,
            mu = 0, # the command uses mu rather than median
            alternative = "greater")

# note that the value of V stated in the Wilcoxon test output is *not* the minimum rank sum

```

```{r 2021 Paper 2 Question 5}

brand_A_n = 11
brand_A_mean = 54
brand_A_st_dev = 5

brand_B_n = 15
brand_B_mean = 47
brand_B_st_dev = 11

# create simulated data sets with correct statistics
simulated_brand_A_times = scale(1:brand_A_n) * brand_A_st_dev + brand_A_mean
simulated_brand_B_times = scale(1:brand_B_n) * brand_B_st_dev + brand_B_mean

t.test(simulated_brand_A_times,
       simulated_brand_B_times,
       paired = FALSE,
       var.equal = TRUE, #this will pool the samples
       alternative = "two.sided")

```

```{r 2021 Paper 2 Question 6(c)}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

time_n = 25
time_mean = 409
time_st_dev = 130

# create simulated data set with correct statistics
simulated_times = scale(1:time_n) * time_st_dev + time_mean

# the output from a z.test gives the confidence interval
z.test(x = simulated_times,
       sigma.x = time_st_dev ,
       alternative = "two.sided",
       conf.level = 0.95)

```

```{r 2021 Paper 2 Question 8}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

time_n = 50
time_mean = 16.1
time_st_dev = 2

# create simulated data set with correct statistics
simulated_times = scale(1:time_n) * time_st_dev + time_mean

# the output from a z.test gives the confidence interval
z.test(x = simulated_times,
       sigma.x = time_st_dev,
       mu = 15,
       alternative = "greater",
       conf.level = 0.95)
```

```{r 2021 Paper 2 Question 10(a)}

selling_beyond = 13
sample_size = 50

# this uses the 2-sample proportion test command syntax, but sets the second sample to have zero 'successes'. This ensures that the generated confidence interval agrees with hand-calculated answers.
prop.test(c(selling_beyond, 0),
          c(sample_size, sample_size),
          alternative = "two.sided",
          conf.level = 0.95,
          correct = FALSE) # to prevent Yate's Continuity Correction

```

```{r 2021 Paper 2 Question 11(b)}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

group_B_n = 70
group_B_mean = 55.4
group_B_st_dev = 10.08

group_C_n = 60
group_C_mean = 51.8
group_C_st_dev = 10.49

# create simulated data sets with correct statistics
simulated_group_B_scores = scale(1:group_B_n) * group_B_st_dev + group_B_mean
simulated_group_C_scores = scale(1:group_C_n) * group_C_st_dev + group_C_mean

z.test(x = simulated_group_B_scores,
       y = simulated_group_C_scores,
       sigma.x = group_B_st_dev,
       sigma.y = group_C_st_dev,
       alternative = "two.sided")

# for interest, by comparison, here is what the t-test would have been
t.test(simulated_group_B_scores,
       simulated_group_C_scores,
       paired = FALSE,
       var.equal = TRUE, #this will pool the samples
       alternative = "two.sided")

```

```{r 2022 Paper 2 Question 1}

table = t(data.frame(
  infected_yes = c(76, 129),
  infected_no = c(399, 332)
))
output = chisq.test(x = table,
                    correct = FALSE) # to prevent Yate's Continuity Correction

output # displays test statistic, degrees of freedom and p-value
print("Expected frequencies:")
output$expected # displays table of expected frequencies

```

```{r 2022 Paper 2 Question 2}

# (b) P(X = 0) where X ~ Po(2.3)
dpois(x = 0,
      lambda = 2.3)

# (c) P(X = 2 and Y = 2) where X ~ Po(2.3) and Y ~ Po(1.7)
dpois(x = 2, lambda = 2.3) * dpois(x = 2, lambda = 1.7)

# (d) P(W > 5) where W ~ Po(4)
ppois(q = 5,
      lambda = 4,
      lower.tail = FALSE)
```

```{r 2022 Paper 2 Question 5(b)}

french = c(67, 83, 71, 59, 49, 89, 42, 55, 77)
german = c(64, 82, 71, 62, 42, 85, 39, 50, 75)

differences = french - german

t.test(differences,
       mu = 0,
       alternative = "two.sided")

```

```{r 2022 Paper 2 Question 7(b)}

successes = 14
sample_size = 50

# this uses the 2-sample proportion test command syntax, but sets the second sample to have zero 'successes'. This ensures that the generated confidence interval agrees with hand-calculated answers.
prop.test(c(successes, 0),
          c(sample_size, sample_size),
          alternative = "two.sided",
          conf.level = 0.99,
          correct = FALSE) # to prevent Yate's Continuity Correction

```

```{r 2022 Paper 2 Question 9(b)}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

width_n = 45
width_mean = 52.6
width_st_dev = sqrt(103.25)

# create simulated data sets with correct statistics
simulated_widths = scale(1:width_n) * width_st_dev + width_mean

z.test(x = simulated_widths,
       sigma.x = width_st_dev,
       mu = 50,
       alternative = "greater")

```

```{r 2022 Paper 2 Question 10(a)}

Sxx = 278.61
Syy = 10.95
Sxy = 46.29
n = 6

r = Sxy / sqrt(Sxx * Syy)
r

t = r * sqrt(n - 2) / sqrt(1 - r^2)
t

p_value = 2 * pt(q = t,
                 df = n - 2,
                 lower.tail = FALSE)
p_value

```

```{r 2022 Paper 2 Question 11(a)}

adults = c(1.3, 2.2, 1.5, 3.5, 0.3, 2.7, 3.5, 2.3, 2.9, 4.0)
juveniles = c(1.1, 4.1, 1.7, 1.3, 0.7, 1.9, 2.9, 1.1, 2.8, 0.9)

# back-to-back stem and leaf diagrams are not built into R by default
# the 'aplpack' package does support them, but technical issues may need to be overcome
stem(adults)
stem(juveniles)

# the process for a Mann-Whitney test is embedded within the wilcox.test command
wilcox.test(x = adults,
            y = juveniles,
            paired = FALSE,
            alternative = "two.sided")
# note that the value of 'W' in the output is *not* the rank sum used in the AH Stats course
```

```{r 2023 Paper 2 Question 2}

steps = c(320, 310, 321, 304, 298, 328, 296, 307, 314, 295)
differences = steps - 300

# remove zeros from differences before performing the Wilxocon test
updated_differences = differences[differences != 0]

wilcox.test(x = updated_differences,
            mu = 0, # the command uses mu rather than median
            alternative = "greater")

# note that the value of V stated in the Wilcoxon test output is *not* the minimum rank sum
```

```{r 2023 Paper 2 Question 4}

observed_offspring = c(78, 90, 152)

# generate probabilities in the ratio 1:1:2
probabilities = c(1/4, 1/4, 2/4)

#conduct the goodness of fit hypothesis test
output = chisq.test(x = observed_offspring,
                    p = probabilities)

# displays expected frequencies, test statistic, degrees of freedom, p-value
print("Expected frequencies:")
output$expected
output 

```

```{r 2023 Paper 2 Question 6}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

length_n = 75
length_mean = 3840 / 75
length_st_dev = sqrt((198240 - 3840^2/75) / (75 - 1))

# create simulated data sets with correct statistics
simulated_lengths = scale(1:length_n) * length_st_dev + length_mean

z.test(x = simulated_lengths,
       sigma.x = length_st_dev,
       mu = 50,
       alternative = "greater")

# for interest, a single sample t-test would have given...
t.test(x = simulated_lengths,
       mu = 50,
       alternative = "greater")

```

```{r 2023 Paper 2 Question 8}

n = 25
r = 0.652

t = r * sqrt(n - 2) / sqrt(1 - r^2)
t

p_value = 2 * pt(q = t,
                 df = n - 2,
                 lower.tail = FALSE)
p_value

```

```{r 2023 Paper 2 Question 9(a)}

with_tracker = c(5.1, 10, 10.8, 7.5, 6.2, 10.2, 5.4, 4.2, 8.1, 11.1, 10.2, 5.3)
without_tracker = c(4, 9.5, 12, 5.5, 5.9, 11, 4.8, 3.5, 6.5, 11.5, 9.4, 5.1)

differences = with_tracker - without_tracker

t.test(differences,
       mu = 0,
       alternative = "greater")

```

```{r 2023 Paper 2 Question 10}

homeless = 23312
sample_size = 37878

# the output of this test gives the correct p-value. All other output can be ignored.
prop.test(x = homeless,
          n = sample_size,
          p = 0.624,
          alternative = "two.sided",
          correct = FALSE) # to prevent Yate's Continuity Correction

```

```{r 2023 Paper 2 Question 11}

# set up a system of two simultaneous equations in the matrix form Ax = B
A = matrix(data = c(1, 1,
                    qnorm(0.1), qnorm(0.95)),
           nrow = 2,
           ncol = 2)
B = c(17, 24)

# solve the system of equations to give mu and sigma
solve(A, B)

```

```{r 2023 Paper 2 Question 12(a)}

in_favour = 55
sample_size = 100

# this uses the 2-sample proportion test command syntax, but sets the second sample to have zero 'successes'. This ensures that the generated confidence interval agrees with hand-calculated answers.
prop.test(c(in_favour, 0),
          c(sample_size, sample_size),
          alternative = "two.sided",
          conf.level = 0.99,
          correct = FALSE) # to prevent Yate's Continuity Correction

```


```{r 2024 Paper 1 Question 1(c)}

table = t(data.frame(
  decade_1980s = c(3, 8, 1, 21, 17),
  decade_1990s = c(2, 7, 5, 20, 16),
  decade_2000s = c(2, 5, 8, 31, 4),
  decade_2010s = c(2, 3, 14, 26, 5)
))
output = chisq.test(x = table,
                    correct = FALSE) # to prevent Yate's Continuity Correction

output # displays test statistic, degrees of freedom and p-value
print("Expected frequencies:")
output$expected # displays table of expected frequencies

```

```{r 2024 Paper 1 Question 1(d)(i)}

age_rating_18 = c(37, 5)
total_all_films = c(150, 50)

p = sum(age_rating_18) / sum(total_all_films) # pooled proportion
p # displays pooled proportion

z = diff(age_rating_18 / total_all_films) / sqrt (p * (1-p) * sum(1 / total_all_films))
z # displays test statistic
```

```{r 2024 Paper 1 Question 2(b)}

standard = c(7500, 8000, 2000, 550, 1250, 1000, 2250, 6800, 3400, 6300, 9100, 970, 1040, 670, 400)
new = c(410, 250, 800, 1400, 7900, 7400, 1020, 6000, 920, 1420, 2700, 4200, 5200, 4100)

# the process for a Mann-Whitney test is embedded within the wilcox.test command
output = wilcox.test(x = standard,
                     y = new,
                     paired = FALSE,
                     alternative = "two.sided")
output # displays Mann-Whitney output

# note that the value of 'W' in the output is *not* the rank sum used in the AH Statistics course
W_output = output$statistic

# the value of 'W' in the standard R output is actually the distance away from the smallest possible rank sum
# the smallest possible rank sum here is 0.5 * m * (m + 1) where m = smallest sample size
m = min(length(new), length(standard))
n = max(length(new), length(standard))

# now refer to the top of page 16 in Statistical Formulae and Tables booklet
Wm = 0.5 * m * (m + 1) + W_output
Wm_reversed_ranks = m * (m + n + 1) - Wm 
W = min(Wm, Wm_reversed_ranks)
W # this is the value of W used in the AH Statistics course
```

```{r 2024 Paper 1 Question 2(c)}

month_3 = c(7500, 8000, 2000, 550, 1250, 1000, 2250, 6800, 3400, 6300, 9100, 970, 1040, 670, 400)
month_6 = c(7480, 8010, 2010, 560, 1210, 980, 2180, 6790, 3370, 6260, 9060, 960, 990, 640, 430)
differences = month_3 - month_6
differences # display differences, and check there are no zeros

ranks = rank(abs(differences))
ranks # display ranks of absolute values of differences

wilcox.test(x = differences,
            mu = 0, # the command uses mu rather than median
            alternative = "two.sided",
            exact = FALSE) # approximate p-value

# note that the value of V stated in the Wilcoxon test output is *not* the minimum rank sum

W_negative = sum((differences < 0) * ranks)
W_negative # display rank sum of negative differences

W_positive = sum((differences > 0) * ranks)
W_positive # display rank sum of positive differences

W = min(W_negative, W_positive)
W # rank sum used in AH Statistics course
```

```{r 2024 Paper 2 Question 1}

wild_birds = c(11, 27, 45, 63, 65, 70, 77, 79, 87, 88, 90, 95, 102, 130)

# boxplot display
boxplot(x = wild_birds,
        horizontal = TRUE,
        las = TRUE)

# the boxplot shows one possible outlier below the lower fence
```

```{r 2024 Paper 2 Question 3(a)}

# calculate P(X = 9) where X ~ B(12, 0.88)
dbinom(x = 9,
       size = 12, # number of trials
       prob = 0.88) # probability of success

```

```{r 2024 Paper 2 Question 3(b)}

# calculate P(X > 36) using Normal Approximation
n = 48
p = 0.88
normal_mean = n * p
normal_st_dev = sqrt(n * p * (1 - p))

pnorm(q = 36.5, # using continuity correction
      mean = normal_mean,
      sd = normal_st_dev,
      lower.tail = FALSE) # so it calculates P(Z > ...) and not P(Z < ...)

# for interesting comparison, calculating P(X > 36) exactly, where X ~ B(48, 0.88)
sum(dbinom(x = c(37:48), # all values of x from 37 to 48 inclusive
           size = 48, # number of trials
           prob = 0.88)) # probability of success

```

```{r 2024 Paper 2 Question 4}

observed_grades = c(185, 170, 197, 163, 155)

# generate probabilities according to discrete uniform distribution, U(5)
probabilities = c(1/5, 1/5, 1/5, 1/5, 1/5)

#conduct the goodness of fit hypothesis test
output = chisq.test(x = observed_grades,
                    p = probabilities)

# displays expected frequencies, test statistic, degrees of freedom, p-value
print("Expected frequencies:")
output$expected
output 

```

```{r 2024 Paper 2 Question 5(b)}

attempt = c(2, 2, 3, 4, 4, 5, 9, 12, 14)
errors = c(11, 10, 7, 6, 5, 4, 2, 1, 1)

# generate scatterplot of original data
plot(x = attempt,
     y = errors)

# generate scatterplot of transformed data
plot(x = attempt,
     y = 1 / errors)

# fit the linear model to transformed data and store it in a variable called 'model'
model = lm((1 / errors) ~ attempt)

# display linear model
model

# calculate prediction interval for rat completing maze on attempt number 7
prediction_interval = predict.lm(model,
                                 newdata = data.frame(attempt = 7),
                                 interval = "prediction",
                                 level = 0.95)

prediction_interval # display prediction interval for 1/errors, for transformed data model

1 / prediction_interval # convert prediction interval from '1/errors' to 'errors'
```

```{r 2024 Paper 2 Question 6(a)}

work_properly = 14
sample_size = 20

# this uses the 2-sample proportion test command syntax, but sets the second sample to have zero 'successes'. This ensures that the generated confidence interval agrees with hand-calculated answers.
prop.test(x = c(work_properly, 0),
          n = c(sample_size, sample_size),
          alternative = "two.sided",
          conf.level = 0.99,
          correct = FALSE) # to prevent Yate's Continuity Correction

```

```{r 2024 Paper 2 Question 8}

area_1_n = 9
area_1_mean = 38.7
area_1_st_dev = 3.42

area_2_n = 11
area_2_mean = 36.9
area_2_st_dev = 3.51

# create simulated data sets with correct statistics
simulated_area_1 = scale(1:area_1_n) * area_1_st_dev + area_1_mean
simulated_area_2 = scale(1:area_2_n) * area_2_st_dev + area_2_mean

t.test(simulated_area_1,
       simulated_area_2,
       paired = FALSE, # non-paired data
       var.equal = TRUE, # this will pool the samples
       alternative = "greater")

```

```{r 2024 Paper 2 Question 9(b)}

mean_rate = 6000 / 2400

# calculate P(X = 0) where X ~ Po(2.5)
dpois(x = 0,
      lambda = mean_rate)

```

```{r 2024 Paper 2 Question 9(c)}

# calculate P(X <= x) where X ~ Po(2.5) and x takes values from 0 to 7
# and look for the first to exceed 0.90
ppois(q = 0:7,
      lambda = 2.5)

# and to check you have the correct values either side of this threshold....

# calculate P(X <= 4) where X ~ Po(2.5)
ppois(q = 4,
      lambda = 2.5)

# calculate P(X <= 5) where X ~ Po(2.5)
ppois(q = 5,
      lambda = 2.5)

```

```{r 2024 Paper 2 Question 10}

z_95 = qnorm(p = 0.95) # value of z such that P(Z < z) = 0.95

# we need to solve for the minimum value of integer n the inequality:
# z_95 * 2.9 / sqrt (n) < 0.7

# we shall use 'brute force' by substituting in lots of values of n, to then find the minimum value of n that works

# first, set up a list of possible values of n, from 1 to 100
n_possible = 1:100

# for each possible value of n, return TRUE or FALSE in a list, depending on whether it satisfies the inequality
evaluations = z_95 * 2.9 / sqrt (n_possible) < 0.7

# extract the minimum value of n that gave a TRUE
min(n_possible[evaluations])

```

```{r 2024 Paper 2 Question 12(a)}

# calculate P(total mass > 25000) using combination of J ~ N(69, 6) and H ~ N(453, 16)
normal_mean = 48 * 69 + 48 * 453
normal_mean

normal_variance = 48 * 6 + 48 * 16
normal_variance

pnorm(q = 25000,
      mean = normal_mean,
      sd = sqrt(normal_variance),
      lower.tail = FALSE)

```

```{r 2024 Paper 2 Question 12(b)}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

mass_n = 10
mass_mean = 527.5
mass_st_dev = 5

# create simulated data sets with correct statistics
simulated_masses = scale(1:mass_n) * mass_st_dev + mass_mean

z.test(x = simulated_masses,
       sigma.x = mass_st_dev,
       mu = 522,
       alternative = "two.sided")

```

```{r 2024 Paper 2 Question 13(a)}

sample_size = 5
sample_mean = 102
sample_st_dev = 0.13
sigma_limits = 2

# calculate both lower and upper limits
sample_mean + c(-1, 1) * sigma_limits * sample_st_dev / sqrt(sample_size)

```

```{r 2025 Paper 1 Question 1(c)}

# x values
xi = c(0:7)

# observed frequencies
Oi = c(0, 3, 6, 6, 5, 5, 2, 1)

# mean rate
sum(Oi * xi) / sum(Oi)
```

```{r 2025 Paper 1 Question 1(d)}

# observed frequencies
Oi = c(0, 3, 6, 6, 5, 8)

# expected frequency using P(X >= 5) = P(X > 4) where X ~ Po(3.5)
sum(Oi) * ppois(q = 4,
                lambda = 3.5,
                lower.tail = FALSE)
```

```{r 2025 Paper 1 Question 1(f)(ii)}

# observed frequencies
Oi = c(3, 6, 6, 5, 8)

# grouped probabilities
Pi = c(ppois(q = 1, lambda = 3.5), # for P(X <= 1)
       dpois(x = 2, lambda = 3.5), # for P(X = 2)
       dpois(x = 3, lambda = 3.5), # for P(X = 3)
       dpois(x = 4, lambda = 3.5), # for P(X = 4)
       ppois(q = 4, lambda = 3.5, lower.tail = FALSE)) # for P(X > 4)

# calculate and display expected frequencies
Ei = sum(Oi) * Pi
Ei

# calculate test statistic, X^2
sum( ((Oi - Ei) ^ 2) / Ei)

```

```{r 2025 Paper 2 Question 1}

x_values = c(1:6)
x_probabilities = c(1/60, 1/6, 1/6, 1/6, 1/6, 19/60)

# calculate and display the mean of X
mean_x = sum(x_values * x_probabilities)
mean_x

# calculate mean of X^2
mean_x2 = sum(x_values^2 * x_probabilities)

# calculate and display the variance of X
variance_x = mean_x2 - (mean_x)^2
variance_x

# variance of y
variance_y = (8^2 - 1) / 12

# variance of X - Y = V(X - Y) = V(X) + V(Y)
variance_x + variance_y

# standard deviation, SD(X - Y)
sqrt(variance_x + variance_y)

```

```{r 2025 Paper 2 Question 2(a)}

distances = c(1, 4, 4, 7, 7, 8, 9, 13, 15, 20, 21, 22, 25, 27, 46)

# boxplot display which shows one possible outlier above the upper fence
boxplot(x = distances,
        horizontal = TRUE)

# calculate Q1 and Q3
quartiles = quantile(x = distances,
                     probs = c(0.25, 0.75),
                     type = 6) # to select the method used by SQA

lower_fence = quartiles[1] - 1.5 * diff(quartiles)

upper_fence = quartiles[2] + 1.5 * diff(quartiles)

# detect if any data is below the lower fence (TRUE or FALSE)
min(distances) < lower_fence

# detect if any data is above the upper fence (TRUE or FALSE)
max(distances) > upper_fence
```

```{r 2025 Paper 2 Question 2(b)}
table = t(data.frame(
  retail = c(36, 26),
  food = c(33, 5)
))

# perform test
output = chisq.test(x = table,
                    correct = FALSE) # to prevent Yate's Continuity Correction

# displays test statistic, degrees of freedom and p-value
output

# display expected frequencies that were used
print("Expected frequencies:")
output$expected
```

```{r 2025 Paper 2 Question 3}

# calculate approximate P(X > 20) where X ~ Po(14)
pnorm(q = 20.5,
      mean = 14,
      sd = sqrt(14),
      lower.tail = FALSE)

# for interest, calculate exact P(X > 20) where X ~ Po(14)
ppois(q = 20,
      lambda = 14,
      lower.tail = FALSE)

```

```{r 2025 Paper 2 Question 4}

# calculate P(X-bar > 5.50) where X-bar ~ N(4.70, 2.80/sqrt(50))
pnorm(q = 5.50,
      mean = 4.70,
      sd = 2.80 / sqrt(50),
      lower.tail = FALSE)

```

```{r 2025 Paper 2 Question 5}

data = c(40.1, 40.3, 39.5, 40.6, 40.8, 39.9, 40.0, 41.3, 38.3)
differences = data - 40

# remove zeros from differences before performing the Wilcoxon test
updated_differences = differences[differences != 0]

# perform Wilcoxon test
wilcox.test(x = updated_differences,
            mu = 0, # the command uses mu rather than median
            alternative = "two.sided",
            exact = FALSE) # due to tied ranks, p-value will be approximate

# note that the value of V stated in the Wilcoxon test output is *not* the minimum rank sum

# calculate and display ranks of absolute values of differences
ranks = rank(abs(updated_differences))
ranks 

# calculate and display rank sum of negative differences
W_negative = sum((updated_differences < 0) * ranks)
W_negative 

# calculate and display rank sum of positive differences
W_positive = sum((updated_differences > 0) * ranks)
W_positive 

# calculate and display the rank sum used in AH Statistics course
W = min(W_negative, W_positive)
W 
```

```{r 2025 Paper 2 Question 6}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

years_00_09_mean = 3.9868
years_00_09_st_dev = 1.3396
years_00_09_n = 990

years_10_18_mean = 4.2510
years_10_18_st_dev = 1.2342
years_10_18_n = 735

# create simulated data sets with correct statistics
simulated_years_00_09_values = scale(1:years_00_09_n) * years_00_09_st_dev + years_00_09_mean
simulated_years_10_18_values = scale(1:years_10_18_n) * years_10_18_st_dev + years_10_18_mean

# perform the z-test
z.test(x = simulated_years_00_09_values,
       y = simulated_years_10_18_values,
       sigma.x = years_00_09_st_dev,
       sigma.y = years_10_18_st_dev,
       alternative = "less")

# for interest, by comparison, here is what the t-test would have been
t.test(simulated_years_00_09_values,
       simulated_years_10_18_values,
       paired = FALSE,
       var.equal = TRUE, #this will pool the samples
       alternative = "less")

```

```{r 2025 Paper 2 Question 8}

# record pass rates, for Centre A then Centre B
successes = c(14, 21)
trials = c(29, 30)

# perform two-sample proportion test
prop.test(x = successes,
          n = trials,
          alternative = "less",
          correct = FALSE) # to prevent continuity correction

# you will see from the above test output the mention of X-squared, which stems from a lot of parallels between two-sample proportion tests and a chi-squared test of association

# the next lines of code perform the same test, but manually.....

# calculate and display the pooled proportion
p = sum(successes) / sum(trials) 
p 

# calculate and display the test statistic
z = -diff(successes / trials) / sqrt (p * (1-p) * sum(1 / trials))
z 

# calculate the p-Value of the test
pnorm(q = z,
      mean = 0,
      sd = 1)

```

```{r 2025 Paper 2 Question 9}

y_mean = 0.254
y_st_dev = 2.54 * 0.25

# calculate the probability of within 1 of 0, P(-1 < Y < 1) = P(Y < 1) - P(Y < -1)
diff(pnorm(q = c(-1,1),
           mean = y_mean,
           sd = y_st_dev))

# calculate expected number of planks
80 * diff(pnorm(q = c(-1,1),
           mean = y_mean,
           sd = y_st_dev))

```

```{r 2025 Paper 2 Question 10}

sigma_x = 27.9
sigma_x2 = 130.2
sample_n = 6

# calculate sample mean and sample standard deviation from summary statistics
sample_mean = sigma_x / sample_n
sample_st_dev = sqrt((sigma_x2 - (sigma_x ^ 2) / sample_n) / (sample_n - 1))

# create a simulated data set that has the correct mean and standard deviation
simulated_masses = scale(1:sample_n) * sample_st_dev + sample_mean

# perform t-test
t.test(x = simulated_masses,
       mu = 4.5,
       alternative = "two.sided")

```

```{r 2025 Paper 2 Question 12}

sample_n = 15
sample_mean = 139.5
sample_st_dev = 0.7

# create a simulated data set that has the correct mean and standard deviation
simulated_heights = scale(1:sample_n) * sample_st_dev + sample_mean

# use the t-test function (without specifying a value for mu) to deliver a 98% confidence interval
t.test(x = simulated_heights,
       conf.level = 0.98)

```

```{r 2026 Paper 1 Question 1b}

table = t(data.frame(
  england = c(7, 12, 8, 4, 1),
  ireland = c(10, 6, 6, 5, 6),
  scotland = c(14, 9, 3, 4, 2),
  wales = c(13, 4, 6, 5, 7)
))

# perform test
output = chisq.test(x = table,
                    correct = FALSE) # to prevent Yate's Continuity Correction

# display expected frequencies that were used
print("Expected frequencies:")
output$expected

```

```{r 2026 Paper 1 Question 2e}

# setup variable values
m = 18
n = 34
W = 404

# calculate mean and variance of normal approximation to W
E_W = m * (m + n + 1) / 2
V_W = m * n * (m + n + 1) / 12

# calculate p-value
pnorm(q = W + 0.5,
      mean = E_W,
      sd = sqrt(V_W))

```

```{r 2026 Paper 2 Question 1b}

# the package called 'BSDA' is needed for the z.test function
if(!require(BSDA)){install.packages("BSDA"); library(BSDA)}

white_mean = 20.2
white_st_dev = 1.38
white_n = 216

asian_mean = 19.8
asian_st_dev = 1.15
asian_n = 87

# create simulated data sets with correct statistics
simulated_white_values = scale(1:white_n) * white_st_dev + white_mean
simulated_asian_values = scale(1:asian_n) * asian_st_dev + asian_mean

# perform the z-test
z.test(x = simulated_white_values,
       y = simulated_asian_values,
       sigma.x = white_st_dev,
       sigma.y = asian_st_dev,
       alternative = "two.sided")

```

```{r 2026 Paper 2 Question 4}

accidental = c(5, 9, 7, 8, 4, 8, 9)
deliberate = c(38, 30, 32, 30, 36, 40, 44)
differences = deliberate - 6 * accidental

# remove any zeros from differences before performing the Wilxocon test
updated_differences = differences[differences != 0]

wilcox.test(x = updated_differences,
            mu = 0, # the command uses mu rather than median
            exact = FALSE, # do not calculate an exact p-value
            alternative = "greater")

```

```{r 2026 Paper 2 Question 6bii}

sample_n = 15
sample_mean = 67
sample_st_dev = 4

# create a simulated data set that has the correct mean and standard deviation
simulated_masses = scale(1:sample_n) * sample_st_dev + sample_mean

# perform t-test
t.test(x = simulated_masses,
       mu = 70,
       alternative = "less")

```

```{r 2026 Paper 2 Question 7}

# 7a) calculate P(X = 8) where X ~ Po(7.5)
dpois(x = 8,
      lambda = 7.5)

# 7c) calculate P(Y > 15) using normal approximation and continuity correction, where Y ~ Po(20)
pnorm(q = 15.5,
      mean = 20,
      sd = sqrt(20),
      lower.tail = FALSE)

```

```{r 2026 Paper 2 Question 8}

# record experimental results for 18-month and 14-month olds
successes = c(38, 6)
trials = c(55, 31)

# perform two-sample proportion test
prop.test(x = successes,
          n = trials,
          alternative = "greater",
          correct = FALSE) # to prevent continuity correction

# you will see from the above test output the mention of X-squared, which stems from a lot of parallels between two-sample proportion tests and a chi-squared test of association

# the next lines of code perform the same test, but manually.....

# calculate and display the pooled proportion
p = sum(successes) / sum(trials) 
p 

# calculate and display the test statistic
z = -diff(successes / trials) / sqrt (p * (1-p) * sum(1 / trials))
z 

# calculate the p-Value of the test
pnorm(q = z,
      mean = 0,
      sd = 1,
      lower.tail = FALSE)

```

```{r 2026 Paper 2 Question 11}

# 11a) calculate P(X >= 1) = P(X > 0), where X ~ B(100, 0.043)
pbinom(q = 0,
       size = 100,
       prob = 0.043,
       lower.tail = FALSE)

# 11c) calculating a 95% confidence interval

late_trains = 6
sample_size = 75

# this uses the 2-sample proportion test command syntax, but sets the second sample to have zero 'successes'. This ensures that the generated confidence interval agrees with hand-calculated answers.
prop.test(c(late_trains, 0),
          c(sample_size, sample_size),
          alternative = "two.sided",
          conf.level = 0.95,
          correct = FALSE) # to prevent Yate's Continuity Correction

```