# Programming with stochastics: simulating random experiments

Let us look briefly at how to simulate samples from random experiments in Python.

Python constructs used: [Modules/import](https://docs.python.org/3/tutorial/modules.html), [Functions/def](https://docs.python.org/3/tutorial/controlflow.html#defining-functions), [Length/len](https://docs.python.org/3/library/functions.html#len), [Sets/set](https://docs.python.org/3/tutorial/datastructures.html#sets), [Generator expression](https://docs.python.org/3/tutorial/classes.html#generator-expressions), [Range](https://docs.python.org/3/tutorial/controlflow.html#the-range-function), [Assertion/assert](https://docs.python.org/3/reference/simple_stmts.html#the-assert-statement), [Counting/Counter](https://docs.python.org/3/library/collections.html#collections.Counter)

## Implementing a probability measure

```python
from fractions import Fraction

def P(A, Omega):
    """The probability of the event A,
       given equally likely outcomes from a sample space Ω."""
    return Fraction(len(A & Omega), len(Omega))
```

Now we use it to program a throw of a die:

```python
Omega = {1, 2, 3, 4, 5, 6}
even = set((x*2 for x in range(1,4)))
print("Even numbers of pips:",even)

print("Probability of an even number of pips:", P(even, Omega))
```

That gives us the probability measure. To draw a sample, though, we still need to get the "randomness" from somewhere.

## Importing randomness

Python's `random` module is an interface to pseudo-random numbers. The function `random.random` is a [Mersenne Twister](https://en.wikipedia.org/wiki/Mersenne_Twister) implemented directly in C. If you need "real" random numbers, for cryptography say, there is the `secrets` module. The seed of the Mersenne Twister can be set, and should be, so that random simulations are reproducible.

```python
from random import random as r
help(r)
print("1", r())
print("2", r())

import random
random.seed(1)
very_random = r()
print("3", very_random)
assert very_random == 0.13436424411240122
print("4", r())
random.seed(1)
print("5", r())
```

![XKCD 221: Random Number](images/random_number.png "RFC 1149.5 specifies 4 as the standard IEEE-vetted random number.")

Randall Munroe, [xkcd 221 "Random Number"](https://xkcd.com/221/) (CC BY-NC 2.5)

From a (pseudo-)random number between $0$ and $1$ (note: possibly $0$, but never $1$) we can generate random throws of a die:

```python
from math import floor
def transform_unit_to_dice(x):
    return floor(1 + 6*x)

assert list(range(1,7)) == [transform_unit_to_dice((x-1)/6)
                            for x in range(1,7)]

N = 100
throws = [transform_unit_to_dice(r()) for n in range(N)]
print(throws)
```

Let us count how often each number of pips came up. The *relative frequency* is the count divided by $N$:

```python
from collections import Counter

frequencies = Counter(throws)
for pips in range(1, 7):
    print(pips, "came up", frequencies[pips], "times, relative frequency", frequencies[pips] / N)
```

To spare you programming such transformations over and over, there are also `random.randint(1,6)`, `random.choice(range(1,7))` and `random.randrange(1,7)`. Here `randint` is shorthand for the corresponding `randrange`, and `choice` is somewhat more general.

But let us put on record what we have done: starting from a uniformly distributed "random variable" $X=$`random.random` with values in $[0,1)$, we constructed and implemented a map $t =$`transform_unit_to_dice` with values in $\Omega := \{1,2,3,4,5,6\}$, and $t(X)$ is uniformly distributed. The mathematical map $t \colon [0,1) \to \Omega$ is a random variable, and we treat the [composition](https://en.wikipedia.org/wiki/Function_composition) $t \circ X$ as the random variable that models the die.

Now you might object that `random.random()` takes no parameter at all, so it is not a mathematical map from a domain into the set $[0,1)$. In fact we have to imagine it as a map $X \colon \Omega' \to [0,1)$, with some probability measure or other defined on $\Omega'$ in such a way that $X$ induces the uniform distribution on $[0,1)$. The set $\Omega'$ plays no concrete role for us – it is where "the randomness comes from", and the notation `random.random()` already tells us that we are not plugging in any particular element of $\Omega'$, but drawing one pseudo-randomly and plugging that into $X$.

## Larger samples

The function `random.sample(population, k)` draws a sample of size $k$ from a population (an urn) **without** replacement, so the elements drawn are pairwise distinct. For $k=1$ this amounts to a uniform distribution on the population.
The function `random.choices(population, weights=None, *, cum_weights=None, k=1)` draws from the population $k$ times **with** replacement (and lets you assign weights instead of assembling a suitable population).

## Looking ahead

The relative frequencies we counted in this notebook for the numbers of pips $1,\dots,6$ are estimates of the die's probabilities $p(1),\dots,p(6)$. In the next chapter we will call the function $p$ the *probability mass function* (pmf) of the random variable 'number of pips', and make the throw of a die our first statistical model.