Lab — Data Science with R

18 practicals, each set out as 1. Question, 2. Aim, 3. Steps, 4. Programme, 5. Execution and Results.

Two versions of every experiment

Every R script is run, and under 5. Execution and Results is what it printed, with the plots it drew. Two cannot simply be run in a terminal: 17_plotly.R builds charts meant for RStudio's Viewer, and 18_shiny_app.R is a web app. Their drivers, _drive_17_plotly.py and _drive_18_shiny_app.py, open them in Chromium, use them as you would, and take the screenshots.

Until October 2026 R could not be installed where this material is verified, and these scripts were desk-checked only, with the numbers in their comments taken from the Python equivalents. Running them found four comments whose numbers R does not give, in experiments 4, 10, 13 and 16; each is corrected, with a note in the file and on this page.

What the Python side buys you: the same calculation, done a second way. When 04_regression.R says the slope is 4.3030, R prints it, and 04_regression.py computes it again from the sums, asserts it, and cross-checks it against Statistical Foundations for Data Science Unit 4, where the same data was worked by hand.

python3 tools/data-science/run_r_equivalents.py

That runs all 14 Python equivalents and 17 of the R scripts (the Shiny app is a server, which capture_lab_outputs.py runs instead), and checks all 18 for balanced delimiters.

The experiments

Experiments 7, 12, 17 and 18 have no Python equivalent: they demonstrate R syntax, ggplot2's grammar, plotly's R interface and the Shiny framework. A translation would teach nothing.

Experiment 1 — Mean, median, mode, variance, SD

1. Question

Find the mean, median, mode, variance and standard deviation of twenty students' marks.

2. Aim

Describe the centre and the spread of a set of marks in R, and know which variance R gives.

3. Steps

In R, 01_descriptive.R:

1. Enter the marks.
2. Find the centre: mean, median and mode.
3. Measure the spread: variance, SD, range and quartiles.
4. Convert to the population variance.
5. Summarise everything at once.

In Python, python/01_descriptive.py:

1. Compute the centre and the spread.
2. Print them, against the R function for each.
3. Check them against the statistics module.

4. Programme

In R, 01_descriptive.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 1: Mean, Median, Mode, Variance, Standard Deviation
# Python equivalent: python/01_descriptive.py

# Step 1: Enter the marks
marks <- c(45, 67, 78, 52, 89, 91, 73, 64, 58, 82,
           76, 69, 71, 85, 60, 55, 93, 48, 79, 66)

# --- Central tendency ---
# Step 2: Find the centre: mean, median and mode
mean(marks)      # 70.05
median(marks)    # 70.00

# R has NO built-in mode() for the statistical mode -- mode() reports the
# storage type. Define one, which the syllabus expects you to know:
statistical_mode <- function(v) {
  freq <- table(v)
  as.numeric(names(freq)[freq == max(freq)])
}
statistical_mode(marks)   # every value occurs once here, so all 20 are returned

# --- Dispersion ---
# Step 3: Measure the spread: variance, SD, range and quartiles
var(marks)       # 202.8921   <- R's var() divides by n-1 (SAMPLE)
sd(marks)        #  14.2440   <- likewise
range(marks)     # 45 93
diff(range(marks))  # 48
IQR(marks)
quantile(marks)

# Step 4: Convert to the population variance
# Population variance, if you need it, must be computed explicitly:
n <- length(marks)
var(marks) * (n - 1) / n     # 192.7475
sqrt(var(marks) * (n - 1) / n)   # 13.8834

# --- Everything at once ---
# Step 5: Summarise everything at once
summary(marks)

# NOTE FOR THE EXAM: R's var() and sd() are the n-1 versions. If a question
# asks for the population variance you must convert, as above. This is the
# same n vs n-1 distinction as Course 4 Unit 1.

In Python, python/01_descriptive.py:

"""Experiment 1 (Python equivalent) -- mean, median, mode, variance, SD.

R version: ../01_descriptive.R
The R script quotes these numbers; this file is what verifies them.
"""

import statistics
from collections import Counter

from _shared import MARKS


# Step 1: Compute the centre and the spread
def describe(values):
    n = len(values)
    mean = sum(values) / n
    ordered = sorted(values)
    median = (statistics.median(values))
    counts = Counter(values)
    top = max(counts.values())
    modes = sorted(v for v, c in counts.items() if c == top)

    # R's var() and sd() are the SAMPLE versions -- divide by n-1.
    sample_var = sum((x - mean) ** 2 for x in values) / (n - 1)
    pop_var = sum((x - mean) ** 2 for x in values) / n

    return {
        "n": n, "mean": mean, "median": median,
        "mode": modes if top > 1 else None,
        "sample_var": sample_var, "sample_sd": sample_var ** 0.5,
        "pop_var": pop_var, "pop_sd": pop_var ** 0.5,
        "range": max(values) - min(values),
    }


if __name__ == "__main__":
    # Step 2: Print them, against the R function for each
    r = describe(MARKS)
    print("EXPERIMENT 1 -- Descriptive statistics")
    print(f"  n            = {r['n']}")
    print(f"  mean         = {r['mean']:.4f}")
    print(f"  median       = {r['median']:.4f}")
    print(f"  mode         = {r['mode'] if r['mode'] else 'none (all values unique)'}")
    print(f"  range        = {r['range']}")
    print(f"  var  (n-1)   = {r['sample_var']:.4f}   <- R's var()")
    print(f"  sd   (n-1)   = {r['sample_sd']:.4f}   <- R's sd()")
    print(f"  var  (n)     = {r['pop_var']:.4f}   <- population")
    print(f"  sd   (n)     = {r['pop_sd']:.4f}")

    # Step 3: Check them against the statistics module
    assert abs(r["mean"] - statistics.mean(MARKS)) < 1e-9
    assert abs(r["sample_var"] - statistics.variance(MARKS)) < 1e-9
    assert abs(r["pop_var"] - statistics.pvariance(MARKS)) < 1e-9
    print("\n  cross-checked against the statistics module ✓")
    print("\n  NOTE: R's var() and sd() use n-1. R has no built-in mode();")
    print("        the R script defines one, as the syllabus expects.")

5. Execution and Results

In R, 01_descriptive.R:

[1] 70.05
[1] 70
 [1] 45 48 52 55 58 60 64 66 67 69 71 73 76 78 79 82 85 89 91 93
[1] 202.8921
[1] 14.24402
[1] 45 93
[1] 48
[1] 20.25
   0%   25%   50%   75%  100%
45.00 59.50 70.00 79.75 93.00
[1] 192.7475
[1] 13.88335
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  45.00   59.50   70.00   70.05   79.75   93.00

In Python, python/01_descriptive.py:

EXPERIMENT 1 -- Descriptive statistics
  n            = 20
  mean         = 70.0500
  median       = 70.0000
  mode         = none (all values unique)
  range        = 48
  var  (n-1)   = 202.8921   <- R's var()
  sd   (n-1)   = 14.2440   <- R's sd()
  var  (n)     = 192.7475   <- population
  sd   (n)     = 13.8834

  cross-checked against the statistics module ✓

  NOTE: R's var() and sd() use n-1. R has no built-in mode();
        the R script defines one, as the syllabus expects.

Every one of the twenty marks occurs once, so the mode function rightly returns all twenty. The Python equivalent gets the same mean, 70.05, and the same sample variance, 202.8921.

Experiment 2 — Binomial, normal, Poisson

1. Question

Plot the Binomial(10, 0.3), Poisson(3) and Normal(100, 15) distributions, and find probabilities from each.

2. Aim

Draw three distributions in R, and find their probabilities with the d, p and q functions.

3. Steps

In R, 02_distributions.R:

1. Plot Binomial(10, 0.3), and find P(X = 3) and P(X <= 3).
2. Plot Poisson(3), and find P(X = 3) and P(X <= 3).
3. Plot Normal(100, 15), and find the areas within 1, 2 and 3 SD.

In Python, python/02_distributions.py:

1. Tabulate Binomial(10, 0.3).
2. Tabulate Poisson(3).
3. Find the normal's areas within 1, 2 and 3 SD.
4. Check the figures the R comments quote.

4. Programme

In R, 02_distributions.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 2: Visualise Binomial, Normal and Poisson distributions
# Python equivalent: python/02_distributions.py

# --- BINOMIAL(n = 10, p = 0.3) ---
# Step 1: Plot Binomial(10, 0.3), and find P(X = 3) and P(X <= 3)
k <- 0:10
pmf <- dbinom(k, size = 10, prob = 0.3)
barplot(pmf, names.arg = k, col = "#1e7fbf",
        main = "Binomial(10, 0.3)", xlab = "k", ylab = "P(X = k)")

dbinom(3, 10, 0.3)   # 0.266828  -- P(X = 3)
pbinom(3, 10, 0.3)   # 0.649611  -- P(X <= 3)
# mean = np = 3.0 ; variance = np(1-p) = 2.1

# --- POISSON(lambda = 3) ---
# Step 2: Plot Poisson(3), and find P(X = 3) and P(X <= 3)
k <- 0:10
barplot(dpois(k, lambda = 3), names.arg = k, col = "#059669",
        main = "Poisson(3)", xlab = "k", ylab = "P(X = k)")

dpois(3, 3)   # 0.224042
ppois(3, 3)   # 0.647232
# mean = variance = lambda = 3   <- the Poisson signature

# --- NORMAL(mu = 100, sigma = 15) ---
# Step 3: Plot Normal(100, 15), and find the areas within 1, 2 and 3 SD
x <- seq(50, 150, by = 0.5)
plot(x, dnorm(x, mean = 100, sd = 15), type = "l", lwd = 2, col = "#0f4c81",
     main = "Normal(100, 15)", ylab = "density")

pnorm(115, 100, 15)                          # 0.841345
pnorm(115, 100, 15) - pnorm(85, 100, 15)     # 0.682689  <- the 68% rule
pnorm(130, 100, 15) - pnorm(70, 100, 15)     # 0.954500  <- 95%
pnorm(145, 100, 15) - pnorm(55, 100, 15)     # 0.997300  <- 99.7%
qnorm(0.95, 100, 15)                         # 124.67 -- the 95th percentile

# NAMING CONVENTION -- worth memorising, it applies to every distribution:
#   d<name>  density / PMF        dbinom, dnorm, dpois
#   p<name>  cumulative (CDF)     pbinom, pnorm, ppois
#   q<name>  quantile (inverse)   qbinom, qnorm, qpois
#   r<name>  random generation    rbinom, rnorm, rpois

In Python, python/02_distributions.py:

"""Experiment 2 (Python equivalent) -- binomial, normal and Poisson.

R version: ../02_distributions.R
Uses the statlib module written for Course 4, so the two courses' distribution
numbers are guaranteed consistent.
"""

import pathlib
import sys

sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[2]
                       / "course-4-stats" / "python"))
import statlib as S     # noqa: E402


def binomial_table(n, p):
    return [(k, S.binomial_pmf(k, n, p), S.binomial_cdf(k, n, p))
            for k in range(n + 1)]


def poisson_table(lam, upto):
    return [(k, S.poisson_pmf(k, lam), S.poisson_cdf(k, lam))
            for k in range(upto + 1)]


if __name__ == "__main__":
    # Step 1: Tabulate Binomial(10, 0.3)
    n, p = 10, 0.3
    print(f"BINOMIAL(n={n}, p={p})           R: dbinom(k, {n}, {p})")
    print(f"  {'k':<4}{'P(X=k)':<12}{'P(X<=k)':<12}")
    for k, pmf, cdf in binomial_table(n, p):
        print(f"  {k:<4}{pmf:<12.6f}{cdf:<12.6f}")
    print(f"  mean = np = {n * p:.2f}    variance = np(1-p) = {n * p * (1 - p):.2f}")

    # Step 2: Tabulate Poisson(3)
    lam = 3
    print(f"\nPOISSON(lambda={lam})            R: dpois(k, {lam})")
    print(f"  {'k':<4}{'P(X=k)':<12}{'P(X<=k)':<12}")
    for k, pmf, cdf in poisson_table(lam, 8):
        print(f"  {k:<4}{pmf:<12.6f}{cdf:<12.6f}")
    print(f"  mean = variance = lambda = {lam}")

    # Step 3: Find the normal's areas within 1, 2 and 3 SD
    mu, sigma = 100, 15
    print(f"\nNORMAL(mu={mu}, sigma={sigma})       R: pnorm(x, {mu}, {sigma})")
    for k in (1, 2, 3):
        lo, hi = mu - k * sigma, mu + k * sigma
        prob = S.normal_cdf(hi, mu, sigma) - S.normal_cdf(lo, mu, sigma)
        print(f"  within {k} sd [{lo:6.1f},{hi:6.1f}] = {prob * 100:6.3f}%")

    # Step 4: Check the figures the R comments quote
    # The empirical rule, to four decimals -- these are the numbers the R
    # script's comments quote.
    assert abs((S.normal_cdf(115, mu, sigma) - S.normal_cdf(85, mu, sigma))
               - 0.682689) < 1e-5
    assert abs(S.binomial_pmf(3, 10, 0.3) - 0.266828) < 1e-5
    assert abs(S.poisson_pmf(3, 3) - 0.224042) < 1e-5
    print("\n  all values cross-checked against Course 4's statlib ✓")

5. Execution and Results

In R, 02_distributions.R:

[1] 0.2668279
[1] 0.6496107
[1] 0.2240418
[1] 0.6472319
[1] 0.8413447
[1] 0.6826895
[1] 0.9544997
[1] 0.9973002
[1] 124.6728

In Python, python/02_distributions.py:

BINOMIAL(n=10, p=0.3)           R: dbinom(k, 10, 0.3)
  k   P(X=k)      P(X<=k)
  0   0.028248    0.028248
  1   0.121061    0.149308
  2   0.233474    0.382783
  3   0.266828    0.649611
  4   0.200121    0.849732
  5   0.102919    0.952651
  6   0.036757    0.989408
  7   0.009002    0.998410
  8   0.001447    0.999856
  9   0.000138    0.999994
  10  0.000006    1.000000
  mean = np = 3.00    variance = np(1-p) = 2.10

POISSON(lambda=3)            R: dpois(k, 3)
  k   P(X=k)      P(X<=k)
  0   0.049787    0.049787
  1   0.149361    0.199148
  2   0.224042    0.423190
  3   0.224042    0.647232
  4   0.168031    0.815263
  5   0.100819    0.916082
  6   0.050409    0.966491
  7   0.021604    0.988095
  8   0.008102    0.996197
  mean = variance = lambda = 3

NORMAL(mu=100, sigma=15)       R: pnorm(x, 100, 15)
  within 1 sd [  85.0, 115.0] = 68.269%
  within 2 sd [  70.0, 130.0] = 95.450%
  within 3 sd [  55.0, 145.0] = 99.730%

  all values cross-checked against Course 4's statlib ✓

The three plots are the binomial's and the Poisson's bar charts and the normal curve. The areas within 1, 2 and 3 SD are 0.6827, 0.9545 and 0.9973: the 68–95–99.7 rule.

Experiment 3 — t-test and chi-square

1. Question

Test whether two groups' mean scores differ, by the t-test, and whether region and purchase type are associated, by the chi-square test.

2. Aim

Run the t-test in its variants and the chi-square test of independence in R, and read their output.

3. Steps

In R, 03_hypothesis_tests.R:

1. Enter the two groups.
2. Run the pooled two-sample t-test.
3. Run Welch's test, and check the equal-variance assumption.
4. Run the one-sample and paired t-tests.
5. Test region against purchase type by chi-square.
6. Check the expected counts and the residuals.

In Python, python/03_hypothesis_tests.py:

1. Run the pooled two-sample t-test.
2. Run the chi-square test.
3. Check against Course 4 Unit 5.

4. Programme

In R, 03_hypothesis_tests.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 3: t-test and Chi-Square test
# Python equivalent: python/03_hypothesis_tests.py
# Same data as Course 4 Unit 5, so the numbers must agree with those notes.

# Step 1: Enter the two groups
group_a <- c(78, 82, 75, 88, 79, 84, 80, 86, 77, 83)
group_b <- c(72, 75, 70, 78, 74, 71, 76, 73, 69, 77)

# --- TWO-SAMPLE t-TEST ---
# Step 2: Run the pooled two-sample t-test
t.test(group_a, group_b, var.equal = TRUE)
#   t = 4.7541, df = 18, p-value = 0.000159
#   -> reject H0: the two group means differ significantly
#   mean of A = 81.20, mean of B = 73.50, pooled variance = 13.1167

# Step 3: Run Welch's test, and check the equal-variance assumption
# var.equal = TRUE gives the POOLED t-test (Course 4 Unit 5).
# Omit it and R runs Welch's t-test instead, which does not assume equal
# variances and reports fractional degrees of freedom. Both are defensible;
# know which one you asked for.
t.test(group_a, group_b)          # Welch -- note the different df

# Check the equal-variance assumption first:
var.test(group_a, group_b)        # F = 1.8618, p = 0.3682 -> variances OK

# Step 4: Run the one-sample and paired t-tests
# One-sample and paired variants:
t.test(group_a, mu = 80)
t.test(group_a, group_b, paired = TRUE)

# --- CHI-SQUARE TEST OF INDEPENDENCE ---
# Step 5: Test region against purchase type by chi-square
observed <- matrix(c(30, 70, 45, 55, 25, 75), nrow = 3, byrow = TRUE,
                   dimnames = list(c("North", "South", "East"),
                                   c("Premium", "Standard")))
observed
chisq.test(observed)
#   X-squared = 9.75, df = 2, p-value = 0.007635
#   -> reject H0: region and purchase type ARE associated

# Step 6: Check the expected counts and the residuals
result <- chisq.test(observed)
result$expected     # all 33.33 / 66.67 -- every one >= 5, assumption satisfied
result$residuals    # which cells contribute most to the statistic

# NOTE byrow = TRUE. R's matrix() fills COLUMN-wise by default, so omitting it
# silently transposes your table and changes the answer.

In Python, python/03_hypothesis_tests.py:

"""Experiment 3 (Python equivalent) -- t-test and chi-square test.

R version: ../03_hypothesis_tests.R  (t.test() and chisq.test())
"""
import pathlib, sys
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[2] / "course-4-stats" / "python"))
import statlib as S     # noqa: E402

GROUP_A = [78, 82, 75, 88, 79, 84, 80, 86, 77, 83]
GROUP_B = [72, 75, 70, 78, 74, 71, 76, 73, 69, 77]

OBSERVED = [[30, 70], [45, 55], [25, 75]]      # region x purchase type


def describe(v):
    m = sum(v) / len(v)
    return m, sum((x - m) ** 2 for x in v) / (len(v) - 1)


def two_sample_t(a, b):
    ma, va = describe(a)
    mb, vb = describe(b)
    na, nb = len(a), len(b)
    pooled = ((na - 1) * va + (nb - 1) * vb) / (na + nb - 2)
    se = (pooled * (1 / na + 1 / nb)) ** 0.5
    t = (ma - mb) / se
    df = na + nb - 2
    return t, df, S.t_sf_two_tailed(t, df), pooled, se


def chi_square(observed):
    rows = [sum(r) for r in observed]
    cols = [sum(observed[i][j] for i in range(len(observed)))
            for j in range(len(observed[0]))]
    total = sum(rows)
    chi2 = 0.0
    expected = []
    for i, row in enumerate(observed):
        exp_row = []
        for j, o in enumerate(row):
            e = rows[i] * cols[j] / total
            exp_row.append(e)
            chi2 += (o - e) ** 2 / e
        expected.append(exp_row)
    df = (len(observed) - 1) * (len(observed[0]) - 1)
    return chi2, df, S.chi2_sf(chi2, df), expected


if __name__ == "__main__":
    # Step 1: Run the pooled two-sample t-test
    t, df, p, pooled, se = two_sample_t(GROUP_A, GROUP_B)
    print("TWO-SAMPLE t-TEST        R: t.test(a, b, var.equal = TRUE)")
    print(f"  mean A   = {sum(GROUP_A)/len(GROUP_A):.4f}")
    print(f"  mean B   = {sum(GROUP_B)/len(GROUP_B):.4f}")
    print(f"  pooled var = {pooled:.4f}    se = {se:.4f}")
    print(f"  t = {t:.4f}   df = {df}   p = {p:.6f}")
    print(f"  -> {'reject' if p < 0.05 else 'fail to reject'} H0 at alpha = 0.05")

    # Step 2: Run the chi-square test
    chi2, cdf, cp, expected = chi_square(OBSERVED)
    print("\nCHI-SQUARE TEST          R: chisq.test(matrix)")
    print("  expected frequencies:")
    for row in expected:
        print("   ", "  ".join(f"{e:7.3f}" for e in row))
    print(f"  chi-square = {chi2:.4f}   df = {cdf}   p = {cp:.6f}")
    print(f"  -> {'reject' if cp < 0.05 else 'fail to reject'} H0")
    print(f"  smallest expected = {min(min(r) for r in expected):.2f} (must be >= 5)")

    # Step 3: Check against Course 4 Unit 5
    # These must agree with Course 4 Unit 5, which uses the same data.
    assert abs(t - 4.754053) < 1e-5, t
    assert abs(chi2 - 9.75) < 1e-9
    print("\n  matches Course 4 Unit 5's worked examples ✓")

5. Execution and Results

In R, 03_hypothesis_tests.R:

    Two Sample t-test

data:  group_a and group_b
t = 4.7541, df = 18, p-value = 0.0001585
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  4.297198 11.102802
sample estimates:
mean of x mean of y
     81.2      73.5


    Welch Two Sample t-test

data:  group_a and group_b
t = 4.7541, df = 16.503, p-value = 0.0001987
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  4.274943 11.125057
sample estimates:
mean of x mean of y
     81.2      73.5


    F test to compare two variances

data:  group_a and group_b
F = 1.8618, num df = 9, denom df = 9, p-value = 0.3682
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.4624493 7.4956691
sample estimates:
ratio of variances
          1.861818


    One Sample t-test

data:  group_a
t = 0.91856, df = 9, p-value = 0.3823
alternative hypothesis: true mean is not equal to 80
95 percent confidence interval:
 78.24473 84.15527
sample estimates:
mean of x
     81.2


    Paired t-test

data:  group_a and group_b
t = 7.4516, df = 9, p-value = 3.886e-05
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
  5.362438 10.037562
sample estimates:
mean difference
            7.7

      Premium Standard
North      30       70
South      45       55
East       25       75

    Pearson's Chi-squared test

data:  observed
X-squared = 9.75, df = 2, p-value = 0.007635

       Premium Standard
North 33.33333 66.66667
South 33.33333 66.66667
East  33.33333 66.66667
         Premium   Standard
North -0.5773503  0.4082483
South  2.0207259 -1.4288690
East  -1.4433757  1.0206207

In Python, python/03_hypothesis_tests.py:

TWO-SAMPLE t-TEST        R: t.test(a, b, var.equal = TRUE)
  mean A   = 81.2000
  mean B   = 73.5000
  pooled var = 13.1167    se = 1.6197
  t = 4.7541   df = 18   p = 0.000159
  -> reject H0 at alpha = 0.05

CHI-SQUARE TEST          R: chisq.test(matrix)
  expected frequencies:
     33.333   66.667
     33.333   66.667
     33.333   66.667
  chi-square = 9.7500   df = 2   p = 0.007635
  -> reject H0
  smallest expected = 33.33 (must be >= 5)

  matches Course 4 Unit 5's worked examples ✓

The pooled test has 18 degrees of freedom and Welch's 16.503: the same t, 4.7541, read against a slightly different distribution. The chi-square test's expected counts are all 33.33 or 66.67, so none is below 5.

Experiment 4 — Correlation and regression

1. Question

Find the correlation between study hours and exam scores, fit a regression line, and predict the score for 7.5 hours.

2. Aim

Fit and read a simple linear regression in R with lm().

3. Steps

In R, 04_regression.R:

1. Enter the hours and scores.
2. Measure the correlation, and plot the points.
3. Fit the regression line.
4. Use the model: coefficients, intervals, a prediction, residuals.
5. Read the ANOVA table.

In Python, python/04_regression.py:

1. Compute r, the line and the ANOVA from the sums.
2. Print the results.
3. Check that R-squared = r^2 and F = t^2.

4. Programme

In R, 04_regression.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 4: Correlation and simple linear regression
# Python equivalent: python/04_regression.py
# Same data as Course 4 Unit 4 -- coefficients must match those notes.

# Step 1: Enter the hours and scores
hours  <- c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
scores <- c(52, 55, 61, 64, 70, 72, 78, 82, 85, 91)
df <- data.frame(hours, scores)

# Step 2: Measure the correlation, and plot the points
cor(hours, scores)                       # 0.997904 -- Pearson
cor(hours, scores, method = "spearman")  # rank correlation
cov(hours, scores)

plot(hours, scores, pch = 19, col = "#1e7fbf",
     main = "Exam score against study hours")

# Step 3: Fit the regression line
model <- lm(scores ~ hours, data = df)
summary(model)
#   (Intercept)  43.0303   Std.Error 0.7011   t 61.38
#   hours         4.3030   Std.Error 0.0987   t 43.615   p 8.43e-11
#   [Corrected: this gave the intercept's standard error as 1.0847 and its
#   t as 39.671. R prints 0.70111 and 61.38, and by hand
#   SE = 0.8961 x sqrt(1/10 + 6.5^2/82.5) = 0.7011.]
#   Multiple R-squared: 0.995812      F: 1902.26 on 1 and 8 DF
#
#   Fitted line:  scores = 43.0303 + 4.3030 * hours
#   Each extra hour of study is ASSOCIATED WITH about 4.3 more marks.

abline(model, col = "red", lwd = 2)

# Step 4: Use the model: coefficients, intervals, a prediction, residuals
coef(model)
confint(model)
predict(model, newdata = data.frame(hours = 7.5))    # 75.30
residuals(model)
par(mfrow = c(2, 2)); plot(model); par(mfrow = c(1, 1))   # diagnostics

# Step 5: Read the ANOVA table
anova(model)     # SS_reg = 1527.58, SS_res = 6.42, F = 1902.26

# TWO FREE ARITHMETIC CHECKS for simple regression (Course 4 Unit 4):
#   R-squared == r^2        0.997904^2 = 0.995812  ✓
#   F         == t^2        43.615^2   = 1902.26   ✓

In Python, python/04_regression.py:

"""Experiment 4 (Python equivalent) -- correlation and simple linear regression.

R version: ../04_regression.R  (cor() and lm())
Uses the same hours/scores pair as Course 4 Unit 4, so the coefficients here
must reproduce that unit's worked example exactly.
"""
import pathlib, sys
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[2] / "course-4-stats" / "python"))
import statlib as S     # noqa: E402
from _shared import HOURS, SCORES


# Step 1: Compute r, the line and the ANOVA from the sums
def regress(x, y):
    n = len(x)
    mx, my = sum(x) / n, sum(y) / n
    sxy = sum((a - mx) * (b - my) for a, b in zip(x, y))
    sxx = sum((a - mx) ** 2 for a in x)
    syy = sum((b - my) ** 2 for b in y)

    r = sxy / (sxx * syy) ** 0.5
    b1 = sxy / sxx
    b0 = my - b1 * mx

    ss_res = sum((b - (b0 + b1 * a)) ** 2 for a, b in zip(x, y))
    ss_reg = syy - ss_res
    r2 = ss_reg / syy
    ms_res = ss_res / (n - 2)
    se_b1 = (ms_res / sxx) ** 0.5
    t = b1 / se_b1
    f = ss_reg / ms_res
    return dict(n=n, r=r, b0=b0, b1=b1, r2=r2, ss_tot=syy, ss_res=ss_res,
                ss_reg=ss_reg, ms_res=ms_res, se_b1=se_b1, t=t, f=f,
                p=S.t_sf_two_tailed(t, n - 2))


if __name__ == "__main__":
    # Step 2: Print the results
    m = regress(HOURS, SCORES)
    print("CORRELATION AND REGRESSION      R: cor(x,y) ; lm(y ~ x)")
    print(f"  Pearson r      = {m['r']:.6f}")
    print(f"  intercept b0   = {m['b0']:.4f}")
    print(f"  slope     b1   = {m['b1']:.4f}")
    print(f"  fitted line: y = {m['b0']:.4f} + {m['b1']:.4f} x")
    print(f"\n  R-squared      = {m['r2']:.6f}   (= r^2 = {m['r']**2:.6f})")
    print(f"  SE(b1)         = {m['se_b1']:.4f}")
    print(f"  t              = {m['t']:.4f}  on {m['n']-2} df,  p = {m['p']:.3e}")
    print(f"  F              = {m['f']:.4f}   (= t^2 = {m['t']**2:.4f})")
    print("\n  ANOVA")
    print(f"    Regression  SS={m['ss_reg']:10.4f}  df=1")
    print(f"    Residual    SS={m['ss_res']:10.4f}  df={m['n']-2}  MS={m['ms_res']:.4f}")
    print(f"    Total       SS={m['ss_tot']:10.4f}  df={m['n']-1}")
    print(f"\n  predict at x=7.5 -> {m['b0'] + m['b1']*7.5:.2f}")

    # Step 3: Check that R-squared = r^2 and F = t^2
    assert abs(m["r2"] - m["r"] ** 2) < 1e-12,  "R2 must equal r squared"
    assert abs(m["f"] - m["t"] ** 2) < 1e-6,    "F must equal t squared"
    assert abs(m["b1"] - 4.3030) < 1e-3
    print("\n  R2 = r^2 and F = t^2 both hold; matches Course 4 Unit 4 ✓")

5. Execution and Results

In R, 04_regression.R:

[1] 0.9979039
[1] 1
[1] 39.44444

Call:
lm(formula = scores ~ hours, data = df)

Residuals:
    Min      1Q  Median      3Q     Max
-1.1515 -0.8409  0.3030  0.6136  1.1515

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 43.03030    0.70111   61.38 5.52e-12 ***
hours        4.30303    0.09866   43.62 8.43e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8961 on 8 degrees of freedom
Multiple R-squared:  0.9958,    Adjusted R-squared:  0.9953
F-statistic:  1902 on 1 and 8 DF,  p-value: 8.425e-11

(Intercept)       hours
   43.03030     4.30303
                2.5 %   97.5 %
(Intercept) 41.413546 44.64706
hours        4.075521  4.53054
       1
75.30303
         1          2          3          4          5          6          7
 0.3636364 -0.9393939  0.7575758 -0.5454545  1.1515152 -1.1515152  0.5454545
         8          9         10
 0.2424242 -1.0606061  0.6363636
Analysis of Variance Table

Response: scores
          Df  Sum Sq Mean Sq F value    Pr(>F)
hours      1 1527.58  1527.6  1902.3 8.425e-11 ***
Residuals  8    6.42     0.8
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

In Python, python/04_regression.py:

CORRELATION AND REGRESSION      R: cor(x,y) ; lm(y ~ x)
  Pearson r      = 0.997904
  intercept b0   = 43.0303
  slope     b1   = 4.3030
  fitted line: y = 43.0303 + 4.3030 x

  R-squared      = 0.995812   (= r^2 = 0.995812)
  SE(b1)         = 0.0987
  t              = 43.6150  on 8 df,  p = 8.425e-11
  F              = 1902.2642   (= t^2 = 1902.2642)

  ANOVA
    Regression  SS= 1527.5758  df=1
    Residual    SS=    6.4242  df=8  MS=0.8030
    Total       SS= 1534.0000  df=9

  predict at x=7.5 -> 75.30

  R2 = r^2 and F = t^2 both hold; matches Course 4 Unit 4 ✓

Corrected: the comment under summary(model) gave the intercept's standard error as 1.0847, and its t as 39.671. R prints 0.70111 and 61.38, and by hand SE = 0.8961 × √(1/10 + 6.5²/82.5) = 0.7011. The slope's figures were right. The second plot is the four diagnostic plots, from plot(model).

Experiment 5 — EDA on a real dataset

1. Question

Explore a dataset: its structure, summary statistics, missing values, categories, distributions, outliers and correlations.

2. Aim

Carry out the steps of exploratory data analysis in R, on the iris data.

3. Steps

In R, 05_eda.R:

1. Load the data.
2. Look at its structure and summary.
3. Check for missing values.
4. Count the categories.
5. Plot the distributions, and find the outliers.
6. Find the correlations.
7. Judge the skew from the mean and median.

In Python, python/05_eda.py:

1. Show the structure.
2. Summarise the numeric columns.
3. Count the missing values and the categories.
4. Draw a text histogram, and find the outliers.
5. Find the correlation.

4. Programme

In R, 05_eda.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 5: Exploratory Data Analysis
# Python equivalent: python/05_eda.py
# Step 1: Load the data
data(iris)                       # or read.csv("yourfile.csv")

# Step 2: Look at its structure and summary
str(iris)                        # structure: 150 obs. of 5 variables
dim(iris); nrow(iris); ncol(iris)
head(iris); tail(iris)
summary(iris)                    # min, Q1, median, mean, Q3, max per column

# Step 3: Check for missing values
colSums(is.na(iris))             # missing values per column
sum(complete.cases(iris))        # rows with no NA at all

# Step 4: Count the categories
table(iris$Species)              # categorical counts
prop.table(table(iris$Species))  # as proportions

# Step 5: Plot the distributions, and find the outliers
hist(iris$Sepal.Length, breaks = 10, col = "#1e7fbf",
     main = "Distribution of sepal length")
boxplot(Sepal.Length ~ Species, data = iris, col = "#059669")
boxplot(iris$Sepal.Length)$out    # the outlier values themselves
pairs(iris[, 1:4], col = iris$Species)   # scatterplot matrix

# Step 6: Find the correlations
cor(iris[, 1:4])                 # correlation matrix -- numeric columns only

# Step 7: Judge the skew from the mean and median
# Skewness needs a package; the sign is what matters
# library(e1071); skewness(iris$Sepal.Length)
# mean > median  -> right-skewed ; mean < median -> left-skewed
mean(iris$Sepal.Length); median(iris$Sepal.Length)

In Python, python/05_eda.py:

"""Experiment 5 (Python equivalent) -- exploratory data analysis.

R version: ../05_eda.R  (str, summary, colSums(is.na()), hist, boxplot)
"""
from _shared import STUDENTS

COLS = ["name", "section", "gender", "hours", "marks", "attendance"]


def to_columns(rows):
    return {c: [r[i] for r in rows] for i, c in enumerate(COLS)}


def summarise_numeric(v):
    s = sorted(v)
    n = len(s)
    def q(p):
        pos = (n - 1) * p
        lo = int(pos)
        hi = min(lo + 1, n - 1)
        return s[lo] + (pos - lo) * (s[hi] - s[lo])
    return dict(min=s[0], q1=q(.25), median=q(.5), mean=sum(s)/n,
                q3=q(.75), max=s[-1])


if __name__ == "__main__":
    d = to_columns(STUDENTS)
    # Step 1: Show the structure
    print("STRUCTURE                      R: str(df)")
    print(f"  {len(STUDENTS)} observations of {len(COLS)} variables")
    for c in COLS:
        kind = "num" if isinstance(d[c][0], (int, float)) else "chr"
        print(f"    {c:<12} {kind}   e.g. {d[c][0]}")

    # Step 2: Summarise the numeric columns
    print("\nSUMMARY                        R: summary(df)")
    print(f"  {'':<12}{'Min':>8}{'1stQu':>9}{'Median':>9}{'Mean':>9}{'3rdQu':>9}{'Max':>8}")
    for c in ("hours", "marks", "attendance"):
        s = summarise_numeric(d[c])
        print(f"  {c:<12}{s['min']:>8.2f}{s['q1']:>9.2f}{s['median']:>9.2f}"
              f"{s['mean']:>9.2f}{s['q3']:>9.2f}{s['max']:>8.2f}")

    # Step 3: Count the missing values and the categories
    print("\nMISSING VALUES                 R: colSums(is.na(df))")
    for c in COLS:
        print(f"    {c:<12} {sum(1 for x in d[c] if x is None)}")

    print("\nCATEGORICAL COUNTS             R: table(df$section)")
    for c in ("section", "gender"):
        counts = {}
        for v in d[c]:
            counts[v] = counts.get(v, 0) + 1
        print(f"    {c}: " + ", ".join(f"{k}={v}" for k, v in sorted(counts.items())))

    # Step 4: Draw a text histogram, and find the outliers
    print("\nHISTOGRAM of marks             R: hist(df$marks)")
    for lo in range(40, 100, 10):
        n = sum(1 for m in d["marks"] if lo <= m < lo + 10)
        print(f"    {lo}-{lo+9}  {'#' * n * 3} {n}")

    print("\nOUTLIERS (1.5 x IQR)           R: boxplot(df$marks)$out")
    s = summarise_numeric(d["marks"])
    iqr = s["q3"] - s["q1"]
    lo, hi = s["q1"] - 1.5 * iqr, s["q3"] + 1.5 * iqr
    out = [m for m in d["marks"] if m < lo or m > hi]
    print(f"    IQR = {iqr:.2f}, fences = [{lo:.2f}, {hi:.2f}]")
    print(f"    outliers: {out if out else 'none'}")

    # Step 5: Find the correlation
    print("\nCORRELATION                    R: cor(df[, c('hours','marks')])")
    x, y = d["hours"], d["marks"]
    n = len(x); mx, my = sum(x)/n, sum(y)/n
    r = (sum((a-mx)*(b-my) for a, b in zip(x, y))
         / (sum((a-mx)**2 for a in x) * sum((b-my)**2 for b in y)) ** .5)
    print(f"    r(hours, marks) = {r:.4f}")

5. Execution and Results

In R, 05_eda.R:

'data.frame':   150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
[1] 150   5
[1] 150
[1] 5
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
145          6.7         3.3          5.7         2.5 virginica
146          6.7         3.0          5.2         2.3 virginica
147          6.3         2.5          5.0         1.9 virginica
148          6.5         3.0          5.2         2.0 virginica
149          6.2         3.4          5.4         2.3 virginica
150          5.9         3.0          5.1         1.8 virginica
  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300
 Median :5.800   Median :3.000   Median :4.350   Median :1.300
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500
       Species
 setosa    :50
 versicolor:50
 virginica :50



Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species
           0            0            0            0            0
[1] 150

    setosa versicolor  virginica
        50         50         50

    setosa versicolor  virginica
 0.3333333  0.3333333  0.3333333
numeric(0)
             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000
[1] 5.843333
[1] 5.8

In Python, python/05_eda.py:

STRUCTURE                      R: str(df)
  10 observations of 6 variables
    name         chr   e.g. Ananya
    section      chr   e.g. A
    gender       chr   e.g. F
    hours        num   e.g. 9
    marks        num   e.g. 85
    attendance   num   e.g. 92

SUMMARY                        R: summary(df)
                   Min    1stQu   Median     Mean    3rdQu     Max
  hours           2.00     4.25     6.50     6.50     8.75   11.00
  marks          41.00    56.75    71.00    69.10    83.50   91.00
  attendance     60.00    72.00    82.50    80.30    89.50   95.00

MISSING VALUES                 R: colSums(is.na(df))
    name         0
    section      0
    gender       0
    hours        0
    marks        0
    attendance   0

CATEGORICAL COUNTS             R: table(df$section)
    section: A=4, B=3, C=3
    gender: F=6, M=4

HISTOGRAM of marks             R: hist(df$marks)
    40-49  ###### 2
    50-59  ### 1
    60-69  ###### 2
    70-79  ###### 2
    80-89  ###### 2
    90-99  ### 1

OUTLIERS (1.5 x IQR)           R: boxplot(df$marks)$out
    IQR = 26.75, fences = [16.62, 123.62]
    outliers: none

CORRELATION                    R: cor(df[, c('hours','marks')])
    r(hours, marks) = 0.9932

The R script uses R's built-in iris data; the Python equivalent, which has no iris without a package, explores the ten students the other experiments use. The four plots are the histogram, the boxplot by species, the boxplot of all sepal lengths, and the scatterplot matrix.

Experiment 6 — Feature engineering

1. Question

Normalise, standardise, encode and bin a column of marks and a column of sections.

2. Aim

Prepare features for modelling in R: scale them, encode categories, and bin a number.

3. Steps

In R, 06_feature_engineering.R:

1. Enter the marks and sections.
2. Normalise to [0, 1] by min-max.
3. Standardise to mean 0 and SD 1.
4. One-hot encode the sections.
5. Encode an ordered category, and avoid the factor trap.
6. Bin the marks into classes.

In Python, python/06_feature_engineering.py:

1. Normalise by min-max.
2. Standardise.
3. One-hot encode.
4. Bin the marks.
5. Check the ranges.

4. Programme

In R, 06_feature_engineering.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 6: Scaling, normalisation and encoding
# Python equivalent: python/06_feature_engineering.py

# Step 1: Enter the marks and sections
marks <- c(85, 62, 91, 55, 74, 79, 48, 88, 68, 41)
section <- factor(c("A","A","B","B","A","C","C","B","A","C"))

# --- MIN-MAX NORMALISATION: x' = (x - min)/(max - min) -> [0, 1] ---
# Step 2: Normalise to [0, 1] by min-max
normalise <- function(x) (x - min(x)) / (max(x) - min(x))
normalise(marks)
range(normalise(marks))    # must be exactly 0 and 1

# --- STANDARDISATION: x' = (x - mean)/sd -> mean 0, sd 1 ---
# Step 3: Standardise to mean 0 and SD 1
scale(marks)               # returns a matrix; use as.vector() for a vector
as.vector(scale(marks))
mean(scale(marks)); sd(scale(marks))    # ~0 and exactly 1

# NOTE: scale() uses sd(), which divides by n-1. Python's
# sklearn StandardScaler divides by n. The values differ slightly --
# harmless for modelling, but do not expect identical numbers.

# --- ONE-HOT ENCODING ---
# Step 4: One-hot encode the sections
model.matrix(~ section - 1)     # the -1 drops the intercept -> k columns
model.matrix(~ section)         # keeps intercept -> k-1 columns (reference = A)

# For MODELLING use the k-1 form: k columns are perfectly collinear
# (they sum to 1), which is the dummy variable trap. lm() and glm() handle
# factors automatically, so you rarely encode by hand.

# --- LABEL / ORDINAL ENCODING (only for genuinely ordered categories) ---
# Step 5: Encode an ordered category, and avoid the factor trap
sizes <- factor(c("small","large","medium"),
                levels = c("small","medium","large"), ordered = TRUE)
as.numeric(sizes)          # 1 3 2 -- correct here, because the order is real

# THE CLASSIC BUG: as.numeric() on an unordered factor returns LEVEL CODES,
# not the values. For a factor of numbers use:
f <- factor(c("10","20","30"))
as.numeric(f)                    # 1 2 3   <- WRONG
as.numeric(as.character(f))      # 10 20 30 <- correct

# --- BINNING ---
# Step 6: Bin the marks into classes
cut(marks, breaks = c(0, 40, 50, 60, 75, 100),
    labels = c("Fail","Pass","Second","First","Distinction"),
    right = FALSE)
table(cut(marks, breaks = c(0, 40, 50, 60, 75, 100),
          labels = c("Fail","Pass","Second","First","Distinction"),
          right = FALSE))

In Python, python/06_feature_engineering.py:

"""Experiment 6 (Python equivalent) -- scaling, normalisation, encoding.

R version: ../06_feature_engineering.R  (scale(), model.matrix(), cut())
"""
from _shared import STUDENTS


def min_max(v):
    lo, hi = min(v), max(v)
    return [(x - lo) / (hi - lo) for x in v]


def standardise(v, sample=True):
    n = len(v)
    m = sum(v) / n
    div = (n - 1) if sample else n
    sd = (sum((x - m) ** 2 for x in v) / div) ** 0.5
    return [(x - m) / sd for x in v], m, sd


def one_hot(values):
    levels = sorted(set(values))
    return levels, [[1 if v == lv else 0 for lv in levels] for v in values]


def bin_values(v, edges, labels):
    out = []
    for x in v:
        for i in range(len(edges) - 1):
            if edges[i] <= x < edges[i + 1]:
                out.append(labels[i]); break
        else:
            out.append(labels[-1])
    return out


if __name__ == "__main__":
    marks = [r[4] for r in STUDENTS]
    sections = [r[1] for r in STUDENTS]

    # Step 1: Normalise by min-max
    print("MIN-MAX NORMALISATION          x' = (x - min)/(max - min)")
    nm = min_max(marks)
    for m, v in list(zip(marks, nm))[:5]:
        print(f"    {m:>3} -> {v:.4f}")
    print(f"    range check: min={min(nm):.4f}  max={max(nm):.4f}   (must be 0 and 1)")

    # Step 2: Standardise
    print("\nSTANDARDISATION                R: scale(x)  -- uses n-1")
    st, mean, sd = standardise(marks)
    print(f"    mean = {mean:.4f}   sd (n-1) = {sd:.4f}")
    for m, v in list(zip(marks, st))[:5]:
        print(f"    {m:>3} -> {v:+.4f}")
    chk_m = sum(st) / len(st)
    chk_s = (sum((x - chk_m) ** 2 for x in st) / (len(st) - 1)) ** 0.5
    print(f"    check: mean={chk_m:.10f} sd={chk_s:.6f}   (must be 0 and 1)")

    # Step 3: One-hot encode
    print("\nONE-HOT ENCODING               R: model.matrix(~ section - 1)")
    levels, encoded = one_hot(sections)
    print(f"    levels: {levels}")
    for s, e in list(zip(sections, encoded))[:5]:
        print(f"    {s} -> {e}")
    print("    NOTE: for a linear model R drops one level as the reference,")
    print("          giving k-1 columns and avoiding the dummy variable trap.")

    # Step 4: Bin the marks
    print("\nBINNING                        R: cut(marks, breaks = ...)")
    labels = ["Fail", "Pass", "Second", "First", "Distinction"]
    binned = bin_values(marks, [0, 40, 50, 60, 75, 101], labels)
    counts = {}
    for b in binned:
        counts[b] = counts.get(b, 0) + 1
    for lab in labels:
        print(f"    {lab:<12} {counts.get(lab, 0)}")

    # Step 5: Check the ranges
    assert abs(min(nm)) < 1e-12 and abs(max(nm) - 1) < 1e-12
    assert abs(chk_m) < 1e-10 and abs(chk_s - 1) < 1e-10
    print("\n  normalisation spans [0,1] and standardisation gives mean 0, sd 1 ✓")

5. Execution and Results

In R, 06_feature_engineering.R:

 [1] 0.88 0.42 1.00 0.28 0.66 0.76 0.14 0.94 0.54 0.00
[1] 0 1
             [,1]
 [1,]  0.91851437
 [2,] -0.41015422
 [3,]  1.26512357
 [4,] -0.81453162
 [5,]  0.28306418
 [6,]  0.57190518
 [7,] -1.21890901
 [8,]  1.09181897
 [9,] -0.06354502
[10,] -1.62328641
attr(,"scaled:center")
[1] 69.1
attr(,"scaled:scale")
[1] 17.31056
 [1]  0.91851437 -0.41015422  1.26512357 -0.81453162  0.28306418  0.57190518
 [7] -1.21890901  1.09181897 -0.06354502 -1.62328641
[1] 3.496986e-16
[1] 1
   sectionA sectionB sectionC
1         1        0        0
2         1        0        0
3         0        1        0
4         0        1        0
5         1        0        0
6         0        0        1
7         0        0        1
8         0        1        0
9         1        0        0
10        0        0        1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$section
[1] "contr.treatment"

   (Intercept) sectionB sectionC
1            1        0        0
2            1        0        0
3            1        1        0
4            1        1        0
5            1        0        0
6            1        0        1
7            1        0        1
8            1        1        0
9            1        0        0
10           1        0        1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$section
[1] "contr.treatment"

[1] 1 3 2
[1] 1 2 3
[1] 10 20 30
 [1] Distinction First       Distinction Second      First       Distinction
 [7] Pass        Distinction First       Pass
Levels: Fail Pass Second First Distinction

       Fail        Pass      Second       First Distinction
          0           2           1           3           4

In Python, python/06_feature_engineering.py:

MIN-MAX NORMALISATION          x' = (x - min)/(max - min)
     85 -> 0.8800
     62 -> 0.4200
     91 -> 1.0000
     55 -> 0.2800
     74 -> 0.6600
    range check: min=0.0000  max=1.0000   (must be 0 and 1)

STANDARDISATION                R: scale(x)  -- uses n-1
    mean = 69.1000   sd (n-1) = 17.3106
     85 -> +0.9185
     62 -> -0.4102
     91 -> +1.2651
     55 -> -0.8145
     74 -> +0.2831
    check: mean=0.0000000000 sd=1.000000   (must be 0 and 1)

ONE-HOT ENCODING               R: model.matrix(~ section - 1)
    levels: ['A', 'B', 'C']
    A -> [1, 0, 0]
    A -> [1, 0, 0]
    B -> [0, 1, 0]
    B -> [0, 1, 0]
    A -> [1, 0, 0]
    NOTE: for a linear model R drops one level as the reference,
          giving k-1 columns and avoiding the dummy variable trap.

BINNING                        R: cut(marks, breaks = ...)
    Fail         0
    Pass         2
    Second       1
    First        3
    Distinction  4

  normalisation spans [0,1] and standardisation gives mean 0, sd 1 ✓

Corrected: this note said the standardised values differ between the R script and its Python equivalent. The Python equivalent divides by n − 1, as R does, so they agree: 85 becomes 0.9185 in both. It is scikit-learn's StandardScaler that would differ.

Experiment 7 — Variables, control structures, functions

1. Question

Use R's variable types, vectors, control structures, apply family and functions.

2. Aim

Write R's basic constructs, and see where R differs from other languages.

3. Steps

1. Make a variable of each type.
2. Index and compute on vectors.
3. Branch with if and ifelse().
4. Loop with for, while and repeat.
5. Apply a function across a matrix or a list.
6. Write functions, with default and variadic arguments.

4. Programme

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 7: Variables, control structures and functions
# No Python equivalent -- this experiment demonstrates R SYNTAX, so a
# translation would teach nothing.

# --- VARIABLES AND TYPES ---
# Step 1: Make a variable of each type
x <- 10L          # integer (note the L)
y <- 3.14         # numeric / double
name <- "Ravi"    # character
flag <- TRUE      # logical
z <- 3 + 4i       # complex

class(x); class(y); class(name); class(flag)
# NOTE: class(10) is "numeric", NOT "integer". Write 10L for an integer.

# --- VECTORS: R's fundamental unit ---
# Step 2: Index and compute on vectors
v <- c(10, 20, 30, 40, 50)
v[1]        # 10  -- R indexes from 1, not 0
v[-1]       # 20 30 40 50  -- negative EXCLUDES; it is not "from the end"
v[v > 25]   # 30 40 50
length(v)

v * 2       # 20 40 60 80 100  -- VECTORISED, no loop needed
c(1,2,3,4) + c(10,20)   # 11 22 13 24 -- RECYCLING of the shorter vector

# --- CONTROL STRUCTURES ---
# Step 3: Branch with if and ifelse()
marks <- 72
if (marks >= 40) {
  print("Pass")
} else if (marks >= 30) {
  print("Supplementary")
} else {
  print("Fail")
}

# ifelse() is VECTORISED -- use it on a whole column, never if()
ifelse(v > 25, "high", "low")

# Step 4: Loop with for, while and repeat
for (i in 1:5) print(i)
for (nm in c("A", "B", "C")) print(nm)

i <- 1
while (i <= 5) { print(i); i <- i + 1 }

repeat { i <- i + 1; if (i > 10) break }

# --- THE apply FAMILY: R's idiomatic alternative to loops ---
# Step 5: Apply a function across a matrix or a list
m <- matrix(1:6, nrow = 2)
apply(m, 1, sum)     # row sums     -- MARGIN 1 = rows
apply(m, 2, mean)    # column means -- MARGIN 2 = columns
sapply(1:5, function(k) k^2)          # 1 4 9 16 25 -- returns a vector
lapply(1:3, function(k) k^2)          # returns a LIST

# --- FUNCTIONS ---
# Step 6: Write functions, with default and variadic arguments
grade <- function(marks, pass_mark = 40) {
  if (marks >= 90) return("A")
  if (marks >= 75) return("B")
  if (marks >= 60) return("C")
  if (marks >= pass_mark) return("D")
  "F"                      # the last expression is returned automatically
}
grade(85)                  # "B"
grade(35, pass_mark = 30)  # "D" -- named argument

total <- function(...) sum(...)     # variadic
total(1, 2, 3, 4)                   # 10

5. Execution and Results

[1] "integer"
[1] "numeric"
[1] "character"
[1] "logical"
[1] 10
[1] 20 30 40 50
[1] 30 40 50
[1] 5
[1]  20  40  60  80 100
[1] 11 22 13 24
[1] "Pass"
[1] "low"  "low"  "high" "high" "high"
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] "A"
[1] "B"
[1] "C"
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1]  9 12
[1] 1.5 3.5 5.5
[1]  1  4  9 16 25
[[1]]
[1] 1

[[2]]
[1] 4

[[3]]
[1] 9

[1] "B"
[1] "D"
[1] 10

Two of R's habits show in the output: indexing starts at 1, so v[1] is 10, and v[-1] leaves the first element out rather than taking the last.

Experiment 8 — CSV, Excel, JSON, XML

1. Question

Read and write data as CSV, Excel, JSON and XML files.

2. Aim

Move a data frame in and out of R in each common file format.

3. Steps

In R, 08_file_io.R:

1. Make a small data frame.
2. Write and read CSV.
3. Load the Excel packages.
4. Write and read JSON.
5. Load the XML package.
6. Save and load R's own formats.

In Python, python/08_file_io.py:

1. Write and read back CSV, JSON and XML.
2. Compare the types that come back.

4. Programme

In R, 08_file_io.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 8: Read and write CSV, Excel, JSON and XML
# Python equivalent: python/08_file_io.py

# Step 1: Make a small data frame
df <- data.frame(name = c("Ananya","Charan","Divya"),
                 section = c("A","B","B"),
                 marks = c(85, 91, 55),
                 stringsAsFactors = FALSE)

# --- CSV ---
# Step 2: Write and read CSV
write.csv(df, "students.csv", row.names = FALSE)
back <- read.csv("students.csv", stringsAsFactors = FALSE)
# row.names = FALSE matters: without it R writes an extra index column and
# re-reading gives you a stray "X" column you did not ask for.

library(readr)                     # faster; returns a tibble; never factorises
write_csv(df, "students2.csv"); read_csv("students2.csv")

# --- EXCEL ---
# Step 3: Load the Excel packages
library(readxl)
# read_excel("students.xlsx", sheet = 1)
# excel_sheets("students.xlsx")    # list the sheet names first
library(writexl)
# write_xlsx(df, "students.xlsx")

# --- JSON ---
# Step 4: Write and read JSON
library(jsonlite)
write_json(df, "students.json", pretty = TRUE)
fromJSON("students.json")          # comes back as a data frame directly
toJSON(df, pretty = TRUE, auto_unbox = TRUE)
# JSON preserves TYPES -- numbers come back as numbers. CSV and XML do not.

# --- XML ---
# Step 5: Load the XML package
library(XML)
# doc <- xmlParse("students.xml")
# xmlToDataFrame(doc)
# Alternative, often easier: library(xml2); read_xml(); xml_find_all()

# --- R's own formats ---
# Step 6: Save and load R's own formats
saveRDS(df, "students.rds"); readRDS("students.rds")   # ONE object
save(df, file = "students.RData"); load("students.RData")  # several, by name

# saveRDS/readRDS is preferred: you choose the variable name on load.
# load() silently overwrites whatever names were saved.

In Python, python/08_file_io.py:

"""Experiment 8 (Python equivalent) -- read/write CSV, JSON, XML.

R version: ../08_file_io.R  (read.csv, jsonlite, XML, readxl)
Excel is omitted here: writing .xlsx needs a third-party library, and the point
of the experiment -- that each format round-trips -- is made by the other three.
"""
import csv, json, tempfile, pathlib
import xml.etree.ElementTree as ET
from _shared import STUDENTS

COLS = ["name", "section", "gender", "hours", "marks", "attendance"]
ROWS = [dict(zip(COLS, r)) for r in STUDENTS]


def csv_roundtrip(d):
    p = d / "students.csv"
    with open(p, "w", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=COLS); w.writeheader(); w.writerows(ROWS)
    with open(p, newline="") as fh:
        back = list(csv.DictReader(fh))
    # everything read from CSV is a string -- the same trap as in Course 3
    return back, all(isinstance(v, str) for v in back[0].values())


def json_roundtrip(d):
    p = d / "students.json"
    p.write_text(json.dumps(ROWS, indent=2))
    back = json.loads(p.read_text())
    return back, isinstance(back[0]["marks"], int)      # JSON preserves types


def xml_roundtrip(d):
    p = d / "students.xml"
    root = ET.Element("students")
    for r in ROWS:
        el = ET.SubElement(root, "student")
        for k, v in r.items():
            ET.SubElement(el, k).text = str(v)
    ET.ElementTree(root).write(p, encoding="utf-8", xml_declaration=True)
    back = [{c.tag: c.text for c in el} for el in ET.parse(p).getroot()]
    return back, all(isinstance(v, str) for v in back[0].values())


if __name__ == "__main__":
    # Step 1: Write and read back CSV, JSON and XML
    with tempfile.TemporaryDirectory() as td:
        d = pathlib.Path(td)
        for name, fn in (("CSV", csv_roundtrip), ("JSON", json_roundtrip),
                         ("XML", xml_roundtrip)):
            back, typed = fn(d)
            same = len(back) == len(ROWS) and back[0]["name"] == ROWS[0]["name"]
            print(f"{name:<6} wrote {len(ROWS)} rows, read {len(back)} back  "
                  f"-> {'round-trips ✓' if same else 'MISMATCH'}")
            if name == "JSON":
                print(f"       numeric types preserved: {typed}")
            else:
                print(f"       all values come back as strings: {typed}")

    # Step 2: Compare the types that come back
    print("\n  KEY POINT: CSV and XML are untyped -- everything returns as text,")
    print("  so numbers must be converted explicitly. JSON preserves numbers and")
    print("  booleans. In R this is why read.csv() has colClasses= and why")
    print("  jsonlite::fromJSON() gives you a usable data frame directly.")

5. Execution and Results

In R, 08_file_io.R:

Rows: 3 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): name, section
dbl (1): marks

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# A tibble: 3 × 3
  name   section marks
  <chr>  <chr>   <dbl>
1 Ananya A          85
2 Charan B          91
3 Divya  B          55
    name section marks
1 Ananya       A    85
2 Charan       B    91
3  Divya       B    55
[
  {
    "name": "Ananya",
    "section": "A",
    "marks": 85
  },
  {
    "name": "Charan",
    "section": "B",
    "marks": 91
  },
  {
    "name": "Divya",
    "section": "B",
    "marks": 55
  }
]
    name section marks
1 Ananya       A    85
2 Charan       B    91
3  Divya       B    55

In Python, python/08_file_io.py:

CSV    wrote 10 rows, read 10 back  -> round-trips ✓
       all values come back as strings: True
JSON   wrote 10 rows, read 10 back  -> round-trips ✓
       numeric types preserved: True
XML    wrote 10 rows, read 10 back  -> round-trips ✓
       all values come back as strings: True

  KEY POINT: CSV and XML are untyped -- everything returns as text,
  so numbers must be converted explicitly. JSON preserves numbers and
  booleans. In R this is why read.csv() has colClasses= and why
  jsonlite::fromJSON() gives you a usable data frame directly.

The Excel and XML reads and writes are comments in the file, as the steps say: only their packages are loaded. JSON brings the marks back as numbers; CSV and XML bring everything back as text, as the Python equivalent shows.

Experiment 9 — dplyr and tidyr

1. Question

Filter, select, create, sort, summarise, join and reshape a table of students.

2. Aim

Wrangle data with dplyr's verbs and tidyr's pivots.

3. Steps

In R, 09_wrangling.R:

1. Make the students data frame.
2. Filter, select, mutate and arrange in one pipe.
3. Summarise by section.
4. Count, find distinct values, take the top three, rename.
5. Join to the teachers table.
6. Reshape from wide to long and back.

In Python, python/09_wrangling.py:

1. Filter, select, mutate and arrange.
2. Group and summarise.
3. Pivot longer and wider.
4. Check the shapes.

4. Programme

In R, 09_wrangling.R:

# =====================================================================
# Run with R 4.3.3 (Rscript --vanilla). What it prints, and the plots it
# draws, are on the lab page, and tools/data-science/run_r_equivalents.py
# runs it again. (Until October 2026 R could not be installed where these
# labs are checked, so this file was desk-checked only; every number in its
# comments has since been checked against R's own output.)
# =====================================================================
# Experiment 9: Data wrangling with dplyr and tidyr
# Python equivalent: python/09_wrangling.py (annotated with these dplyr calls)

library(dplyr); library(tidyr)

# Step 1: Make the students data frame
students <- data.frame(
  name    = c("Ananya","Bhavana","Charan","Divya","Eshwar",
              "Fiona","Gopal","Harika","Ismail","Jyothi"),
  section = c("A","A","B","B","A","C","C","B","A","C"),
  hours   = c(9, 5, 11, 4, 7, 8, 3, 10, 6, 2),
  marks   = c(85, 62, 91, 55, 74, 79, 48, 88, 68, 41),
  stringsAsFactors = FALSE
)

# --- THE FIVE VERBS, chained with the pipe ---
# Step 2: Filter, select, mutate and arrange in one pipe
students %>%
  filter(marks > 60) %>%                    # rows      -- SQL WHERE
  select(name, section, marks) %>%          # columns   -- SQL SELECT
  mutate(grade = case_when(                 # new column
    marks >= 90 ~ "A",
    marks >= 75 ~ "B",
    marks >= 60 ~ "C",
    TRUE        ~ "F")) %>%                 # TRUE ~ is the else branch
  arrange(desc(marks))                      # sort      -- SQL ORDER BY

# --- GROUPED SUMMARY -- SQL's GROUP BY ---
# Step 3: Summarise by section
students %>%
  group_by(section) %>%
  summarise(n        = n(),
            avg      = mean(marks),
            highest  = max(marks),
            pass_pct = mean(marks >= 40) * 100,   # mean of a logical!
            .groups  = "drop") %>%
  arrange(desc(avg))

# mean() of a logical vector gives a PROPORTION, because TRUE counts as 1.
# .groups = "drop" ungroups the result -- omit it and later operations
# silently stay grouped, which is a common source of confusion.

# --- USEFUL EXTRAS ---
# Step 4: Count, find distinct values, take the top three, rename
students %>% count(section)
students %>% distinct(section)
students %>% slice_max(marks, n = 3)
students %>% rename(score = marks)

# --- JOINS: the same seven as Course 5 ---
# Step 5: Join to the teachers table
sections <- data.frame(section = c("A","B","C","D"),
                       teacher = c("Rao","Devi","Kumar","Reddy"))
inner_join(students, sections, by = "section")
left_join (students, sections, by = "section")
anti_join (sections, students, by = "section")   # section D has no students

# --- RESHAPING with tidyr ---
# Step 6: Reshape from wide to long and back
wide <- data.frame(name = c("A","B"), maths = c(85,72),
                   science = c(78,88), english = c(92,65))

long <- wide %>%
  pivot_longer(cols = c(maths, science, english),
               names_to = "subject", values_to = "marks")
long           # 6 rows: one per student-subject pair

long %>% pivot_wider(names_from = subject, values_from = marks)   # back again

# ggplot2 WANTS LONG DATA. That is the practical reason this matters:
# to draw one bar per subject, subject must be a COLUMN, not three columns.

In Python, python/09_wrangling.py:

"""Experiment 9 (Python equivalent) -- dplyr/tidyr operations in pandas.

R version: ../09_wrangling.R
Every pandas call is annotated with its dplyr counterpart, which is the point
of this file: the two libraries do the same six things.
"""
import pandas as pd
from _shared import STUDENTS

COLS = ["name", "section", "gender", "hours", "marks", "attendance"]
df = pd.DataFrame(STUDENTS, columns=COLS)

if __name__ == "__main__":
    pd.set_option("display.width", 100)

    # Step 1: Filter, select, mutate and arrange
    print("filter()   dplyr: filter(df, marks > 70)")
    print(df[df.marks > 70][["name", "section", "marks"]].to_string(index=False))

    print("\nselect()   dplyr: select(df, name, marks)")
    print(df[["name", "marks"]].head(3).to_string(index=False))

    print("\nmutate()   dplyr: mutate(df, grade = case_when(...))")
    df["grade"] = pd.cut(df.marks, bins=[0, 40, 60, 75, 90, 101],
                         labels=["F", "D", "C", "B", "A"], right=False)
    print(df[["name", "marks", "grade"]].head(5).to_string(index=False))

    print("\narrange()  dplyr: arrange(df, desc(marks))")
    print(df.sort_values("marks", ascending=False)[["name", "marks"]]
            .head(3).to_string(index=False))

    # Step 2: Group and summarise
    print("\ngroup_by + summarise")
    print("  dplyr: group_by(section) %>% summarise(n=n(), avg=mean(marks))")
    g = (df.groupby("section")
           .agg(n=("marks", "size"), avg=("marks", "mean"),
                highest=("marks", "max"),
                pass_pct=("marks", lambda s: (s >= 40).mean() * 100))
           .reset_index())
    print(g.to_string(index=False))

    # Step 3: Pivot longer and wider
    print("\npivot_longer()  tidyr: pivot_longer(cols = c(hours, marks))")
    long = df.melt(id_vars=["name", "section"], value_vars=["hours", "marks"],
                   var_name="measure", value_name="value")
    print(long.head(4).to_string(index=False))
    print(f"  wide {df.shape} -> long {long.shape}")

    print("\npivot_wider()   tidyr: pivot_wider(names_from, values_from)")
    wide = long.pivot_table(index=["name", "section"], columns="measure",
                            values="value").reset_index()
    print(wide.head(3).to_string(index=False))

    # Step 4: Check the shapes
    assert len(long) == len(df) * 2, "melt must double the rows for two measures"
    assert set(g.section) == {"A", "B", "C"}
    print("\n  long form has 2x the rows, as pivot_longer would produce ✓")