EnE Modelling and Advanced Data Science Pre-work: Solutions

This notebook contains detailed solutions to the corresponding Exercises notebook.

Tip

Try each question first in the Exercises notebook. Use these solutions to check your reasoning and fill gaps.

---

1. Functions & graphs

Exercise 1 — Solution

Let the radial speed be \(v\) (m/day). Then the radius after \(t\) days is

\[ r(t)=vt. \]

Assuming the infected region is approximately circular, the area is

\[ A(t)=\pi r(t)^2=\pi (vt)^2 = \pi v^2 t^2. \]

This is a polynomial in \(t\) of degree 2 (constant \(\pi v^2\) times \(t^2\)).

For \(v=3\),

\[ A(t)=\pi(3t)^2 = 9\pi t^2. \]

So

\[ A(2)=9\pi(4)=36\pi,\quad A(4)=9\pi(16)=144\pi,\quad A(8)=9\pi(64)=576\pi\ \text{m}^2. \]

If time is measured in weeks \(t_w=t/7\rightarrow t=7t_w\), then

\[ A(t_w)=\pi v^2 (7t_w)^2 = 49\pi v^2 t_w^2. \]

If you instead express speed in m/week, \(v_w=7v\), then \(A(t_w)=\pi (v_w t_w)^2\) as usual.

Exercise 2 — Solution

A power law \(Y=kX^\alpha\) increases with \(X\) when \(\alpha>0\) and decreases when \(\alpha<0\).

• Here \(L=kD^{1.84}\) has exponent \(1.84>0\), so leaf area increases with stem diameter.

• \(S=cT^{-0.49}\) has exponent \(-0.49<0\), so volume fraction decreases with thickness.

Log–log transform:

\[ \log Y = \log k + \alpha \log X. \]

Thus plotting \(\log Y\) vs \(\log X\) gives a straight line whose slope is \(\alpha\) and intercept is \(\log k\).

With \(k=1\):

\[ L(1)=1^{1.84}=1,\quad L(2)=2^{1.84},\quad L(4)=4^{1.84}. \]

With \(c=1\):

\[ S(1)=1^{-0.49}=1,\quad S(2)=2^{-0.49},\quad S(4)=4^{-0.49}. \]

(Use Python to evaluate the non-integer powers.)

Exercise 3 — Solution

Start with

\[ N(t)=N_0\left(\tfrac12\right)^{t/h}. \]

Use \(a^x=e^{x\ln a}\):

\[ \left(\tfrac12\right)^{t/h} = e^{(t/h)\ln(1/2)} = e^{-(t/h)\ln 2}. \]

So

\[ N(t)=N_0 e^{-\mu t}\quad\text{with}\quad \mu = \frac{\ln 2}{h}. \]

For \(h=6\) h, \(\mu=\ln2/6\). Then:

• \(t=3\): \(N(3)=N_0(1/2)^{1/2}=N_0/\sqrt{2}\)

• \(t=6\): \(N(6)=N_0(1/2)=N_0/2\)

• \(t=12\): \(N(12)=N_0(1/2)^2=N_0/4\)

Biologically, \(\mu\) is the instantaneous mortality/decay rate (per hour here). Larger \(\mu\) means faster decline.

Exercise 4 — Solution

In

\[ Z(t)=Z_0 + A\cos(\omega t + \phi), \]

• \(Z_0\): baseline (mean) abundance

• \(A\): amplitude (half the peak-to-trough range)

• \(\omega\): angular frequency (how fast it oscillates)

• \(\phi\): phase shift (when peaks occur)

Period \(T\) and \(\omega\) satisfy \(\omega=2\pi/T\). For \(T=12\) months:

\[ \omega = \frac{2\pi}{12}=\frac{\pi}{6}\ \text{rad/month}. \]

Doubling \(\omega\) halves the period (more cycles per time), so the series oscillates twice as fast. The plot illustrates this.

---

2. Calculus (limits, differentiation, integration, Taylor series)

370623

Exercise 1 — Solution

By definition,

\[ r(t)=\lim_{\Delta t\to 0}\frac{1}{N(t)}\frac{N(t+\Delta t)-N(t)}{\Delta t} = \frac{1}{N(t)}\frac{dN}{dt}, \]

since the limit is exactly the derivative \(dN/dt\).

If \(N(t)=N_0e^{rt}\), then

\[ \frac{dN}{dt}=rN_0e^{rt}=rN(t). \]

Thus

\[ r(t)=\frac{1}{N(t)}\frac{dN}{dt}=\frac{1}{N(t)}(rN(t))=r, \]

a constant: per-capita growth does not change through time in exponential growth.

If the population is decreasing, \(dN/dt<0\) so \(r(t)<0\).

Exercise 2 — Solution

Given

\[ U(r)=\frac{U_{\max}r}{K+r}, \]

use the quotient rule (or rewrite \(U=U_{\max}\frac{r}{K+r}\)). Then

\[ \frac{d}{dr}\left(\frac{r}{K+r}\right) =\frac{(K+r)\cdot 1 - r\cdot 1}{(K+r)^2} =\frac{K}{(K+r)^2}. \]

So

\[ \frac{dU}{dr}=U_{\max}\frac{K}{(K+r)^2}>0, \]

since all parameters are positive.

As \(r\to\infty\),

\[ U(r)=U_{\max}\frac{r}{K+r}\to U_{\max}\frac{1}{1+0}=U_{\max}. \]

For half-saturation,

\[ \frac{U_{\max}r}{K+r}=\frac{U_{\max}}{2} \rightarrow \frac{r}{K+r}=\frac12 \rightarrow 2r=K+r \rightarrow r=K. \]

So \(K\) is the resource level at which uptake is half maximal.

Exercise 3 — Solution

Integrate over one year:

\[ \int_0^1 (a+b\cos(\omega t))\,dt = a\int_0^1 dt + b\int_0^1 \cos(\omega t)\,dt = a + b\left[\frac{1}{\omega}\sin(\omega t)\right]_0^1 = a + \frac{b}{\omega}\big(\sin\omega-\sin 0\big). \]

If the period is one year, \(\omega=2\pi\), then \(\sin(2\pi)=0\) and the cosine term contributes 0, so total production is \(a\).

Average over the year is total divided by 1 year, so it is \(a\) (when the cosine completes an integer number of periods). In general, any pure cosine/sine integrates to zero over whole periods because positive and negative lobes cancel.

Exercise 4 — Solution

For small \(N/K\), the factor \(1-\frac{N}{K}\) is close to 1, so the logistic equation becomes approximately

\[ \frac{dN}{dt}\approx rN, \]

i.e. early-time growth is approximately exponential (density dependence negligible).

For logs, use the Taylor series

\[ \ln(1+x)=x-\frac{x^2}{2}+\frac{x^3}{3}-\cdots,\quad |x|<1. \]

Set \(x=-N/K\) (small in magnitude when \(N\ll K\)):

\[ \ln\left(1-\frac{N}{K}\right) \approx -\frac{N}{K}-\frac{1}{2}\left(\frac{N}{K}\right)^2-\cdots. \]

Biologically: valid at low density relative to carrying capacity, when feedbacks like competition are weak.

---

3. Linear algebra (matrices, eigenvalues/eigenvectors)

370623

Exercise 1 — Solution

The update

\[\begin{split} \begin{pmatrix} J_{t+1}\\ A_{t+1}\end{pmatrix} = \begin{pmatrix} 0 & f\\ s & g\end{pmatrix} \begin{pmatrix} J_t\\ A_t\end{pmatrix} \end{split}\]

means:

\[ J_{t+1}=0\cdot J_t + fA_t = fA_t,\qquad A_{t+1}=sJ_t + gA_t. \]

Interpretations:

• \(f\): fecundity (new juveniles per adult per time step)

• \(s\): maturation/survival from juvenile to adult per step

• \(g\): adult survival (adults remaining adults)

With \(f=3\), \(s=0.2\), \(g=0.8\), \(J_t=10\), \(A_t=5\):

\[ J_{t+1}=3\cdot 5=15,\quad A_{t+1}=0.2\cdot 10 + 0.8\cdot 5 = 2+4=6. \]

Python should return \([15,6]\).

Exercise 2 — Solution

Eigenvalues \(\lambda\) solve \(\det(M-\lambda I)=0\). For

\[\begin{split} M=\begin{pmatrix}0&3\\0.2&0.8\end{pmatrix}, \end{split}\]

\[\begin{split} \det\begin{pmatrix}-\lambda&3\\0.2&0.8-\lambda\end{pmatrix} = (-\lambda)(0.8-\lambda) - 0.6 = \lambda^2 - 0.8\lambda - 0.6=0. \end{split}\]

Solve:

\[ \lambda=\frac{0.8\pm\sqrt{0.8^2+4\cdot 0.6}}{2} =\frac{0.8\pm\sqrt{3.04}}{2}. \]

The larger (dominant) eigenvalue determines long-run growth because \(M^t\) applied to almost any initial vector aligns with the dominant eigenvector and scales like \(\lambda_{\max}^t\).

An eigenvector \(v\) for \(\lambda_{\max}\) gives the stable stage proportions (after normalising so entries sum to 1).

Exercise 3 — Solution

Write as \(A\mathbf{z}=\mathbf{b}\) with

\[\begin{split} A=\begin{pmatrix}3&-7\\1&7\end{pmatrix},\quad \mathbf{z}=\begin{pmatrix}x\\y\end{pmatrix},\quad \mathbf{b}=\begin{pmatrix}4\\10\end{pmatrix}. \end{split}\]

Compute determinant:

\[ \det(A)=3\cdot 7 - (1)(-7)=21+7=28\neq 0, \]

so \(A\) is invertible. For a \(2\times2\) matrix,

\[\begin{split} A^{-1}=\frac{1}{\det(A)}\begin{pmatrix}7&7\\-1&3\end{pmatrix} =\frac{1}{28}\begin{pmatrix}7&7\\-1&3\end{pmatrix}. \end{split}\]

Then

\[\begin{split} \mathbf{z}=A^{-1}\mathbf{b} =\frac{1}{28}\begin{pmatrix}7&7\\-1&3\end{pmatrix}\begin{pmatrix}4\\10\end{pmatrix} =\frac{1}{28}\begin{pmatrix}28+70\\-4+30\end{pmatrix} =\begin{pmatrix}98/28\\26/28\end{pmatrix} =\begin{pmatrix}7/2\\13/14\end{pmatrix}. \end{split}\]

So \(x=3.5\), \(y\approx 0.9286\).

Exercise 4 — Solution

Euler’s formula:

\[ e^{i\beta t}=\cos(\beta t)+i\sin(\beta t). \]

So if a linear system has solutions involving \(e^{(\alpha+i\beta)t}\),

\[ e^{(\alpha+i\beta)t}=e^{\alpha t}e^{i\beta t} =e^{\alpha t}\big(\cos(\beta t)+i\sin(\beta t)\big), \]

whose real part is \(e^{\alpha t}\cos(\beta t)\): an exponentially growing/decaying oscillation.

For \(x''+\beta^2x=0\), try \(x=e^{\lambda t}\):

\[ \lambda^2 e^{\lambda t}+\beta^2 e^{\lambda t}=0 \rightarrow \lambda^2+\beta^2=0 \rightarrow \lambda=\pm i\beta. \]

Thus the general real solution is \(x(t)=C_1\cos(\beta t)+C_2\sin(\beta t)\).

---

4. Differential equations (equilibria, stability, phase-plane intuition)

370623

Exercise 1 — Solution

Solve \(dN/dt=rN\):

\[ \frac{1}{N}dN = r\,dt\rightarrow \ln N = rt + C \rightarrow N(t)=N_0e^{rt} \]

using \(N(0)=N_0\).

Doubling time \(T_d\) satisfies \(N(T_d)=2N_0\):

\[ 2N_0=N_0e^{rT_d}\rightarrow 2=e^{rT_d}\rightarrow T_d=\frac{\ln 2}{r}. \]

For \(r=0.4\), \(T_d=\ln2/0.4\approx 1.7329\) hours.

Exercise 2 — Solution

Equilibria satisfy \(dN/dt=0\):

\[ rN\left(1-\frac{N}{K}\right)=0 \rightarrow N=0 \quad\text{or}\quad 1-\frac{N}{K}=0\rightarrow N=K. \]

Stability via sign analysis:

• If \(0<N<K\), then \(1-N/K>0\) so \(dN/dt>0\): trajectories increase away from 0 ⇒ \(N=0\) is unstable.

• If \(N>K\), then \(1-N/K<0\) so \(dN/dt<0\): trajectories decrease toward \(K\).
  Thus \(N=K\) is stable.

Sketch: \(dN/dt\) is a concave-down parabola in \(N\) crossing zero at \(0\) and \(K\), positive between them, negative above \(K\).

Exercise 3 — Solution

Compute eigenvalues of

\[\begin{split} A=\begin{pmatrix}0&1\\-2&-0.4\end{pmatrix}. \end{split}\]

They solve \(\det(A-\lambda I)=0\):

\[\begin{split} \det\begin{pmatrix}-\lambda&1\\-2&-0.4-\lambda\end{pmatrix} =\lambda(0.4+\lambda)+2 =\lambda^2+0.4\lambda+2=0. \end{split}\]

Discriminant \(\Delta=0.4^2-8<0\), so eigenvalues are complex with negative real part (\(-0.2\)). Therefore the equilibrium is a stable spiral: trajectories oscillate while decaying toward the origin.

Optional simulation would show spiralling trajectories in the \((x,y)\) plane.

Exercise 4 — Solution

Solve

\[ \frac{dN}{dt}-rN=-H. \]

Integrating factor: \(I(t)=e^{-rt}\). Multiply both sides:

\[ e^{-rt}\frac{dN}{dt}-re^{-rt}N = -He^{-rt} \rightarrow \frac{d}{dt}\big(e^{-rt}N\big)=-He^{-rt}. \]

Integrate:

\[ e^{-rt}N = \int -He^{-rt}\,dt + C = \frac{H}{r}e^{-rt}+C. \]

Multiply by \(e^{rt}\):

\[ N(t)=\frac{H}{r}+Ce^{rt}. \]

Use \(N(0)=N_0\): \(N_0=H/r + C\rightarrow C=N_0-H/r\). So

\[ N(t)=\frac{H}{r}+\left(N_0-\frac{H}{r}\right)e^{rt}. \]

Equilibrium: set \(dN/dt=0\rightarrow rN-H=0\rightarrow N^*=H/r\). Feasible if \(H/r>0\) (always, given \(H,r>0\)), but note the meaning differs: this is the level where growth balances harvest.

Non-crossing of zero: because \(e^{rt}\) grows, if \(N_0<H/r\) then \(C<0\) and \(N(t)\) decreases without bound (eventually negative in the model, signalling extinction/invalidity). For \(N(t)\ge 0\) for all \(t\ge 0\) you need \(C\ge 0\rightarrow N_0\ge H/r\).

---

5. Probability & distributions

370623

Exercise 1 — Solution

For a discrete random variable \(X\) with pmf \(P(X=x)\),

\[ E[X]=\sum_x x\,P(X=x). \]

If \(X\sim\text{Poisson}(\lambda)\), then

\[ E[X]=\lambda,\qquad \mathrm{Var}(X)=\lambda. \]

With \(\lambda=1.5\),

\[ P(X=0)=e^{-\lambda}=e^{-1.5}, \]

and

\[ P(X\ge2)=1-P(0)-P(1)=1-e^{-1.5}-1.5e^{-1.5}. \]

Exercise 2 — Solution

If \(Y\sim\text{Binomial}(n,p)\),

\[ E[Y]=np,\qquad \mathrm{Var}(Y)=np(1-p). \]

For \(n=50\), \(p=0.1\):

\[ E[Y]=5,\quad \mathrm{Var}(Y)=50(0.1)(0.9)=4.5. \]

Also

\[ P(Y=0)=(1-p)^n=0.9^{50}. \]

Interpretation: probability a survey of 50 hosts finds no infections even though prevalence is 10%—useful for planning sample sizes.

Exercise 3 — Solution

Conditional density:

\[ L_m\mid L \sim \mathcal{N}(L,\sigma^2), \]

so

\[ f(L_m\mid L)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(L_m-L)^2}{2\sigma^2}\right). \]

If you average \(k\) independent replicates, the average error has variance \(\sigma^2/k\). Thus replication reduces variance and improves precision by averaging out noise.

Exercise 4 — Solution

A joint distribution specifies probabilities for the four outcomes \((A,B)\in\{0,1\}^2\) that sum to 1, e.g.

\[ P(0,0),P(0,1),P(1,0),P(1,1)\ge 0,\quad \sum P=1. \]

Marginal:

\[ P(A=1)=P(1,0)+P(1,1). \]

Conditional:

\[ P(A=1\mid B=1)=\frac{P(1,1)}{P(0,1)+P(1,1)}. \]

Interpretation: compares presence of A in patches where B is present; deviations from \(P(A=1)\) suggest positive/negative association (though causality requires more modelling).

---

6. Optimisation & likelihood

370623

Exercise 1 — Solution

Write

\[ S(\beta_0,\beta_1)=\sum_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2. \]

Differentiate:

For \(\beta_0\):

\[ \frac{\partial S}{\partial \beta_0} =\sum 2(y_i-\beta_0-\beta_1x_i)(-1) =-2\sum (y_i-\beta_0-\beta_1x_i). \]

Set to 0:

\[ \sum (y_i-\beta_0-\beta_1x_i)=0 \rightarrow n\beta_0+\beta_1\sum x_i=\sum y_i. \]

For \(\beta_1\):

\[ \frac{\partial S}{\partial \beta_1} =\sum 2(y_i-\beta_0-\beta_1x_i)(-x_i) =-2\sum x_i(y_i-\beta_0-\beta_1x_i). \]

Set to 0:

\[ \sum x_i(y_i-\beta_0-\beta_1x_i)=0 \rightarrow \beta_0\sum x_i + \beta_1\sum x_i^2 = \sum x_i y_i. \]

These are the normal equations (a linear system) giving the minimiser. Geometrically: search over a 2D parameter space \((\beta_0,\beta_1)\) for the point where the “error surface” is lowest.

Exercise 2 — Solution

Poisson likelihood for i.i.d. data \(x_1,\dots,x_n\):

\[ L(\lambda)=\prod_{i=1}^n \frac{e^{-\lambda}\lambda^{x_i}}{x_i!} = e^{-n\lambda}\lambda^{\sum x_i}\prod\frac{1}{x_i!}. \]

Log-likelihood:

\[ \ell(\lambda)=\log L(\lambda)= -n\lambda + \left(\sum x_i\right)\log\lambda - \sum \log(x_i!). \]

Differentiate:

\[ \ell'(\lambda)=-n + \frac{\sum x_i}{\lambda}. \]

Set to 0:

\[ -n + \frac{\sum x_i}{\lambda}=0 \rightarrow \hat\lambda=\frac{1}{n}\sum x_i = \bar{x}. \]

The “height” of the likelihood curve at a given \(\lambda\) indicates how plausible that parameter value makes the observed data (relative to other \(\lambda\)’s). For the sample \([0,1,0,2,1]\), \(\hat\lambda=\bar{x}=0.8\).

Exercise 3 — Solution

For Normal(\(\mu,\sigma^2\)) i.i.d.:

\[ L(\mu,\sigma)=\prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y_i-\mu)^2}{2\sigma^2}\right). \]

Log-likelihood:

\[ \ell(\mu,\sigma)= -n\log\sigma - \frac{1}{2\sigma^2}\sum (y_i-\mu)^2 + \text{constant}. \]

mles:

\[ \hat\mu=\bar{y},\qquad \hat\sigma^2=\frac{1}{n}\sum (y_i-\bar{y})^2 \]

(for mle; note unbiased estimator uses \(n-1\) in the denominator).

Because \(\ell\) depends on two parameters, you can view it as a surface over the \((\mu,\sigma)\) plane; contours are “slices” of equal log-likelihood, helping you see ridges/curvature (identifiability, uncertainty, etc.).

Exercise 4 — Solution

Examples of positive parameters: mortality rates, growth rates, diffusion coefficients, variances, concentrations.

If a decay model \(N(t)=N_0e^{-\mu t}\) is fitted with \(\mu<0\), then \(-\mu t\) is positive and the model predicts exponential growth, contradicting “decay” and likely indicating a sign/parameterisation issue.

reparameterise \(\mu=e^\theta\) with unconstrained \(\theta\in\mathbb{R}\). Then \(\mu>0\) automatically. Optimisation becomes over \(\theta\) rather than \(\mu\); the likelihood or error function is composed with \(\mu(\theta)=e^\theta\), which can improve stability and ensures biological feasibility.