21/4/2023 20:58MAST90045 Assignment 2
第1/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
MAST90045 Assignment 2
Due on 1 May 2023
Question 1 (50 pts)
Imagine you are an entomologist studying new species on two remote islands (island A
and island B) to develop new taxonomy of related species. You hire local people to catch
insects and take measurements of their body lengths. Each measurement is assumed to
be independently collected.
You have received 100 measurements of a new species from the island A. A histogram
looks as follows.
21/4/2023 20:58MAST90045 Assignment 2
第2/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
Given the histogram, you have decided to use a normal distribution as
a probabilistic model to describe the body length distribution. Let the data be denoted by
, which is stored in X.A in the answer source Ple.
In settings like this, the standard approach to specifying the model (i.e., specifying
and in this case) is maximum likelihood estimation (MLE). The idea of MLE is
that, given the parameterised probabilistic model, we search for a parameter that most
likely gives rise to the observed data. Formally, the best parameter is a solution to
the following optimisation problem:
where
is the PMF/PDF for each of independent and identically distributed (IID)
random variables parameterised by , and
is the corresponding observed data.
∼ N( , )XA μ1 σ21
{xAn }100n=1
μ1
σ1
( )θ∗
∈ f ( ; θ),θ∗ argmax
θ∈Θ ∏n=1
N
xn
f N
{Xn}Nn=1 θ
{xn}Nn=1
21/4/2023 20:58MAST90045 Assignment 2
第3/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
The objective function
is a function of and called the likelihood function. It is the product of function values
obtained by evaluating the PMF/PDF at for . The multiplication is justiPed by
the independence of . Note that
and the left-hand side is typically easier to work with in both analytical and numerical
optimisation.
In this question, you are asked to practise using R built-in optimiser optim with
method="L-BFGS-B" , a popular choice among applied researchers. Read its help page
and learn how to use it, especially how to use lower and upper parameters to restrict
the search space .
1-1 (3 pts)
Back to the island A. Let where (i.e., ).
Now, based on the modelling assumption , formulate a maximisation
problem.
1-2 (4 pts)
Find a local maximiser and round it to one decimal place. Playing with par , lower ,
and upper based on the histogram helps you choose good values.
1-3 (3 pts)
Moving on to the island B. You have received 300 measurements of new species from
the island B. A histogram looks as follows.
f ( ; θ)∏
n=1
N
xn
θ N
{xn}Nn=1 θ
{Xn}Nn=1
log( f ( ; θ)) = f ( ; θ),argmaxθ∈Θ ∏n=1
N
xn argmax
θ∈Θ ∏n=1
N
xn
Θ
= ( , ) ∈θA μ1 σ1 ΘA = ℝ ×ΘA ℝ++ > 0σ1
∼ N( , )XA μ1 σ21
θÂ
21/4/2023 20:58MAST90045 Assignment 2
第4/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
While you are certain that there are two new species (species 2 and 3) whose size
distributions are distinct, local people were unable to tell the difference and recorded
only measurements of lengths. Encouraged by the above histogram, your conviction that
the data came from two distinct species leads to the following model.
The new species’ body length distribution consists of two normal distributions—it is
drawn from with probability and drawn from with probability
where . In other words, the PDF for this random variable is
As for the island A model, using the above PDF, your approach to Pnding the best
parameter is MLE. Note that where
. Formulate a maximisation problem.
N( , )μ2 σ22 q N( , )μ3 σ23
1 − q q ∈ (0, 1) XB
f ( ) = exp(− ( − ) + exp(− ( − ).xB q2πσ22‾ ‾‾‾‾√
1
2σ22
xB μ2)2
1 − q
2πσ23‾ ‾‾‾‾√
1
2σ23
xB μ3)2
θB∗ = ( , , , , q) ∈θB μ2 σ2 μ3 σ3 ΘB
= ℝ × × ℝ × × (0, 1)ΘB ℝ++ ℝ++
21/4/2023 20:58MAST90045 Assignment 2
第5/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
1-4 (6 pts)
Try to Pnd a local maximiser in the above problem using optim . Depending on your
speciPcation and creativity, you may or may not be able to solve it. If you make it, report
rounded to one decimal place. If you cannot, make a sensible conjecture about the
cause of the challenge. Either case, remember to explain your answer.
1-5 (3 pts)
You have got an idea about an alternative approach to the problem. Although is a
deterministic parameter, it is for the probability that a body length is drawn from species
2 (and from species 3). So, why don’t we treat the species selection as another
random event? SpeciPcally, dePne a new random variable such that
where is the PMF of . Then, we have the size distribution conditional on as
follows:
Find an expression for , the joint PDF for a mix of continuous and discrete .
Remember to be explicit about where and where .
1-6 (4 pts)
Marginalise over to recover the marginal PDF and marginalise
over to recover the marginal PMF .
1-7 (7 pts)
θB̂
θB̂
q
1 − q
S
p(2) = q and p(3) = 1 − q,
p(s) S S = s
f ( |s) = exp(− ( − ).xB 12πσ2s‾ ‾‾‾‾√
1
2σ2s
xB μs)2
f ( , s)xB XB S
f ( , s) > 0xB f ( , s) = 0xB
f ( , s)xB S f ( )xB f ( , s)xB
XB p(s)
21/4/2023 20:58MAST90045 Assignment 2
第6/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
At this point, you come to regret not telling the local people how to distinguish two
species so that they can also record the species for each measurement. Why? Explain in
detail how the knowledge of species for each measurement could make your MLE
easier.
Note that, the independence of each observation implies that the likelihood of given
by the joint PDF for all 300 pairs of random variables would be simply the
product of each component PDF :
1-8 (5 pts)
You are tenacious and start thinking along the following lines. Now, the species selection
is a random variable, and the knowledge of species for each measurement could make
the MLE easier. Although is unmeasured, we can derive its distribution as a function of
. So, to eliminate the dependency on , why don’t we take the expectation of
the MLE objective function with respect to and then maximise it. In particular,
to make use of the measurement data , we should take the expectation with
respect to the PMF of conditional on ,
where and .
Derive an expression for .
1-9 (8 pts)
Given the derived , here is the formal description of the above algorithm. In
search of ,
0. Begin with some initial .
1. With , specify .
2. Find a local maximiser by solving
θB
{( , )XBn Sn }300n=1
f ( , ; )xBn sn θB
f ( , ; ).∏
n=1
300
xBn sn θB
S
S
θB {sn}300n=1
{Sn}300n=1
{xBn }300n=1
{Sn}300n=1 {xBn }300n=1
p(s| ; ),xB θB
s = ( , ,… , )s1 s2 s300 = ( , ,… , )xB xB1 xB2 xB300
p(s| ; )xB θB
p(s| ; )xB θB
θB∗
θBold
θBold p(s| ; )xB θBold
θBnew
피[log( f ( , ; ))] ,
300
n
B
21/4/2023 20:58MAST90045 Assignment 2
第7/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
where the expectation is taken with respect to .
3. If , go back to #1 and repeat the process using in place
of ; Otherwise, you are done.
Before implementing the algorithm, you need to make Step #2 more explicit. Using
found in Q1-8, derive an expression for the objective
. Your expression should be in terms of , , , , ,
, and and refect your answer to Q1-7. It may be a long chain of
equalities, so make sure to explain key steps of your derivation to demonstrate your
understanding.
1-10 (7 pts)
By completing the given code, implement the algorithm and report rounded to one
decimal place. As usual, to demonstrate your understanding, make sufPcient comments
on your code.
Question 2 (50 pts)
In this question, you are asked to compare performance of two ODE solvers, RK4 and
adaptive RK4 (ARK4), in long-term prediction using the predator-prey model studied in
the lecture. Long-term prediction simply means large time horizon . By working on a
problem that consists of many component processes and requires non-trivial
computation on the whole, you will get a feel of large simulations used in practice.
To measure performance of a solver, we dePne closeness of a simulated trajectory to the
true orbit. First, recall that the orbit is dePned by the equation
where , , , , and . You may assume
the orbit is contained in .
피[log( f ( , ; ))] ,argmax∈θB ΘB ∏n=1
300
xBn Sn θB
p(s| ; )xB θBold
‖ − > ϵθBnew θBold‖∞ θBnew
θBold
p(s| ; )xB θB
피 [log(∏ f ( , ; ))]xBn Sn θB μ2 σ2 μ3 σ3 q
θBold , , … ,xB1 xB2 xB300
θB̂
T
c = gb − e log( ) + b − a log( ),y1 y1 y2 y2
a = 0.05 b = 2 × 10−4 g = 0.8 e = 0.03 c = −0.3
R = [10, 810] × [20, 830]
21/4/2023 20:58MAST90045 Assignment 2
第8/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
The closeness of a trajectory is dePned to be the average of all deviations computed at
every visited . The deviation of arbitrary location from the orbit is dePned as
the shortest distance to reach the orbit by moving from to one of four directions: up,
down, left, or right. For example, in the following Pgure, the arrows indicate the
deviations from and the deviation from .
By computing the deviation of for every recorded , you can get the
average deviation of a trajectory.
2-1 (10 pts)
The algorithm to compute the deviation of arbitrary location is constructed as
follows. First, based on the orbit equation
dePne
(t)ŷ ∈ Ryc
yc
= (600, 730)yd = (130, 400)ye
(t)ŷ t ∈ [0, T]
∈ Ryc
c = gb − e log( ) + b − a log( ),y1 y1 y2 y2
21/4/2023 20:58MAST90045 Assignment 2
第9/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
Then, observe:
1. Finding a distance from to the orbit in either vertical or horizontal
direction is equivalent to Pnding a root of with Pxed or
respectively.
2. Let and . Then, both and are strictly
convex functions.
3. For a strictly convex function, regarding roots, there are three possible scenarios:
one root if the global minimum is equal to 0; no root if the global minimum is
greater than 0; and two roots otherwise.
4. In this particular problem, both and are strictly positive if they are evaluated at
a boundary point.
Now, complete the given code that implements deviation , which takes arbitrary
location y contained in and returns the deviation of y . Use bisection provided in
the answer source Ple for numerical root-Pnding.
(While you could use other root-Pnding methods, the mathematical structure of this
particular problem makes it easy to use the bisection method
(https://yujisaikai.com/MAST90045/root_Pnding.html#bisection-method), which has the
advantage of guaranteed convergence.)
2-2 (4 pts)
Compared to the implementation presented in the lecture, RK4 here uses
slightly different parameterisation,
separate K4 function for slope computation, and
solution matrix Y with three columns for and length-2 vector .
Both K4 and RK4 are found in the answer source Ple.
f ( , ) = gb − e log( ) + b − a log( ) − c.y1 y2 y1 y1 y2 y2
= ( , )yc yc1 yc2
f =y1 yc1 =y2 yc2
( ) = f ( , )f1 y1 y1 yc2 ( ) = f ( , )f2 y2 yc1 y2 f1 f2
f1 f2
R
t ŷ
21/4/2023 20:58MAST90045 Assignment 2
第10/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
The trajectory generated by RK4 with , , and
looks as follows. (n.b. the choice of , and is arbitrary; it
produces an illustrative Pgure.)
K4 <- function(f, y, h) {
k1 <- f(y)
k2 <- f(y + h/2*k1)
k3 <- f(y + h/2*k2)
k4 <- f(y + h*k3)
return(k1/6 + k2/3 + k3/3 + k4/6)
}
RK4 <- function(f, y0, T, N) {
h <- T/N
Y <- matrix(c(0,y0), nrow=1) # initialised with y0 at t=0
for (n in 1:N) {
y.hat <- Y[n,2:3] + h*K4(f, Y[n,2:3], h)
Y <- rbind(Y, c(h*n, y.hat)) # append for t=h*n
}
return(Y)
}
y(0) = = (150, 803.9929)y0 T = 1176
N = 100 T = 1176 N = 100
21/4/2023 20:58MAST90045 Assignment 2
第11/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
Write code that computes the average deviation of this trajectory. Print what you get.
2-3 (10 pts)
Similar to adaptive quadrature
(https://yujisaikai.com/MAST90045/numerical_integration.html#adaptive-quadrature),
ARK4 continues to halve subintervals until the stopping criterion is locally met for each
subinterval. SpeciPcally, Prst, initialise ARK4 by dividing into evenly-space
subintervals. Then, for each subinterval (n.b., values of and depend on a
speciPc subinterval to work on), repeat the following.
1. Compute by taking a single step of size from .
2. Compute by taking double steps (i.e. two steps) of size
from .
3. Examine whether .
4. If the condition is met, stop halving, take the double steps, and go to the next
subinterval. Otherwise, go back to #1 for each of two subintervals
T N0
[ , ]ta tb ta tb
( )ŷs tb h = −tb ta ( )ŷ ta
( )ŷd tb h = ( − )/2tb ta
( )ŷ ta
‖ ( ) − ( ) ≤ ϵŷd tb ŷs tb ‖∞
[ , ( + )/2]ta ta tb
21/4/2023 20:58MAST90045 Assignment 2
第12/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
and as new .
Note:
As for RK4 above, use K4 function to compute slopes.
When taking a double step, remember to update the slope after the Prst half step
and use it for the second half step.
When halving is halted, make sure to append and from both steps to the
solution table Y .
Now, by completing the given code, implement ARK4 using ARK4.recur recursively.
To help you check your code and gain insight into how ARK4 works, two plotting
functions are provided in the answer source Ple: ploty and ploty_num . These
functions take a solution matrix Y and plots a trajectory. In particular, ploty_num plots
a trajectory with explicit at each . It will be legible only for , but nothing
stops you from experimenting different and . While is typically part of the
problem, and are algorithm parameters balancing accuracy and computational
requirement. (We Px in this assignment.)
2-4 (14 pts)
1. Run ARK4 for , , and .
2. Compute the average deviation of the resulting trajectory.
3. Report the smallest and largest steps in the trajectory as well as the locations
where these steps were taken. Note that the length of an actual step/movement is
the Euclidean norm of h*K4 , the step size h multiplied by the slope vector
computed using K4 .
4. Find how many steps RK4 needs to expend to achieve a similar average deviation
(use < 0.1% difference between two deviations as the dePnition of being similar).
Comment on your Pnding.
Note that is just a reasonable value chosen based on the following
heuristics. We know a single round roughly takes . As a crude initialisation,
start with expending four evenly-spaced step to make a single round (i.e., taking four
turns). Hence, divide 30000 by 200 and multiply it by 4 to get 600.
2-5 (12 pts)
[( + )/2, ]ta tb tb [ , ]ta tb
t (t)ŷ
t (t)ŷ T ≤ 200
T N0 T
N0 ϵ
ϵ = 0.05
y(0) = = (150, 803.9929)y0 T = 30000 = 600N0
= 600N0
Δt = 200
21/4/2023 20:58MAST90045 Assignment 2
第13/13⻚file:///Users/bonnie/Downloads/asmt2_q.html
Finding optimal inputs to simulators is a common problem in many applications. In our
case, a key input is .
Come up with an algorithm to Pnd the best that minimises the average
deviation of ARK4 trajectory for , and if starting with .
Write code that implements the algorithm. Your code should outputs and the
corresponding average deviation. Explain your search algorithm and its implementation.
(Remember you have the time budget of 10 minutes for the entire assignment, so the
search effort should not be overly exhaustive.)
∈ Ry0
∈ Ry0∗
T = 30000 = 600N0 y(0) = y0∗
y0∗