MATH40006 Practice Paper B Solutions

Notes — cell 1

# MATH40006 - Introduction to Computation

## Practice Paper Paper B - Executed solutions

**REVIEW CANDIDATE - NOT PUBLISHED**

- Practice duration: 90 minutes
- Total: 75 marks; answer all three questions
- Work in this single Jupyter notebook.
- Use code cells for programs, tests, figures and requested output.
- Use Markdown answer cells for explanations, invariants, derivations and comments.
- No external files, network access, `input()` or package installation are required.
- Before submission, restart the kernel and run all cells from top to bottom.

The 90-minute/75-mark structure is a review assumption based on the nearest
alternative-assessment precedent and must be checked against Tony's own notice.

Code — cell 2

%matplotlib inline
import sys
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sympy as sp
import pandas as pd
print('Environment ready:', sys.version.split()[0], 'NumPy', np.__version__, 'Matplotlib', matplotlib.__version__)
print('SymPy', sp.__version__, 'pandas', pd.__version__)

Environment ready: 3.12.10 NumPy 2.3.5 Matplotlib 3.11.1
SymPy 1.14.0 pandas 3.0.1

Notes — cell 3

## Question 1 (30 marks): Vectorised absorbing random walks

A walk occupies an integer in `{0,1,...,L}`. From an interior point it moves
right with probability `p` and left otherwise, stopping at `0` or `L`. The
valid contract requires positive integer `n_walks,L`, integer `x0` in `[0,L]`,
`p` in `[0,1]`, and non-negative integer `max_steps`.

**(a) [3]** With `np.random.default_rng(seed)`, demonstrate a Boolean array
encoding right/left moves for several walkers.

**(b) [10]** Write
`absorbing_walks(n_walks,L,x0,p=0.5,max_steps=10000,seed=0)`. Vectorise over
walkers and loop only over time. Return Boolean `hit_right` and integer
absorption times. Raise a clear error if active walks remain at `max_steps`.

**(c) [5]** Test shapes, fixed-seed reproducibility, non-negative times,
immediate absorption at `x0=0,L`, and invalid inputs.

**(d) [5]** For 50000 fair walks with `L=20,x0=7,seed=20260829`, estimate the
right-hit probability and mean absorption time. Compare with `x0/L` and
`x0*(L-x0)`.

**(e) [4]** Construct the stated normal-approximation 95% interval. In
Markdown, check whether `7/20` lies inside and interpret cautiously. Random
agreement must not be a pass/fail assertion.

**(f) [3]** Plot an absorption-time histogram. In Markdown, comment on its
shape and explain the `max_steps` safety contract.

Notes — cell 4

### Q1 code answers for parts (a)-(f) - solution

Code — cell 5

import numpy as np
import matplotlib.pyplot as plt


demo_rng = np.random.default_rng(20260829)
demo_right_moves = demo_rng.random(8) < 0.5
assert demo_right_moves.dtype == bool and demo_right_moves.shape == (8,)


def absorbing_walks(n_walks, L, x0, p=0.5, max_steps=10000, seed=0):
    """Simulate vectorised walks absorbed at 0 or L."""
    integer_inputs = (n_walks, L, x0, max_steps)
    if any(not isinstance(value, (int, np.integer))
           or isinstance(value, (bool, np.bool_)) for value in integer_inputs):
        raise TypeError("n_walks, L, x0 and max_steps must be integers")
    if n_walks < 1 or L < 1 or not (0 <= x0 <= L):
        raise ValueError("invalid size or starting point")
    if not (0 <= p <= 1) or max_steps < 0:
        raise ValueError("invalid probability or max_steps")
    rng = np.random.default_rng(seed)
    positions = np.full(n_walks, x0, dtype=int)
    steps = np.zeros(n_walks, dtype=int)
    active = (positions > 0) & (positions < L)
    elapsed = 0
    while active.any() and elapsed < max_steps:
        move_right = rng.random(n_walks) < p
        positions[active] += np.where(move_right[active], 1, -1)
        steps[active] += 1
        active = (positions > 0) & (positions < L)
        elapsed += 1
    if active.any():
        raise RuntimeError("max_steps exceeded")
    hit_right = positions == L
    return hit_right, steps


right1, steps1 = absorbing_walks(500, 10, 3, seed=20260829)
right2, steps2 = absorbing_walks(500, 10, 3, seed=20260829)
assert np.array_equal(right1, right2)
assert np.array_equal(steps1, steps2)
assert right1.shape == steps1.shape == (500,)
assert np.all(steps1 >= 1)
assert absorbing_walks(4, 5, 0, seed=1)[1].tolist() == [0, 0, 0, 0]
right_at_L, steps_at_L = absorbing_walks(4, 5, 5, seed=1)
assert right_at_L.all() and steps_at_L.tolist() == [0, 0, 0, 0]
try:
    absorbing_walks(10, 5, 6)
except ValueError:
    pass
else:
    raise AssertionError("invalid starting point was not rejected")
try:
    absorbing_walks(10.5, 5, 2)
except TypeError:
    pass
else:
    raise AssertionError("non-integer n_walks was not rejected")

hit_right, durations = absorbing_walks(
    50000, L=20, x0=7, p=0.5, max_steps=10000, seed=20260829
)
p_hat = hit_right.mean()
mean_duration = durations.mean()
se = np.sqrt(p_hat * (1 - p_hat) / len(hit_right))
interval_95 = (p_hat - 1.96 * se, p_hat + 1.96 * se)
contains_theory = interval_95[0] <= 7/20 <= interval_95[1]

fig, ax = plt.subplots(figsize=(7.0, 4.2))
ax.hist(durations, bins=45, color="#2563A7", edgecolor="white")
ax.set(xlabel="absorption time", ylabel="number of walks",
       title="Fair absorbing walk: L=20, x0=7")
fig.tight_layout()
plt.show()

<Figure size 700x420 with 1 Axes>

Code — cell 6

print(f"B1 p_hat={p_hat:.6f}, mean={mean_duration:.5f}, 95% CI={interval_95}, contains 7/20={contains_theory}")

B1 p_hat=0.350960, mean=91.16992, 95% CI=(np.float64(0.34677654009130854), np.float64(0.35514345990869145)), contains 7/20=True

Notes — cell 7

### Q1 written answers for parts (d)-(f) - solution

For the fair walk, the theoretical benchmarks are `P(hit L)=x0/L=7/20` and `E[T]=x0(L-x0)=91`. The displayed interval is a random, approximate diagnostic rather than proof, so agreement is recorded as a Boolean result and not asserted. The histogram is right-skewed with a long tail. `max_steps` guarantees controlled failure if inputs or a run behave unexpectedly.

Notes — cell 8

## Question 2 (35 marks): Horner evaluation, derivatives and Newton iteration

`[c0,...,cn]` represents `p(x)=c0*x**n+...+cn` in descending powers.

**(a) [6]** Write `horner(coefficients,x)` using `value=value*x+c`. Accept a
real scalar or real NumPy array, reject an empty list and do not mutate input.

**(b) [5]** Test a constant polynomial, vector `x`, repeated/zero coefficients
and agreement with `np.polyval`.

**(c) [4]** In Markdown, give the exact multiplication/addition counts and
Theta time; contrast with rebuilding every power by repeated multiplication.

**(d) [8]** Derive in Markdown and implement
`horner_with_derivative(coefficients,x)`, using the old polynomial accumulator
when updating the derivative.

**(e) [6]** For coefficients `[1,1,-7,-1,6]` and 601 points on `[-3,3]`,
verify `p` with `np.polyval`, verify `p'` with SymPy `diff` and `lambdify`, and
plot both.

**(f) [6]** Write guarded Newton iteration using extended Horner. Stop on a
step below `tol`; reject a near-zero derivative and exhausted `max_iter`.
Starting at `0.8`, verify the root and residual.

Notes — cell 9

### Q2 code answers for parts (a), (b), (d), (e) and (f) - solution

Code — cell 10

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp


def horner(coefficients, x):
    """Evaluate descending-power coefficients by Horner's method."""
    coeffs = list(coefficients)
    if not coeffs:
        raise ValueError("at least one coefficient is required")
    scalar_input = np.ndim(x) == 0
    value = np.full(np.shape(x), coeffs[0], dtype=float)
    for coefficient in coeffs[1:]:
        value = value * x + coefficient
    return value.item() if scalar_input else value


def horner_with_derivative(coefficients, x):
    """Evaluate p and p' simultaneously by extended Horner."""
    coeffs = list(coefficients)
    if not coeffs:
        raise ValueError("at least one coefficient is required")
    scalar_input = np.ndim(x) == 0
    value = np.full(np.shape(x), coeffs[0], dtype=float)
    derivative = np.zeros(np.shape(x), dtype=float)
    for coefficient in coeffs[1:]:
        derivative = derivative * x + value
        value = value * x + coefficient
    if scalar_input:
        return value.item(), derivative.item()
    return value, derivative


def newton_horner(coefficients, x0, tol=1e-12, max_iter=50):
    """Newton iteration using extended Horner for p and p'."""
    x = float(x0)
    for iteration in range(max_iter):
        value, derivative = horner_with_derivative(coefficients, x)
        if abs(derivative) < 1e-14:
            raise RuntimeError("derivative too small")
        new_x = x - value / derivative
        if abs(new_x - x) < tol:
            return new_x, iteration + 1
        x = new_x
    raise RuntimeError("maximum iterations exceeded")


coeffs = [1, 1, -7, -1, 6]
saved_coeffs = coeffs.copy()
x_grid = np.linspace(-3, 3, 601)
values, derivatives = horner_with_derivative(coeffs, x_grid)
assert coeffs == saved_coeffs
assert np.allclose(values, np.polyval(coeffs, x_grid))
assert horner([7], np.array([-1.0, 2.0])).tolist() == [7.0, 7.0]
zero_case = [2, 0, 0, -3]
assert np.allclose(horner(zero_case, x_grid), np.polyval(zero_case, x_grid))

s = sp.symbols("s", real=True)
symbolic = sum(c * s**(len(coeffs)-1-i) for i, c in enumerate(coeffs))
symbolic_derivative = sp.diff(symbolic, s)
derivative_function = sp.lambdify(s, symbolic_derivative, "numpy")
assert np.allclose(derivatives, derivative_function(x_grid))

root, iterations = newton_horner(coeffs, 0.8)
assert abs(root - 1.0) < 1e-12
assert abs(horner(coeffs, root)) < 1e-12

fig, ax = plt.subplots(figsize=(7.0, 4.2))
ax.plot(x_grid, values, label="p(x)")
ax.plot(x_grid, derivatives, label="p'(x)")
ax.axhline(0, color="black", linewidth=0.7)
ax.set(xlabel="x", ylabel="value", title="Horner evaluation and derivative")
ax.legend()
fig.tight_layout()
plt.show()

<Figure size 700x420 with 1 Axes>

Code — cell 11

print(f"B2 root={root:.12f}, iterations={iterations}, residual={abs(horner(coeffs, root)):.3e}")

B2 root=1.000000000000, iterations=4, residual=0.000e+00

Notes — cell 12

### Q2 written answers for parts (c) and (d) - solution

Degree-n Horner evaluation makes exactly n multiplications and n additions, hence Theta(n); separately rebuilding each power by repeated multiplication costs Theta(n**2). Differentiating the nested update gives `derivative=derivative*x+value`, which must use the old `value` before `value=value*x+c`. Newton checks a nearly zero derivative and maximum iterations; both final step size and residual are reported.

Notes — cell 13

## Question 3 (10 marks): A pandas report for a supplied simulation log

Run the non-assessed **GIVEN SETUP** cell below unchanged. It creates a
reproducible synthetic run log, so this question has no external-file or
Question 1 dependency.

**(a) [3]** Build a DataFrame with one row per record and columns `side`
(`'left'/'right'`) and `steps`.

**(b) [4]** With `groupby` and `agg`, report count, mean, median and maximum
duration for each side. Verify counts total `report_n`.

**(c) [3]** Find the 99th percentile of `steps`, use `query` to select records
at or above it, and count by side. In Markdown, explain why this conditional
subset cannot estimate the unconditional right-side probability.

Code — cell 14

import numpy as np
import pandas as pd


# GIVEN SETUP - run this cell unchanged before Question 3.
report_rng = np.random.default_rng(4000603)
report_n = 50000
report_hit_right = report_rng.random(report_n) < 0.35
report_durations = report_rng.geometric(1 / 91.0, size=report_n)

Notes — cell 15

### Q3 code answers for parts (a)-(c) - solution

Code — cell 16

walks = pd.DataFrame({
    "side": np.where(report_hit_right, "right", "left"),
    "steps": report_durations,
})
summary = walks.groupby("side")["steps"].agg(
    count="size", mean="mean", median="median", maximum="max"
)
cutoff = walks["steps"].quantile(0.99)
long_walks = walks.query("steps >= @cutoff")
long_counts = long_walks.groupby("side").size()

assert len(walks) == report_n
assert summary["count"].sum() == report_n
assert set(summary.index) == {"left", "right"}
assert len(long_walks) >= report_n // 100
assert long_counts.sum() == len(long_walks)

Code — cell 17

display(summary); print("99% cutoff:", cutoff); display(long_counts.rename("count"))

count mean median maximum
side 
left 32523 91.053685 62.0 954
right 17477 90.714196 63.0 920

99% cutoff: 417.0

side
left 334
right 169
Name: count, dtype: int64

Notes — cell 18

### Q3 written answer for part (c) - solution

The 99th-percentile subset was deliberately selected using duration, so its right/left mix is conditional on being a long record. It cannot replace the unconditional right-side proportion computed from all records. Discrete ties may also leave slightly more than one percent of rows.

Code — cell 19

print('FINAL CHECK: all three questions executed without an exception.')

FINAL CHECK: all three questions executed without an exception.