1 Monte Carlo simulations 1: the basics

Note

Hi there student risk manager! If you haven’t yet done so, please first have a quick glance through About this course and About the organisation and format of these notes to set the stage for this course and some useful info about how this set of online course materials works.

1.1 Why??

You as risk manager are facing a risk i.e. a loss of uncertain size, modelled by a (non-negative) random variable. Or maybe even multiple risks/random variables. One thing you at least want to be able to do, is to compute stuff like (joint) probabilities, expectations etc. You’ll remember from earlier Prob courses the standard tools for this: an integral involving the pdf for continuous random variables, a summation involving the pmf for discrete random variables — see e.g. Section A.1.2 and Section A.1.3 for a reminder of these key formulas if you’d like. But what if this brings you to an integral or summation that is (way) too complicated to evaluate algebraically? Or what if you’re dealing with a random variable that is in some way constructed in terms of other random variables (very common in practice) making it very hard/impossible to figure out what its pdf/pmf is?

In such cases, Monte Carlo simulations can be a great tool: a pretty simple yet incredibly useful technique that is very commonly employed in many real world risk management situations!

1.2 Key idea

How do you find out whether a coin that you have found is fair i.e. whether it shows Heads and Tails with equal probability? Well if you’re very bored then you would probably just toss it an awful lot of times, the more the better, and you would count how often you got Heads and Tails. If that results in about $50\%$ Heads and Tails then on balance you’d probably be happy to accept that it’s indeed a fair coin. Buy why actually? Well because we all have this deeply engrained feeling that it should be true that for ‘large’ $n$ \[\begin{align} \P(\text{Heads}) &\approx \frac{\text{Nr of Heads in $n$ tosses}}{n} \quad \text{and} \\ \P(\text{Tails}) &\approx \frac{\text{Nr of Tails in $n$ tosses}}{n}. \end{align} \tag{1.1}\]

Let’s formulate this principle in a slightly more general way. Suppose that we have a rv $X$ that denotes the outcome of some experiment, and that we would repeat this experiment a large $n$ number of times (independently) with observed values $x_1, x_2, \ldots, x_n$ of $X$. Further take some set $B \subseteq \R$, say an interval. Then by exactly the same feeling as above, we would expect that the theoretical probability that $X$ ends up in $B$ is approximately equal to how many of the observed outcomes out of $n$ have ended up in $B$. That is: \[\P(X \in B) \approx \frac{\text{Nr of $x_i$'s that satisfy $x_i \in B$}}{n} \quad \text{for 'large' $n$.} \tag{1.2}\] For instance, with $X$ denoting the random outcome of rolling a fair die, on the one hand we know that $\P(X \leq 2)=1/3$, and on the other hand if we execute this experiment a large $n$ number of times with observed outcomes $x_1, \ldots, x_n$ then we would expect the number of outcomes less than or equal than $2$ makes up about $1/3$ of all outcomes i.e. there should be about $n/3$ no?

Here is another example of that same principle. How to figure out what the expected profit for some new betting game is? Well if you could play it a large $n$ number of times, then your expected profit per game should be about equal to your average profit i.e. your total profit after the $n$ games divided by $n$ no? Or, slightly more generally formulated: if we again denote by the rv $X$ the (random) outcome of some experiment and we repeat it $n$ times (independently) to observe the samples/values $x_1, x_2, \ldots, x_n$ of $X$, then our gut tells us that the theoretical mean of $X$ should be approximately equal to the average of the observed values: \[\E[X] \approx \frac{1}{n} \sum_{i=1}^n x_i \quad \text{for 'large' $n$.} \tag{1.3}\] Or, slightly more general, if rather than $X$ itself we look at a rv $g(X)$ for some (half decent) function $g: \R \to \R$, just because it allows us some extra flexibility in our applications) then the analogue approximation reads \[\E[g(X)] \approx \frac{1}{n} \sum_{i=1}^n g(x_i) \quad \text{for 'large' $n$.} \tag{1.4}\]

So, now what? All of the above seems nothing more than stating some rather obvious facts — let’s accept that these approximations are indeed valid for the moment, we’ll discuss justification and error bounds in Section 1.3 below. Well here is one beautiful use of them you may not have anticipated. Suppose that we find ourselves in the following situation:

it is, for whatever reason, inefficient/hard/impossible to compute $\P(X \in B)$ or $\E[g(X)]$ using your standard tools from earlier Probability courses,
but we can think of a way to produce samples from $X$,

then we can obtain an approximation of $\P(X \in B)$ resp. $\E[g(X)]$ by applying Equation 1.2 resp. Equation 1.4!

The underlying big picture idea is here as follows. If we have many samples $x_1,\ldots,x_n$ from some random variable $X$ available, then that gives us a pretty good idea what the distribution of $X$ is i.e., roughly speaking, what values $X$ tends to take and with what probabilities. For example, if each sample is an integer between $1$ and $6$ and there are about equally many of each of these integers, then that suggests that $X$ has the range $\{1,\ldots,6\}$ and takes all these values with probability $1/6$ (i.e. it could be the outcome of rolling a fair die). Or if these samples look like real numbers and seem pretty much evenly scattered throughout the interval $[0,1]$, then that suggests that $X \sim \text{Unif}(0,1)$. Or, if the samples do not give such an immediately obvious picture, we could for instance put them in a histogram to try and visualise them. If that histogram has a shape that reminds us strongly of the famous bell curve, then that suggests that $X$ has a Normal distribution. Etc¹.

¹ Just for a bit more context (not relevant for us): as you may well have encountered in a past course, in Statistics there exist techniques to make this all a lot more precise. In particular there is a whole battery of “distribution goodness-of-fit” tests that apply the principles of hypothesis testing to determine how much evidence the samples provide that they come from a certain distribution

Further, if these samples allow us to learn (approximately) what the distribution of $X$ is, then it only makes sense that we can also use these samples to (approximately) compute any distributional property of $X$ i.e. any probability or expectation. The only question is then how you would use these samples exactly for the probability/expectation you are interested in, and the obvious way to do that is expressed in Equation 1.2 resp. Equation 1.4.

Definition 1.1 The Monte Carlo simulation technique entails that for a given random variable $X$, we generate $n$ samples (i.e. independent realisations) $x_1, \ldots, x_n$ of $X$ and then apply the approximation \[\P(X \in B) \approx \frac{\text{Nr of $x_i$'s that satisfy $x_i \in B$}}{n} \tag{1.5}\] for a (half decent) set $B \subseteq \R$ and/or \[\E[g(X)] \approx \frac{1}{n} \sum_{i=1}^n g(x_i), \tag{1.6}\] for a (half decent) $g: \R \to \R$.

Note: “half decent” refers to the fact that that set resp. function need to satisfy some technical conditions for the probability resp. expectation to be well defined, however that’s not something you need to worry about for the purposes of this course.

Further, this technique remains equally valid when we are dealing with a random vector i.e. a vector of random variables rather than a single random variable, see Remark 1.1.

On a practical note, this technique is not something that is very suitable for a poor mathematician with only pen and paper: not only do you need to somehow find these samples somewhere, but also actually working out the summations in Equation 1.2 and Equation 1.4 by hand if we’re using say a million samples is not very pleasant! That’s where computers come in, and that also explains why this technique has become really prominent in recent decades only as sufficient electronic computing capacity started to become generally available. And also for producing samples from a given distribution computers are very useful, using software like R or Python for instance.

Remark 1.1. An important category of applications for Monte Carlo simulations is in situations where rather than a single random variable $X$ we are facing a random vector: a vector of random variables, say $(X^{(1)},X^{(2)},\ldots,X^{(k)})$. For instance, a company looking ahead at the coming financial year may very well have to deal with with multiple sources of uncertainty: income per client, outgo per client, number of new clients, etc.

Let’s discuss how the ideas from Definition 1.1 extend from working with a single random variable to working with a random vector $(X^{(1)},X^{(2)},\ldots,X^{(k)})$.

1. Samples of random vectors. Note that where a sample from a single random variable is a real number, a sample of a random vector $(X^{(1)},X^{(2)},\ldots,X^{(k)})$ such a vector is a vector $(x^{(1)},x^{(2)},\ldots,x^{(k)}) \in \R^k$, where for each $i=1,\ldots,k$, $x^{(i)} \in \R$ is a sample from the random variable $X^{(i)}$. Warning: if $X^{(1)},X^{(2)},\ldots,X^{(k)}$ are mututally independent then you can simply produce each $x^{(i)}$ separately and put them together to form a vector $(x^{(1)},x^{(2)},\ldots,x^{(k)})$, however if they are not mututally independent then that is not a valid approach! Studying random vectors with some (non-trivial) dependency structure is a (challenging) topic all on its own, which we’ll return to later in the course when we discuss copulas.

2. Probabilities for random vectors. Because this random vector takes its values in $\R^k$ (each sample/realisation is a point in $\R^k$ i.e. a length $k$ vector of real numbers after all), probabilities (always) take the form of “what is the probability that the random vector ends up in $B \subseteq \R^k$” i.e. you specify some subset $B \subseteq \R^k$ and wonder what \[\P\left( (X^{(1)},X^{(2)},\ldots,X^{(k)}) \in B \right)\] is. Some examples, for $k=2$:

The probability that both random variables are larger than $1$ is $\P((X^{(1)},X^{(2)}) \in B)$ with \[B=\{ (x^{(1)},x^{(2)}) \, | \, x^{(1)}>1 \text{ and } x^{(2)}>1 \}=(1,\infty) \times (1,\infty) \subseteq \R^2\] (recall that the “$\times$” notation is a Cartesian product);
The probablity that $X^{(1)}$ is larger than $X^{(2)}$ is $\P((X^{(1)},X^{(2)}) \in B)$ with \[B=\{ (x^{(1)},x^{(2)}) \, | \, x^{(1)}>x^{(2)} \} \subseteq \R^2;\]
The probability that $(X^{(1)},X^{(2)})$ ends up in the unit circle is $\P((X^{(1)},X^{(2)}) \in B)$ with \[B=\{ (x^{(1)},x^{(2)}) \, | \, (x^{(1)})^2+(x^{(2)})^2 \leq 1 \} \subseteq \R^2,\]

etc. etc. etc. It can be helpful to visualise such a situation: take a plane with an $x^{(1)}$-axis and $x^{(2)}$-axis representing $\R^2$ as the space where $(X^{(1)},X^{(2)})$ takes its values, draw your $B$ in there, and then realise that each sample/realisation of $(X^{(1)},X^{(2)})$ is some point in this space, generally sometimes inside $B$ and sometimes outside of it, and that $\P((X^{(1)},X^{(2)}) \in B)$ is nothing but the probability that a sample/realisation ends up in $B$. See also Figure A.3 e.g.

Now, Equation 1.5 from Definition 1.1 translates as follows: if we produce $n$ samples $(x^{(1)}_1,x^{(2)}_1,\ldots,x^{(k)}_1)$, $(x^{(1)}_2,x^{(2)}_2,\ldots,x^{(k)}_2)$, …, $(x^{(1)}_n,x^{(2)}_n,\ldots,x^{(k)}_n)$ from $(X^{(1)},X^{(2)},\ldots,X^{(k)})$, and we have chosen the set $B \subseteq \R^k$ we are interested in, then we approximate \[\P\left( (X^{(1)},X^{(2)},\ldots,X^{(k)}) \in B \right) \approx \frac{\text{Nr of these } n \text{ samples that end up in $B$}}{n}. \tag{1.7}\]

3. Expectations of random vectors. The natural analogue of the expectation of a single random variable for a vector of random variables is simply a vector filled with expected values of the component random variables. However this concept is of no particular use to us, so we won’t further discuss this.

4. A function of a random vector. Another common application involving a random vector is when you’re interested in a random variable $X$ that is in some way made up from say $k$ other random variables/sources of randomness: \[X=f \left( X^{(1)},X^{(2)},\ldots,X^{(k)} \right)\] for some (known) function $f: \R^k \to \R$. In such cases Equation 1.5 and Equation 1.6 from Definition 1.1 can just directly be applied, with the observation that a sample $x$ from $X$ can be generated by producing a sample $(x^{(1)},x^{(2)},\ldots,x^{(k)})$ from the random vector $(X^{(1)},X^{(2)},\ldots,X^{(k)})$ and then setting \[x=f \left( x^{(1)},x^{(2)},\ldots,x^{(k)} \right).\]

Example 1.1 A real life example of point 4 in Remark 1.1: imagine that you are running some business, say an investment business where you invest on the financial markets on behalf of your clients (but almost any other type of business would do as well!). To do some proper risk analysis, you want to investigate your profit $X$ after say a year. Naturally your profit would depend on many uncertain factors: investment returns in a number of different markets/product categories, number of new clients, how much money these new clients bring in, etc. If we model each of these factors as a random variable (and decide on a distribution for them, typically by using past data and some suitable technique from Stats), then the profit is indeed some function that we can work out, say $f$, of these factors/random variables.

Note that a sample from all of these random variables, leading to a certain profit value by plugging them into our $f$, represents a possible way in which the coming year could unfold, i.e. a possible scenario if you will. By creating many such samples/possible scenarios, you basically create a “distribution” of many possible profit values, which allows you to assess e.g. what (approximately) your expected profit is, but also how likely it (approximately) is that you end the year with a negative profit i.e. a loss etc. — questions at the very heart of financial risk management!

1.3 Justification and error in the approximation

Before we further look at how to implement and use Monte Carlo simulations, let’s discuss the mathematical justification for the approximations in Definition 1.1 as well as whether we can come up with any error bounds. For simplicity it is enough to consider the approximation \[\E[X] \approx \frac{1}{n} \sum_{i=1}^n x_i \tag{1.8}\] where $X$ is some random variable and $x_1,\ldots,x_n$ are samples from it. Indeed the second formula in Definition 1.1 is just Equation 1.8, but applied to the random variable $Y:=g(X)$. For the first formula in Definition 1.1, the probability one, given some $B \subseteq \R$ we can just apply Equation 1.8 to the random variable $\mathbf{1}_A$, where $\mathbf{1}_A$ is the indicator function of the event $A=\{X \in B\}$ (I’ll leave the further details up to you, not so important for us).

As you may have guessed, the mathematical underpinning for Equation 1.8 consists of our lovely friends the Law of Large Numbers (LLN) and the Central Limit Theorem (CLT). You’ll probably remember these very well from your earlier Prob courses, but in case you’d like a refresher, have a look in Section A.7.

First for a bit of precise framing. Let $X_1, \ldots, X_n$ be a sequence of i.i.d. copies of $X$. Then obtaining $n$ samples $x_1,\ldots,x_n$ from $X$ is equivalent to obtaining a single sample $(x_1,\ldots,x_n)$ from the random vector $(X_1,\ldots,X_n)$, and in particular obtaining the number $\sum_{i=1}^n x_i$ is equivalent to obtaining a single sample from the random variable $S_n := \sum_{i=1}^n X_i$. So Equation 1.8 can equivalently be stated as follows: \[\E[X] \approx \text{a sample from the random variable } \frac{S_n}{n}. \tag{1.9}\] This formulation also stresses the following. You’ll be quite used to approximating some unknown number by some other number, you’ve done that zillions of times. But here we are rather approximating the number $\E[X]$ by a sample from a random variable! The crucial difference is that the right hand side is subject to uncertainty/randomness: if, in a given situation, we do the approximation Equation 1.9 once today and once tomorrow, then normally we will find two different samples and hence two different approximations of $\E[X]$! That’s a bit weird compared to those earlier approximations you’ve done so often, and very important to be well aware of. It’s not an actual problem — after all if we could somehow guarantee that any sample from $S_n/n$ is close enough (whatever meaning we want to give that exactly in our given situation) to $\E[X]$ then we would be happy to go forward with whatever sample we get. We won’t be able to get such a guarantee entirely, but we can get something close to it.

The first bit of justification for Equation 1.9 we immediately get from the LLN. Indeed it states that with probability $1$, we have that \[\lim_{n \to \infty} S_n/n = \E[X]. \tag{1.10}\] That is, if we would be able to produce a sample of infinite length from $X$ then, with probability $1$, Equation 1.9 would hold with an equality even². That’s very good news, because even though there is no such thing as an infinite length sample in the real world, it does strongly suggest that we are on the right track: namely, that for ‘large enough’ $n$, we should be able to get pretty close to $\E[X]$.

² The reason that we can’t say that Equation 1.10 happens with certainty but only the slightly weaker “with probability $1$” is due to unavoidable complications around uncountable sample spaces — something which you may or may not have/will encounter(ed) in a nice Prob course, not something we worry about in this course

So let’s investigate a little bit further what happens for ‘large’ $n$ exactly. For notational simplicity, let’s write \[\mu_X:=\E[X] \quad \text{and} \quad \sigma_X:=\sqrt{\var(X)}\] (assuming that $X$ has non-zero and finite variance) so that Equation 1.9 reads as \[\mu_X \approx \text{a sample from the random variable } \frac{S_n}{n}. \tag{1.11}\] Further we can easily compute, from the def of $S_n$ \[\begin{align*} \E \left[ \frac{S_n}{n} \right] &= \frac{1}{n} \sum_{i=1}^n \E[X_i] \\ &= \frac{1}{n} \sum_{i=1}^n \mu_X \\ &= \frac{1}{n} n \mu_X=\mu_X \end{align*} \] and \[\begin{align*} \var \left( \frac{S_n}{n} \right) &= \frac{1}{n^2} \sum_{i=1}^n \var(X_i) \\ &= \frac{1}{n^2} \sum_{i=1}^n \sigma_X^2 \\ &= \frac{1}{n^2} n \sigma_X^2 = \frac{\sigma_X^2}{n}. \end{align*} \] Recall that intuitively, the variance is a measure for how ‘spread out’ the values that a random variable takes are around its expected value. So we see that for large $n$, $S_n/n$ has quite a small variance and hence it tends to take values quite close to its expected value, which is exactly $\mu_X$. In other words, this intuition confirms that Equation 1.11 (and hence also Equation 1.9 and Equation 1.8) indeed holds!

To express this hopeful conclusion in a more explicit/measurable way, we have two options. Firstly we could turn to Chebyshev’s inequality (cf. Section A.8) which reads (with the choice $Z=S_n/n$ and plugging in the above expressions for the expectation and variance of $S_n/n$) as follows: for any $z>0$ we have that \[\P \left( \left\lvert \frac{S_n}{n} -\mu_X \right\rvert \geq z \frac{\sigma_X^2}{n} \right) \leq \frac{1}{z^2}\] or equivalently for $z \geq 1$ \[\P \left( \left\lvert \frac{S_n}{n} -\mu_X \right\rvert < z \frac{\sigma_X^2}{n} \right) \geq 1-\frac{1}{z^2}.\] That is, with a probability of at least $1-1/z^2$, the distance between a sample of $S_n/n$ and $\mu_X$ is less than $z \sigma_X^2/n$. Or, slightly reworded: if we produce a sample from $S_n/n$ and use that as our estimate for the unknown $\mu_X$, then with probability $1-1/z^2$ we make an error of \[z \sigma_X^2/n \tag{1.12}\] at most! If this doesn’t yet look very impressive to you, plug in a value for $z$, say $z=\sqrt{10}$ to get: if we produce a sample from $S_n/n$ and use that as our estimate for the unknown $\mu_X$, then with probability $0.9$ we make an error of at most $\sqrt{10} \sigma_X^2/n$. Here is our error bound in the approximation Equation 1.11 (and hence also Equation 1.9 and Equation 1.8)! Further, for any probability $p \in (0,1)$ (by choosing $z$ solving $1-1/z^2=p$), we can make the error bound Equation 1.12 as small as we want if only we make $n$, the number of samples that we’re using, large enough!

Remark 1.2.

Note how the error bound Equation 1.12 comes with the qualifier ‘with probability $p$’ for some $p \in (0,1)$ (depending on what value for $z$ you choose). This implies that there is also a positive probability, namely $1-p$, that a sample from $S_n/n$ does not respect that error bound i.e. is further away from $\mu_X$, in which case that sample makes for a relatively poor or even wildly wrong approximation of $\mu_X$! Unfortunately this is an intrinsic weakness of the Monte Carlo simulation approach.

And this is also not a real surprise actually, if you think about it. If we roll a fair die $n$ times and take the average outcome as estimate for the exact expected value of $3.5$, then it is possible that we get quite unlucky and roll only ⚀s for instance, leading to an estimate of $1$ — quite an error indeed! The larger we make $n$, the more unlikely such an unfavourable/unrepresentative sequence of samples is to occur, but it is never entirely ruled out. In practice that means that a Monte Carlo based estimate should always be approached with this possibility in mind, and requires sanity checking.
You could very rightly argue that in the above error bound $z \sigma_X^2/n$, we may very well not know what $\sigma_X$ is (after all, it may well be the case that we don’t know anything about $X$ other than how to produce samples from it) and hence in such a situation we can’t actually compute it. True! However, there aren’t too many situations in which working out that bound explicitly is very useful. If it is needed for some reason, then your best bet is typically to estimate it using the sample variance that you’ll remember from a previous Stats course.

We promised above that we would have two options to measure the error, indeed the second one uses that $S_n$ has a specific form: a sum of i.i.d. random variables. This means that we can look at our good old friend the CLT (recall from Section A.7) to analyse the distribution of $S_n$. Let’s fix some $h>0$ and consider the probability that in Equation 1.11 we make an error of at most $h$ i.e. \[\P \left( \frac{S_n}{n} \in [\mu_X-h,\mu_X+h] \right).\] Doing the usual standardisation \[\begin{align*} \frac{S_n}{n} \in [\mu_X-h,\mu_X+h] &\iff \frac{S_n-n\mu_X}{n} \in [-h,h] \\ &\iff \frac{S_n-n\mu_X}{\sqrt{n} \sigma_X} \in \left[ -h \frac{\sqrt{n}}{\sigma_X},h \frac{\sqrt{n}}{\sigma_X} \right] \end{align*} \] we can apply the usual CLT approximation (cf. Equation A.63) to derive that for ‘large’ $n$ \[\begin{align*} \P \left( \frac{S_n}{n} \in [\mu_X-h,\mu_X+h] \right) &= \P \left( \frac{S_n-n\mu_X}{\sqrt{n} \sigma_X} \in \left[ -h \frac{\sqrt{n}}{\sigma_X},h \frac{\sqrt{n}}{\sigma_X} \right] \right) \\ &\approx \Phi \left( h \frac{\sqrt{n}}{\sigma_X} \right) - \Phi \left( -h \frac{\sqrt{n}}{\sigma_X} \right) \\ &= 2 \Phi \left( h \frac{\sqrt{n}}{\sigma_X} \right) -1, \end{align*} \] where $\Phi$ is the cdf of a $\mathcal{N}(0,1)$ distribution. I.e. if we produce a sample from $S_n/n$ for ‘large’ $n$ and use that as our estimate for the unknown $\mu_X$, then with approximately probability $2 \Phi(h \sqrt{n}/\sigma_X)-1$ we make an error of at most $h$.

The advantage of the CLT approximation over Chebyshev’s inequality we used above is that it typically gives sharper bounds, but obviously the downside is that it only gives an approximate probability of respecting that bound and that it requires ‘large enough’ $n$ which is of course a bit vague³.

³ It is possible to make this more precise, see e.g. the discussion in Section A.7, but we leave that outside the course

Let’s put the key conclusions of this section together:

Proposition 1.1 Assume that $X$ has finite and non-zero variance. Then the Monte Carlo approximation Equation 1.8 (resp. Equation 1.9 and Equation 1.11) is generally a good one (though keep the first bullet point in Remark 1.2 in mind). Error bounds can be obtained using Chebyshev’s inequality: for $z \geq 1$ it holds that \[\P \left( \left\lvert \frac{S_n}{n} -\mu_X \right\rvert < z \frac{\sigma_X^2}{n} \right) \geq 1-\frac{1}{z^2}\] and using the CLT: for any $h>0$ and ‘large’ $n$ it holds that \[\P \left( \left\lvert \frac{S_n}{n} -\mu_X \right\rvert \leq h \right) \approx 2 \Phi \left( h \frac{\sqrt{n}}{\sigma_X} \right) -1,\] where $\Phi$ is the cdf of a $\mathcal{N}(0,1)$ distribution.

1.4 Some experiments in R

Finally in this chapter, let’s take a first look at the practical side of things: how to work with the Monte Carlo simulation technique from Definition 1.1 in the software package R.

1.4.1 Computer “randomness”

So in order to implement the Monte Carlo method, we’ll need to have our computer produce samples from a certain given random variable — essentially, execute a random experiment. However computers (well at least the non-quantum ones!) are fundamentally deterministic machines: they take some input and execute a series of actions as prescribed by whatever code it is running to produce some output. Nothing random there, all strictly following a set of rules.

But then here is a quite smart insight: it turns out that there exist carefully crafted deterministic algorithms that are able to generate real numbers⁴ from the interval $(0,1)$ that look and feel an awful lot like genuine samples from a $\text{Unif}(0,1)$ distribution: as they are generated there are no apparent patterns or connections between them (emulating independence) and it looks like they are at random scattered around the interval, gradually filling it up in a uniform way. Good enough for us! These are called pseudorandom numbers and the algorithm/piece of code producing them a pseudorandom number generator (PRNG). If you’re interested, a prominent class of PRNGs are the Mersenne Twister variations.

⁴ Or, more precisely, the computer implementation of those

So that’s a good first step, but of course we’re interested in samples from many other distributions as well, not just from the $\text{Unif}(0,1)$ distribution. That’s also a topic where mathematicians have worked hard on in recent decades: coming up with algorithms/“recipes” to turn samples from a $\text{Unif}(0,1)$ distribution into samples from some other distribution. We’ll discuss some of these in more detail in Chapter 2. The good news for us is that any half decent scientific/mathematical software package, including R, has a good number of such algorithms built in so that it can produce (pseudorandom) samples from a lot of “well known” distributions without further ado at our end, all we need to do is use the correct function call in that software.

1.4.2 To seed or not to seed

Typically a PRNG as introduced in Section 1.4.1 to produce (pseudorandom) samples from some distribution needs an input to start the algorithm with, called the seed (normally some integer value). This gives you two options:

Specify a seed value explicitly in your code. If you do so, then each time you run your code, the algorithm will generate the same sequence of samples (after all, under the hood it is a deterministic algorithm, so same input=same output!). This can be very useful for instance if you’re debugging your code: if you make a small change to your code and run it again, then with a fixed seed you know for sure that the change in output you get to see is actually due to your code changes and not due to a different sequence of samples being generated!
Do not specify a seed in your code. In that case, the algorithm still produces a seed for itself, and it tries to do that in a “random” way thereby also using the system time (in milliseconds) from your computer, so that in different runs of your code normally a different seed is generated and (hence) also a different sequence of samples.

You can set a seed value by using the set.seed() function. If you have set a seed but later in your code want to get rid of it (so go back to the situation in which no seed is specified) you can call set.seed(NULL). This can also be handy at the start of your code in case you set a seed later in your code, because if you set a seed then this is stored and preserved if you do a second run of your code.

Note that setting a seed “resets” the algorithm, in the sense that if you call set.seed() twice in your code with the same seed value, then after the second call the algorithm will restart with the same sequence of samples you got after the first call. Below code uses the sample() function to randomly pick 3 elements from the sequence $1,2,\ldots,25$ (which also uses the PRNG), just as a trivial example — this function will be discussed further in Section 1.4.3 below.

Listing 1.1: Playing with seeds

1.4.3 General discrete distributions (with finite range)

Firstly a word about discrete distributions with finite range. It is pretty straightforward to produce samples from such guys using the R function sample (R doc). Its first argument is a vector containing the possible values of your random variable, and the prob argument a vector of equal length giving the probabilities for each of the possible values.

For instance in the coin example we started this chapter with, in order to simulate ten tosses with a fair coin you could run the code

Listing 1.2: Simulating tosses of a fair coin (without a seed)

(obviously, since we didn’t set a seed the output will generally look differently each time you run this code).

Another example we’ve alluded to a few times is rolling a fair die. If we denote by $X$ the outcome of that experiment then, as we know well, $\E[X]=3.5$. Implementing the Monte Carlo simulation technique from Definition 1.1 to estimate $\E[X]$ could for instance be done as follows (recall the seq (R doc) and rep (R doc) functions for easily producing longer vectors):

Listing 1.3: Estimating the expected value for rolling a fair die

1.4.4 Built-in distributions

In addition to distributions with a finite range, those that are “built-in” into R (and have a large overlap with the well known ones i.e. the examples you have enocountered in your earlier Prob & Stats courses) are particularly simple to produce samples with. Indeed the ones appearing in Appendix B are all available: binom, geom, pois; beta, exp, gamma, norm, unif. For each of these distributions, four useful functions are available whose names are formed by one of the letters d, p, q and r followed by the distribution name:

d[name](x=.,...): the pmf or pdf evaluated at x,
p[name](q=.,...): the cdf evaluated at q,
q[name](p=.,...): the inverse cdf/quantile function evaluated at p,
r[name](n=.,...): returns a length n vector of samples.

The R doc gives a full list of all distributions available (by default). Note that the extra function arguments, in the above list represented by ..., are the parameter values that a distribution needs.

Remark 1.3. As a general rule, for a number of distributions there are (notoriously) different parametrisations in common use and so it is always important to double check in the documentation of the software package that you’re using whether their parametrisation is the same as the one you’re using!

Conveniently, for all distributions appearing in Appendix B R uses the same parametrisation as is used in that appendix, with one notable exception: for the Normal distribution, R uses as parameters the mean and the standard deviation (in that order) rather than the mean and the variance. So, if you e.g. want the value of the pdf of a $\mathcal{N}(5,4)$ evaluated at $1$, you use dnorm(1,5,2) rather than dnorm(1,5,4).

For example, suppose that $X=UV$, where $U$ and $V$ are independent random variables with $U \sim \mathcal{N}(1,4)$ and $V \sim \text{Beta}(1,1)$. This is a situation as in point 4 in Remark 1.1: we have the form \[X=f(X^{(1)},X^{(2)})\] with $X^{(1)}=U$, $X^{(2)}=V$, and $f(x^{(1)},x^{(2)})=x^{(1)} x^{(2)}$ to translate it into the notation used there. As these are both built-in distributions (and they are independent), we can simply use the rnorm resp. rbeta function to generate samples from $U$ resp. $V$, take their products to obtain samples from $X$, and then we could for instance approximate $\E[X]$ (with exact value $\E[U] \E[V]=1 \cdot 1/2=1/2$) and $\P(X>2)$ (less obvious!) using the formulas from Definition 1.1:

Listing 1.4: A slightly more interesting/useful example

Note that for estimating the probability in lines 9–18, we used sapply (R doc) to apply the function indicator to each element of the vector samples_X. The resulting vector samples_X_indicated contains (hence) zeros and ones, $0$ resp. $1$ if the corresponding element in samples_X was at most $2$ resp. larger than $2$. In order to use the probability formula/approximation in Definition 1.1, we should count how many elements in samples_X_indicated equal $1$ and divide that number by n. We use the little trick that in a vector of zeros and ones, the number of entries that equal $1$ is exactly equal to the sum of the elements.

Now, the above lines 9–17 could be written a lot more compactly, for instance be collapsed into the single line est_prob_X<-length(which(samples_X>2))/n, where the which function (R doc) is used to collect the indices of the elements in samples_X that are larger than $2$ in a new vector, and then the length function (R doc) gives (indeed) the length of that vector. If you’d prefer to use something like this then that’s fine of course. However the above approach using sapply is a more generally useful construction, and so good to be aware of!

As a final comment for this example: slightly more advanced applications of the structure we used in this example may well require that we apply a function of multiple variables to multiple vectors (so, the first vector gives the values for the first function argument, the second vector gives the values for the second function argument, etc.). Such a situation is nicely handled by the mapply function (R doc).

1.5 Some exercises

About the exercises

Each exercise has a (rough) indication of its difficulty, as follows:

*	easier: can be solved by (almost) only using relevant definitions/results,
**	medium: in addition to relevant definitions/results, needs a limited amount of work/creativity,
*******	harder: in addition to relevant definitions/results, needs a larger amount of work/serious creativity,
💀	warning: might make your brain hurt! These are mainly intended to provide some extra challenge for those of you keen on that and are generally quite hard. You don't need to worry about these too much for exam purposes.

The exam consists of mostly ** and *** level questions, some *, and possibly at most a few marks worth of 💀.

A bit of preaching: it is an incredibly important part of the study process to try and work on the exercises as much as possible. To become a better mathematician/learn new maths (and also to get a good exam mark ;)), above all you need to do it. And yes, of course that includes falling over things, and making mistakes, and getting stuck, and getting frustrated — all part of the game and what you’re supposed to be doing! Your lecturers have done that as well and still do it. What matters is that you don’t let that discourage you and that you make good use of the help and resources available to help you develop your skills. As part of that, many exercises have a hint in a block like this:

Hint (click to view)

Hint!

These are trying to help you on your way if you don’t know where to start or to provide some ideas if you get stuck. In spirit of the above, always have a look at these first and try again before you look at the full solution. (These hints are an extra service that won’t be available in the exam I’m afraid ;).)

Of course, we have our classes and there’s office hours, email etc. as well — I’m at any time very happy to help you with any questions you may have, and you should please never feel that any question is “too dumb” to ask!

Full/detailed solutions for the exercises will become available, just immediately below the exercises, immediately after our tutorial hour (you may have to refresh the page).

Exercise 1.1 [*] Let $X$ be the result of rolling a fair five-sided die (so the possible outcomes are $1,\ldots,5$). In R, simulate rolling this die $n$ times, where $n=10^k$ for $k=1,2,\ldots,6$, and each time use the samples generated to approximate $\E[X]$ as well as compute the relative error the approximation makes. How does the relative error change as $n$ grows?

Listing 1.5: You can do your coding here

Hint (click to view)

Using that $X$ is a simple discrete random variable, you shouldn’t have too many problems working out that $\E[X]=3$ (if you’re unsure, review Section A.1.2). For the R code, we need an only small variation of the R code in Listing 1.3 so have another look at that one if you’d like some ideas.

Full solution (click to view)

First note that since $X$ takes values $1,\ldots,5$ with probability $1/5$ each, we can easily compute the expectation of this discrete random variable as (cf. Equation A.10) \[\E[X]=1 \cdot \frac{1}{5}+2 \cdot \frac{1}{5}+3 \cdot \frac{1}{5}+4 \cdot \frac{1}{5}+5 \cdot \frac{1}{5}=3.\]

For the R code, since $X$ is a discrete random variable with finite range, we look at Section 1.4.3, in particular Listing 1.3. We can essentially copy that code, only adjusting the parameters in below line 5 accordingly. I’m using a loop to walk through the different values of $k$, but of course you could also just run the code multiple times for different values of $k$.

Listing 1.6: Solution code – click “Run Code” to see the ouput

Overall, the trend in your output should be that the larger $n$, the smaller the relative error (recall the discussion in Section 1.3) but you may very well also observe a case where a step to a larger $n$ results in a larger relative error! This is not some mistake but just a consequence of the probabilistic/random nature of the Monte Carlo approach, as also discussed in the first bullet point of Remark 1.2!

Exercise 1.2 [*/**] Consider two random variables: $X \sim \text{Gamma}(2.1,1)$ and $Y \sim \text{Gamma}(0.0021,0.001)$, so that $\E[X]=\E[Y]=2.1$. Use Monte Carlo simulations to approximate $\E[X]$ and $\E[Y]$ (separately). Run the code(s) a few times, with the same number of samples for both, keeping an eye on the error that is made. In which case is the approximation generally (much) better? Can you explain why?

Listing 1.7: You can do your coding here

Hint (click to view)

The R code should not be too complicated — recall from Section 1.4.4 that the Gamma distribution is one of the built-in ones, and so producing samples is a one-line for each random variable. With these samples available, we can estimate the expectation in two lines, see e.g. how that is done in Listing 1.3 if you’re not sure.

You should find that the approximation of $Y$ generally seems to be considerably less accurate. For an explanation, have a look at the error bound Equation 1.12 for instance, and wonder how that bound compares for $X$ vs for $Y$?

Full solution (click to view)

Writing the code is not very complicated: the Gamma distribution is one of the built-in ones (cf. Section 1.4.4) and hence producing samples is a one-liner for each, and then we use Equation 1.6 to approximate the expectation as we also did in Exercise 1.1 above. So for instance:

Listing 1.8: Solution code – click “Run Code” to see the ouput

If you run this code multiple times, you’ll see that the trend is that the approximation of $\E[Y]$ is (much) poorer than that of $\E[X]$. A sense for what is behind this are the error bounds we found in Section 1.3: they all involve the variance of the random variable we’re sampling from. Take for instance Equation 1.12: all else equal, a random variable with a larger variance results in a larger error bound. This suggests that a large variance results in a lower accuracy. That’s indeed what we see here: we can compute (using Appendix B as needed): \[\var(X)=\frac{2.1}{1.1^2} \approx 1.7 \quad \text{and} \quad \var(Y)=\frac{0.0021}{0.001^2}=2100.\] So indeed, from this perspective it is not surprising that the approximation of $\E[Y]$ performs worse.

Some intuition here: as discussed in Section 1.2, the key idea behind Monte Carlo simulations is that by producing many samples, we get a pretty good idea of the kind of values a rv tends to take and with what probabilities. The more “spread out” over the real line the values of a random variable tend to be (for which the variance is a key measure as you’ll remember), the more samples we need to properly explore that whole range. Or, equivalently: for a given number of samples, the information they provide us with about the distribution of the random variable becomes less reliable/accurate the more “spread out” the distribution is/the larger the variance is.

Exercise 1.3 [**] Let $X$ and $Y$ be independent, $\mathcal{N}(0,1)$ distributed random variables. Consider the pair $(X,Y)$ as a random point in $\R^2$, and let \[R=\sqrt{X^2+Y^2}\] be the (random) distance of that point to the origin. Using Monte Carlo simulations, approximate $\E[R]$.

Hint: recall the mapply function we mentioned in Section 1.4.4.

Listing 1.9: You can do your coding here

Hint (click to view)

Observe that we are dealing with a random variable $R$ here which is defined as a function, say $f$, of the pair of random variables $(X,Y)$. For applying Equation 1.6 to $R$ we need to produce samples $r_1, \ldots, r_n$ from $R$, and (following the discussion in Remark 1.1) we can do this by producing samples $(x_1,y_1)$, $(x_2,y_2)$, …, $(x_n,y_n)$ from $(X,Y)$ and then setting $r_1=f(x_1,y_1)$, $r_2=f(x_2,y_2)$, …, $r_n=f(x_n,y_n)$. For this final step, consider using the mapply function as we are dealing with a function of two variables here (as also mentioned in Section 1.4.4).

Full solution (click to view)

Note that we are dealing with a pair of random variables here rather than a single one. Recall that we discussed this in general in Remark 1.1. Following that remark, we use Equation 1.6 in the following form: \[\E[R]=\E[f(X,Y)] \approx \frac{1}{n} \sum_{i=1}^n f(x_i,y_i),\] where $f$ is the function $f(x,y)=\sqrt{x^2+y^2}$ and $(x_1,y_1)$, $(x_2,y_2)$, …, $(x_n,y_n)$ are samples from $(X,Y)$ i.e. $r_1=f(x_1,y_1)$, $r_2=f(x_2,y_2)$, …, $r_n=f(x_n,y_n)$ are samples from $R$. Since $X$ and $Y$ are independent, we can produce such samples simply by creating a sample from $X$, one from $Y$, put them together as a pair, and feed the pair to $f$.

For the implementation, we can do essentially the same thing as in the first few lines of Listing 1.4, except that line 4 in that code needs adjusting: rather than taking the product of samples we now need to put them through the function $f$ defined above. A general strategy for that is: implement the relevant $f$ in R, and then make use of mapply. The R function sapply used in Listing 1.4 can only handle applying a function of one variable to a vector of numbers, while here we want to apply a function of two variables to two vectors of numbers (one providing the first argument and one the second argument). That’s why we turn to mapply as the hint suggests.

This brings us to the following:

Listing 1.10: Solution code – click “Run Code” to see the ouput

Exercise 1.4 [**] Imagine an experiment where we scatter points uniformly throughout the square \[S=\{ (x,y) \, | \, x \in [-1,1] \text{ and } y \in [-1,1] \} \subseteq \R^2.\] Let $C$ be the unit circle in $\R^2$: \[C=\{ (x,y) \, | \, x^2+y^2 \leq 1 \} \subseteq \R^2.\] Due to the uniform nature, the probability that a point ends up in $C$ equals the area of $C$ divided by the area of $S$ i.e. $\pi/4$. Or: $\pi$ equals $4$ times the probability that a point ends up in $C$.

That means that by approximating the probability that a point ends up in $C$, we can obtain an approximation of $\pi$. Now approximate this probability using Monte Carlo simulations, and check what approximation of $\pi$ you get. Note that to scatter points uniformly troughout $S$, it is natural to use the random vector $(U,V)$, where $U$ and $V$ are independent, $\text{Unif}(-1,1)$ distributed random variables.

Listing 1.11: You can do your coding here

Hint (click to view)

Observe that the probability we need to approximate is $\P((U,V) \in C)$, i.e. a probability involving a random vector rather than (just) a random variable, and recall that we discussed this in general in Remark 1.1. So in particular you want to use Equation 1.7: the approximation is given by the number of samples of $(U,V)$ that end up in $C$ divided by the total number of samples. To implement this, consider following the same structure as in Listing 1.4. Use runif twice, once to produce samples from $U$ and once to produce samples from $V$ (they are independent random variables after all). Then write a function, analogue to indicator in Listing 1.4, that returns value $1$ if a sample $(u,v)$ is in $C$ and $0$ otherwise. A difference with Listing 1.4 is that this function will take two arguments, but that’s no problem if we use mapply rather than sapply as we also did in Exercise 1.3 above.

Full solution (click to view)

The probability we’re interested in is $\P((U,V) \in C)$, and we follow Equation 1.7 to approximate it given $n$ samples from $(U,V)$: \[\P((U,V) \in C) \approx \frac{\text{Nr of samples in } C}{n}.\] Here is one possible implementation. This code follows the same structure as in Listing 1.4: the logic of using an indicator function to test which points/samples are inside the unit circle and then compute the approximate probability is the same. The key difference is that now we want to apply that indicator function to a sample $(u,v)$ from $(U,V)$ rather than to a real valued sample from a single random variable, so rather than using sapply we turn to mapply again as we also did in Exercise 1.3 above.

Listing 1.12: Solution code – click “Run Code” to see the ouput

Exercise 1.5 [**] In the English Premier League there are 20 football teams, so each team plays $2 \cdot 19=38$ games and can earn at most $3 \cdot 38=114$ points. Suppose that the number of points that Manchester United resp. Manchester City earn during a season can be modelled as $X \sim \text{Binomial}(114,0.75)$ resp. $Y \sim \text{Binomial}(114,0.8)$, and we assume independence. Using Monte Carlo simulatons, find the approximate probability that United earns more points than City.

Listing 1.13: You can do your coding here

Hint (click to view)

Note that the prob we’re after here is $\P(X>Y)$ i.e. $\P((X,Y) \in B)$ where \[B=\{(x,y) \, | \, x>y \} \subseteq \R^2.\] Follow the same ideas as in Exercise 1.4 and Exercise 1.3 above and you’ll be fine!

Full solution (click to view)

Again (as in Exercise 1.4 and Exercise 1.3 above) we’re dealing here with a probability involving more than one random variable, in this case the pair $(X,Y)$. We’re asked to approximate $\P(X>Y)$ i.e. $\P((X,Y) \in B)$ where \[B=\{(x,y) \, | \, x>y \} \subseteq \R^2.\] Following Equation 1.7 again, we approximate \[\P((X,Y) \in B) \approx \frac{\text{Nr of samples in } B}{n}.\] Following the same ideas as in Exercise 1.4 and Exercise 1.3 above, we can e.g. do this as follows:

Listing 1.14: Solution code – click “Run Code” to see the ouput

So, despite the fact that the model seems to indicate a relatively small difference between the two teams (United earns every available point with prob $0.75$ and City with prob $0.8$), apparently this small difference per point translates into a massive advantage for City over a whole season!

Exercise 1.6 [**/***] Your home made pizzas are very popular with your fellow students, and you’ve decided to start a little business selling them to make some extra pocket money. With your fancy maths skills you’ve built the following model (all amounts are in £):

The number of pizzas you’re selling during a year can be modelled as a $\text{Poisson}(100)$ distribution.
Your cost per pizza is random as well, depending on price fluctuations of your ingredients etc. You model it as a $\mathcal{N}(8,1.5)$ distribution.
You sell your pizzas for £10 each.

Using Monte Carlo simulations, find approximations for the following:

The probability that a pizza is sold for a (positive) profit. Also compute the exact probability (using R).
The probability that during a given year, you make a (positive) profit. Does the result make sense to you?
Your expected profit per year.

Hint: for parts ii & iii, you may find the rep (R doc) function and a for loop (R doc) helpful to use.

Listing 1.15: You can do your coding here

Hint (click to view)

Part i is a pretty straightforward application of Equation 1.5 in Definition 1.1 and should be not much of a problem if you could do the previous questions.

Parts ii & iii are a bit trickier! The key question is: how do we produce samples from the profit? At the end of the day, this is nothing but (yet another) application of point 4 in Remark 1.1, though a bit more challenging than in the previous questions.

If you’re not sure where to start, consider seeing a sample of the profit as the result of a three step process:

Generate a sample $m$ from $M \sim \text{Poisson}(100)$,
Generate samples $x_1, \ldots, x_m$ from $X \sim \mathcal{N}(8,1.5)$,
Compute the resulting sample of the profit.

and implement this in R.

Full solution (click to view)

For part i, if we denote by $X \sim \mathcal{N}(8,1.5)$ the random cost for a pizza, then the probability we’re asked to look at is $\P(X<10)$. Approximating this using Monte Carlo simulations is a direct application of Equation 1.5 in Definition 1.1. To compute the eaxct probability, note that \[\P(X<10)=\P(X \leq 10)\] is nothing but the cdf of $X$ evaluated in $10$ (cf. Equation A.4) and we know from Section 1.4.4 how to compute this in R. So we get (for instance) — don’t forget about Remark 1.3!:

Listing 1.16: Solution code – click “Run Code” to see the ouput

Obviously, we would expect both the exact and approximate probability to be very close!

For parts ii & iii, note that if you sell $m$ pizzas with cost per pizza $x_1, \ldots, x_m$, then your profit $P$ is given by \[P=\sum_{k=1}^m 10-x_k. \tag{1.13}\] or by \[P=10m - \sum_{k=1}^m x_k. \tag{1.14}\] So to generate a sample from the (random) profit $P$, we could follow the following steps:

Generate a sample $m$ from $M \sim \text{Poisson}(100)$,
Generate samples $x_1, \ldots, x_m$ from $X \sim \mathcal{N}(8,1.5)$,
Compute the resulting sample of the profit using Equation 1.13 or Equation 1.14.

Once we have these samples of the profit, the required probability and expectation can be done in the (by now) usual way. Note that this is in principle nothing but the “function of a random vector” idea from point 4 in Remark 1.1 — just a bit more advanced than we have seen so far!

In the R code below, we make use of the rep (R doc) function and a for loop (R doc).

Listing 1.17: Solution code – click “Run Code” to see the ouput

Of course the results you’ll find are random, but very likely you’ll find that the estimated probability equals $1$. Hmm… This does (of course) not mean that the exact probability also equals $1$ — it is definitely not (essentially) guaranteed that you will make a positive profit. It’s just that the exact probablity is very close to $1$.