STA510 Statistical modelling and simulation, autumn 2020.
Mandatory assignment 1
Deadline: Friday September 19th at 12:00 (noon, Norwegian time).
Read carefully through the information about the page on Mandatory assignments on Canvas.
Notice in particular that the assignments should be solved individually.
Hand in on Canvas. Submit two files, one pdf-file with a report containing the answers to the
theory questions, and one file including the R-code. The first line of the R-code file should be:
rm(list=ls()) . Check that the R-code file runs before you submit it. Use comments in the
R-code to clearly identify which question each part of the R-code belong to. Also try to add some
comments to explain important parts of the code. The file ending of the R-code file should be .R
or .r. The report can be handwritten and scanned to pdf-file, or written in your choice of text
editor and converted to pdf. Cite the sources you use.
Problems marked with an R should be solved in R, the others are theory questions.
Problem 1:
In this problem, we consider a random variable X with probability density function given by
f(x) = (12 −x4, if − 1 ≤ x ≤ 10, otherwise
(1)
a) Verify that this is a proper probability density function.
Find the mean E(X) and the variance V ar(x).
Also find E.
b) Find the cumulative distribution function F(x) associated with X, and use this to calculate
P(X < 0) and P(−0.5 < X < 0.5).
Find also the inverse cumulative distribution function F
−1
(u) (Hint: You should get a
quadratic equation which generally have two solutions, but only one of these solutions
(which?) gives sensical results here.)
c)R Write a function (with single argument n) which produces n random numbers with density
(1) using the inverse transform method. Ideally, the function should not contain any forloops.
Check that your function produces correct results by comparing a histogram (use argument
probability=TRUE) of many random variables to the density (1).
d) Using a uniform(-1,1) proposal distribution g, do the calculations required to obtain an
accept-reject method for drawing random numbers with density (1).
1
e)R Write a function (with single argument n) which produces n random numbers with density
(1) using the accept-reject method.
Check the results in the same way as in point c.
Which method would you prefer in this case - accept-reject or inverse transform? The
function system.time is useful in this regard.
f)R Check the answers you got for expectations, variance and probabilities in points a) and b)
above by simulation. For E(X), also find a simulation-based 95% confidence interval. Use
either the inverse transform or accept-reject algorithm for generating random numbers with
density (1).
Problem 2:
In this problem, we let X1, X2, . . . , Xn be iid N(µ, σ2
) random variables.
a) Which distribution does A =Pn
i=1 Xi have?
Which distribution does B =1nPn
i=1 Xi have?
Which distribution does C =√n(B−µ)σ
have?
b)R Check the results you got in a) by simulating many independent replications of each of
A, B, C. You may use e.g. n = 100, µ = 1, σ = 2.
Now, think of X1, X2, . . . , Xn as random outcomes of some stochastic simulation, and you interested
in knowing the mean µ which you estimate using B above. Through initial runs of the
simulation, you have worked out that σ = 2 (at least approximately, and we will use this number
in subsequent calculations). To get the desired precision in the estimate B, you choose n so that
a 95% confidence interval is no longer than 0.1.
c) How many simulations n do you need to obtain such a level of precision in the mean estimate?
d)R For the n obtained in c) and µ = 1, simulate many outcomes of X1, X2, . . . , Xn. Based on
these outcomes, compute approximate sample mean confidence intervals on the form
where S is the sample standard deviation (e.g. obtained using R function sd()). How often
does the interval contain the true mean? Redo this calculation with n = 10. What do you
find?
Problem 3:
In this last problem, we consider the Kissler et al (2020) SEIR model. Please read the documentation
(available on Canvas) of the model carefully, and download the appropriate R-file
seir model.R before starting to solve this exercise. Throughout this problem, we consider a
population of one million persons - be sure to scale results appropriately.
2
a)R Using the default parameters, i.e. by running the model using
> f <- seir sim()
find the following quantities from the output f:
-The number of persons in state S (susceptible) on March 15th 2021.
-The maximum number of persons in critical care (CC) over the period, and also on which
date this occurs.
-The average number of persons in critical care in year 2021.
b)R In this point, we turn of social distancing by using the argument upper thresh=10000, i.e.
> f2 <- seir sim(upper thresh=10000)
-Find the maximum number of persons in critical care over the period.
-Using the definition that an epidemic ends when fewer than 1 per 1 million persons are in
the E-state, how long does it take to reach the end of the epidemic? (be careful here, as
the E-state is also zero before the epidemic starts...)
-Which proportion of the population has been infected (i.e. in RR, RH or RC) when the
epidemic ends in this model (herd immunity)?
Now suppose you are uncertain of the reproduction number R0 and the influence of seasonality
on the infection dynamics. To incorporate this uncertainty, you model the relevant parameters as
follows:
• R0max ∼ uniform(1.8, 3.0).
• season ampliude ∼ uniform(0, 0.4).
R0max and season ampliude are assumed to be independent. The remaining parameters are set
equal to their default values.
c)R Simulate many (at least 100) SEIR models with different random values of R0max and
season ampliude drawn from the above distributions. Based on these simulations; find the
following:
-The mean and standard deviation, the 90%, 95% and 99% quantiles of the maximum
number of required critical care beds during 2020.
-Calculate the probability that more than 100 critical care beds per 1 million persons are
required.
-Make a graphical representation of the maximum number of required critical care beds
during 2020. Is this distribution well approximated by a normal distribution?
d)R -Make scatterplots of the maximum number of critical care beds required in 2020 against
R0max and season ampliude.
-Would you say there is a linear relation between the required critical care beds and the
two parameters.
-Which of these parameters are influencing the required number of critical care beds?
-Calculate the correlations between the required number of critical care beds and the two
parameters.