# Coursework 2: probabilistic ranking
_[Original](http://mlg.eng.cam.ac.uk/teaching/4f13/1718/cw/coursework2.pdf) by Carl Rasmussen and Manon Kok for [CUED course 4f13](http://mlg.eng.cam.ac.uk/teaching/4f13/1718/). This version adapted by Damon Wischik._

In this assignment, you'll be using the (binary) results of the 2011 ATP men's tennis singles for 107 players in a total of 1801 games (which these players played against each other in the 2011 season), to compute probabilistic rankings of the skills of these players.

**What to submit.**
Your answers should contain an explanation of what you do, and
2&ndash;4 central commands to achieve it. Complete listings are
unnecessary. The focus of your answer should be
_interpretation:_ explain what the numerical values and graphs
you produce _mean,_ and why they are as they are.  The text of
your answer to each question should be no more than a paragraph or
two.


```python
import numpy as np
import pandas
import matplotlib
import matplotlib.pyplot as plt
import scipy.io, scipy.stats, scipy.linalg
from numpy.linalg import solve
import requests, io
from sympy.functions.special.delta_functions import DiracDelta
%matplotlib inline
import time
```

## Data import

The match data is provided as `http://mlg.eng.cam.ac.uk/teaching/4f13/1920/cw/tennis_data.mat`. It contains 
a vector $W$ of length 107 whose $i$th entry is the name of player $i$; and 
an array $G$ of dimension $1801\times2$, one row per game, where the first column contains the identity of 
the player who won and the second column contains the identify of the player who lost. Note that this convention means that the variable $y_g$ (the game outcome) in lecture notes is always $+1$, and can consequently be ignored. Some rows will appear more than once (corresponding to two players having played each other several times with the same outcome).


```python
r = requests.get('http://mlg.eng.cam.ac.uk/teaching/4f13/1920/cw/tennis_data.mat')
with io.BytesIO(r.content) as f:
    data = scipy.io.loadmat(f)
    W = np.concatenate(data['W'].squeeze())
    G = data['G'] - 1   # W[G[i,0]] is winner of game i, W[G[i,1]] is loser
    M = W.shape[0]      # number of players M = 107
    N = G.shape[0]      # number of games N = 1801


```

The appendix at the bottom of this document shows a crude ranking of the players, based on the fraction of their games that they win.

## Question (a)

Complete the following code for Gibbs sampling, by adding the lines required to sample from the conditional distributions needed for Gibbs sampling for the ranking model discussed in lectures. Run the Gibbs sampler, e.g. for 1100 iterations. Plot some of the sampled player skills as a function of the Gibbs iteration (see the [appendix](#Appendix) for some tips on multiple plotting).

Does it look as though the Gibbs sampler is able to move around the posterior distribution? What are the burn-in and mixing times for the Gibbs sampler? It may be helpful to look at autocorrelation plots; see the [appendix](#Appendix).


```python
pv = 0.5 * np.ones(M)   # prior variance
w = np.zeros(M)         # skills, initialized to be the prior mean μ0 = 0
    
Σinv_1 = np.zeros((M,M),dtype='float')


for g in range(N):
    Σinv_1[G[g][0],G[g][0]] += 1
    Σinv_1[G[g][1],G[g][1]] += 1
    Σinv_1[G[g][1],G[g][0]] -= 1
    Σinv_1[G[g][0],G[g][1]] -= 1

Σinv_0 = np.diag(np.ones((M,), dtype=int)*(1/0.5))

Σinv = Σinv_0 + Σinv_1

w_1 = np.array([])
w_2 = np.array([])
w_3 = np.array([])
iterations = np.array([])


for _ in range(1100):
    # Sample performance differences (t) given skills (w) and outcomes (G)
    s = w[G[:,0]] - w[G[:,1]]
    σ = 1
    t = s + σ * scipy.stats.norm.ppf(1 - np.random.uniform(size=N)*(1-scipy.stats.norm.cdf(-s/σ)))
            
    # Sample skills given performance differences
    Σ = np.linalg.inv(Σinv)
    μtilde = np.zeros((M,),dtype='float')
    for g in range(N):
        μtilde[G[g][0]] += t[g]
        μtilde[G[g][1]] -= t[g]

    μ = Σ @ np.transpose(μtilde)
    w = np.random.multivariate_normal(mean=μ, cov=Σ)
    
    if _ % 10 == 0:
        w_1 = np.append(w_1, [w[0]])
        w_2 = np.append(w_2, [w[1]])
        w_3 = np.append(w_3, [w[2]])
        iterations = np.append(iterations, [_])


fig = plt.figure(figsize = [10, 10])
plt.subplots_adjust(hspace = 0.5)
ax = fig.add_subplot(3, 1, 1)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_1, linewidth=4)
ax.set_title(W[0], fontsize = 16)
ax.grid()
ax.set_ylim([-0.5, 2.5])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)
ax = fig.add_subplot(3, 1, 2)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_2, linewidth=4)
ax.set_title(W[1], fontsize = 16)
ax.grid()
ax.set_ylim([-0.5, 2.5])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)
ax = fig.add_subplot(3, 1, 3)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_3, linewidth=4)
ax.set_title(W[2], fontsize = 16)
ax.grid()
ax.set_ylim([-0.5, 2.5])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)

w_1 = np.cumsum(w_1)
w_2 = np.cumsum(w_2)
w_3 = np.cumsum(w_3)
for i in range(w_1.shape[0]):
    if i != 0:
        w_1[i] = w_1[i]/i 
        w_2[i] = w_2[i]/i 
        w_3[i] = w_3[i]/i 

fig = plt.figure(figsize = [10, 10])
plt.subplots_adjust(hspace = 0.5)
ax = fig.add_subplot(3, 1, 1)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_1, linewidth=4)
ax.set_title(W[0], fontsize = 16)
ax.grid()
#ax.set_ylim([-0.5, 2.5])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)
ax = fig.add_subplot(3, 1, 2)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_2, linewidth=4)
ax.set_title(W[1], fontsize = 16)
ax.grid()
#ax.set_ylim([-0.5, 0.5])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)
ax = fig.add_subplot(3, 1, 3)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_3, linewidth=4)
ax.set_title(W[2], fontsize = 16)
ax.grid()
#ax.set_ylim([-0.5, 2.5])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)

print(w)
```

## Question (b)

Do inference in the model instead, by running message passing and expectation propagation, using the code below. For the same players you plotted in part (a), show the mean and standard deviation of their skill levels as a function of the iteration. How long does it take to converge?

(See the Python programming note in the appendix, about the `yield` statement and how to use `gaussian_ep`.)


```python
def gaussian_ep(G, M):
    def Ψ(x): return scipy.stats.norm.pdf(x) / scipy.stats.norm.cdf(x)
    def Λ(x): return Ψ(x) * (Ψ(x) + x)
    N = len(G)

    μ_s, p_s = np.empty(M), np.empty(M)
    μ_gs, p_gs = np.zeros((N,2)), np.zeros((N,2))
    μ_sg, p_sg = np.empty((N,2)), np.empty((N,2))
    
    while True:
        # 1. Compute marginal skills
        # Let skills be N(μ_s, 1/p_s)
        p_s = np.ones(M) * 1/0.5
        μ_s = np.zeros(M)
        for j,(winner,loser) in enumerate(G):
            p_s[winner] += p_gs[j,0]
            p_s[loser] += p_gs[j,1]
            μ_s[winner] += μ_gs[j,0] * p_gs[j,0]
            μ_s[loser] += μ_gs[j,1] * p_gs[j,1]
        μ_s = μ_s / p_s

        # 2. Compute skill -> game messages
        # winner's skill -> game: N(μ_sg[,0], 1/p_sg[,0])
        # loser's skill -> game: N(μ_sg[,1], 1/p_sg[,1])
        p_sg = p_s[G] - p_gs
        μ_sg = (p_s[G]*μ_s[G] - p_gs*μ_gs) / p_sg

        # 3. Compute game -> performance messages
        v_gt = 1 + np.sum(1/p_sg, 1)
        σ_gt = np.sqrt(v_gt)
        μ_gt = μ_sg[:,0] - μ_sg[:,1]

        # 4. Approximate the marginal on performance differences
        μ_t = μ_gt + σ_gt * Ψ(μ_gt/σ_gt)
        p_t = 1 / v_gt / (1-Λ(μ_gt/σ_gt))

        # 5. Compute performance -> game messages
        p_tg = p_t - 1/v_gt
        μ_tg = (μ_t*p_t - μ_gt/v_gt) / p_tg

        # 6. Compute game -> skills messages
        # game -> winner's skill: N(μ_gs[,0], 1/p_gs[,0])
        # game -> loser's skill: N(μ_gs[,1], 1/p_gs[,1])
        p_gs[:,0] = 1 / (1 + 1/p_tg + 1/p_sg[:,1])  # winners
        p_gs[:,1] = 1 / (1 + 1/p_tg + 1/p_sg[:,0])  # losers
        μ_gs[:,0] = μ_sg[:,1] + μ_tg
        μ_gs[:,1] = μ_sg[:,0] - μ_tg
        
        yield (μ_s, np.sqrt(1/p_s))

w_1_mean = np.array([])
w_2_mean = np.array([])
w_3_mean = np.array([])
w_1_std = np.array([])
w_2_std = np.array([])
w_3_std = np.array([])
iterations = np.array([])

for i, element in enumerate(gaussian_ep(G, M)):
    if i > 10:
        if abs(w_1_mean[w_1_mean.shape[0]-1] - element[0][0]) < 0.001:
            break
    w_1_mean = np.append(w_1_mean, element[0][0])
    w_2_mean = np.append(w_2_mean, element[0][1])
    w_3_mean = np.append(w_3_mean, element[0][2])
    w_1_std = np.append(w_1_std, element[1][0])
    w_2_std = np.append(w_2_std, element[1][1])
    w_3_std = np.append(w_3_std, element[1][2])
    iterations = np.append(iterations, [i])

   

fig = plt.figure(figsize = [10, 10])
plt.subplots_adjust(hspace = 0.5)
ax = fig.add_subplot(3, 1, 1)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_1_mean, linewidth=4, label='Mean')
ax.fill_between(iterations.squeeze(), w_1_mean.squeeze()-w_1_std, w_1_mean.squeeze()+w_1_std, alpha=.2, label='Standard deviation')
ax.set_title(W[0], fontsize = 16)
ax.grid()
ax.set_ylim([-0.5, 2])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)
ax.legend(loc = 'lower right')
ax = fig.add_subplot(3, 1, 2)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_2_mean, linewidth=4, label='Mean')
ax.fill_between(iterations.squeeze(), w_2_mean.squeeze()-w_2_std, w_2_mean.squeeze()+w_2_std, alpha=.2, label = 'Standard deviation')
ax.set_title(W[1], fontsize = 16)
ax.set_ylim([-0.5, 2])
ax.grid()
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)
ax.legend(loc = 'best')
ax = fig.add_subplot(3, 1, 3)  # 2 rows of panels, 3 columns
ax.plot(iterations, w_3_mean, linewidth=4, label='Mean')
ax.fill_between(iterations.squeeze(), w_3_mean.squeeze()-w_3_std, w_3_mean.squeeze()+w_3_std, alpha=.2, label='Standard deviation')
ax.set_title(W[2], fontsize = 16)
ax.grid()
ax.set_ylim([-0.5, 2])
ax.set_xlabel("Iterations", fontsize = 14)
ax.set_ylabel("Player skill", fontsize = 14)
ax.legend(loc = 'lower right')



```

## Question (c)

Explain the concept of _convergence_ for both the Gibbs sampler and the message passing algorithms. What type of object are we converging to in the two cases, and how do you judge convergence? For each of the players you plotted in (a), find the mean and standard deviation of their skill as found by the converged algorithm in (b), and annotate the plot from part (a) to show those values. Explain what you see.

## Question (d)

Who is better out of `W[0]` Rafael Nadal and `W[15]` Novak Djokovic? There are several ways we could answer this.

* Using the Gibbs sampler (and taking account of burn-in and mixing time), first compute the marginal distribution of skill for each of these players, then work out the probability that Nadal's skill is higher.
* Using the Gibbs sampler, use joint samples from the Gibbs sampler to directly estimate the probability.
* Using the message passing algorithm, read off the distribution of skills of these two players, and work out the probability that Nadal's skill is higher.
* Repeat these three methods, but this time to work out the probability that Nadal would beat Djokovic in a match.

Calculate all these six probabilities. Explain how they relate to each other. Which answer is best? Why? Support your explanation with appropriate plots.


```python
pv = 0.5 * np.ones(M)   # prior variance
w = np.zeros(M)         # skills, initialized to be the prior mean μ0 = 0
    
Σinv_1 = np.zeros((M,M),dtype='float')


for g in range(N):
    Σinv_1[G[g][0],G[g][0]] += 1
    Σinv_1[G[g][1],G[g][1]] += 1
    Σinv_1[G[g][1],G[g][0]] -= 1
    Σinv_1[G[g][0],G[g][1]] -= 1

Σinv_0 = np.diag(np.ones((M,), dtype=int)*(1/0.5))

Σinv = Σinv_0 + Σinv_1

w_1 = np.array([])
w_2 = np.array([])
w_3 = np.array([])
iterations = np.array([])


for _ in range(1100):
    # Sample performance differences (t) given skills (w) and outcomes (G)
    s = w[G[:,0]] - w[G[:,1]]
    σ = 1
    t = s + σ * scipy.stats.norm.ppf(1 - np.random.uniform(size=N)*(1-scipy.stats.norm.cdf(-s/σ)))
            
    # Sample skills given performance differences
    Σ = np.linalg.inv(Σinv)
    μtilde = np.zeros((M,),dtype='float')
    for g in range(N):
        μtilde[G[g][0]] += t[g]
        μtilde[G[g][1]] -= t[g]

    μ = Σ @ np.transpose(μtilde)
    w = np.random.multivariate_normal(mean=μ, cov=Σ)

    if _ % 1 == 0:
        w_1 = np.append(w_1, [w[0]])
        w_2 = np.append(w_2, [w[15]])
        iterations = np.append(iterations, [_])

print(w_1[500:600])
print(w_2[500:600])
```

    [1.55487112 1.41998191 1.51073664 1.25296892 1.37212907 1.20484217
     1.06091822 1.20986535 1.3472423  1.55660573 1.5817697  1.90669306
     1.81189225 1.56300879 1.58256473 1.7232599  1.55447544 1.75191393
     1.77999671 1.27124817 1.41429748 1.42490659 1.48637662 1.54709257
     1.62029956 1.50602068 1.74333146 1.53385639 1.69202299 1.6482965
     1.40348136 1.4321042  1.19463614 1.39867811 1.161957   1.36712037
     1.63737991 1.65011926 1.59995507 1.68312856 1.68933975 1.61445527
     1.86172248 1.95199603 1.57069288 1.56157303 1.65445182 1.4762276
     1.58276378 1.76319391 1.99854725 2.05947297 1.58288166 1.49220072
     1.3394894  1.45482831 1.58251882 1.63997027 1.66523021 1.59100603
     1.64172296 1.60901921 1.44782576 1.70705828 1.71774222 1.3651469
     1.37836906 1.43912102 1.4828025  1.52328383 1.53726591 1.52450137
     1.84461729 1.54409358 1.46054353 1.21726762 1.26477711 1.49029226
     1.32992047 1.53519848 1.37633123 1.41419537 1.17851943 1.60745878
     1.43195279 1.72966086 1.31583328 1.53664918 1.51957318 1.49697634
     1.44670705 1.53994912 1.51736767 1.44240859 1.24326778 1.68196966
     1.3858646  1.53769087 1.34521931 1.3093682 ]
    [1.99175284 2.18171449 1.97628417 2.04680433 2.22804569 1.91735068
     1.76821901 1.84780697 2.01081663 2.18749236 2.12607603 2.37215767
     2.2414722  2.25881482 2.2900471  2.16208446 1.92105467 2.01203609
     1.92491431 1.5932097  1.5423015  1.81278844 1.80284134 1.67290648
     1.61609232 1.47196921 1.737638   1.56600453 1.66832421 1.96186392
     1.50100271 1.73808619 1.63263494 1.92262471 1.9782193  2.02673454
     2.43148896 2.09237299 2.07558058 2.13837829 2.24345468 1.90780855
     1.80748489 2.02237271 1.97287309 1.92543681 1.84639478 1.63760813
     1.68005753 1.70328329 2.047697   2.01511555 1.91578126 2.13266812
     2.30674951 2.37997954 2.20028298 2.32609989 2.24570691 2.46692475
     2.33866894 2.19752714 2.03934303 2.04086433 1.89998046 2.10306589
     2.21814003 1.96541133 1.82504005 1.90427884 1.89092214 1.84391452
     2.29690592 2.1859445  1.86425523 1.73348657 1.87092915 1.97796014
     1.96501385 2.08244451 1.922393   2.43681365 2.06187936 2.06214927
     1.99001657 2.23155935 1.96886502 1.99179602 2.37642703 2.23950582
     1.86832784 1.79150089 1.79928328 1.82539083 1.80314857 2.11607047
     1.64499637 1.65978288 1.72054028 1.54277598]


## Question (e)

Jarrett (Biometrika, 1979) compiled a dataset of the dates of coal mining explosions, from the _Colliery Year Book and Coal Trades Directory_. An initial plot suggests that there was a step change in the frequency of explosions, some time between 1880 and 1900. Consider the model
$$
\text{num.explosions in year }i \sim \operatorname{Poisson}(\lambda 1_{i\leq\theta} + \mu 1_{i>\theta})
$$
with prior distributions
$$
\lambda \sim \operatorname{Gamma}(a,b)\,,\quad
\mu \sim \operatorname{Gamma}(a,b)\,,\quad
\theta \sim \operatorname{Uniform}[1851,1852,\dots,1961]
$$
Use Gibbs sampling to estimate $\theta$. It's up to you to pick appropriate $a$ and $b$.

_This question asks you to invent a Gibbs sampler from scratch. This model is much simpler than the ranking model
studied in lectures. In your answer, you should explain how your Gibbs sampler works, but you shouldn't spend time worrying about burn-in and mixing time._


```python
url = 'https://raw.githubusercontent.com/damonjw/probml/master/data/coal.txt'
coal = pandas.read_csv(url, names=['date'], parse_dates=['date'])

# Count number of explosions by year, and fill in missing years with count=0
counts = coal['date'].groupby(coal['date'].dt.year).apply(len)
y,Y = np.min(counts.index.values), np.max(counts.index.values)
df = pandas.DataFrame({'count': np.zeros(Y-y+1, dtype=int)}, index=np.arange(y,Y+1))
df.index.name = 'year'
df.loc[counts.index, 'count'] = counts
# Discard first and last year (for which counts might be incomplete)
df = df.reset_index()
df = df.loc[(df['year']>y) & (df['year']<Y)]

with matplotlib.rc_context({'figure.figsize': [12,2]}):
    plt.bar(df['year'].values, df['count'].values)
    plt.title('mine explosions per year')
plt.show()

y = df['count'].values
years = np.arange(0,110,1)
last_year = years[years.shape[0]-1]
first_year = 0

a = np.mean(y)**2/np.var(y)
b = np.var(y)/np.mean(y)
a_lamda = a
b_lamda = 1/b
a_mu = a
b_mu = 1/b
p_theta = np.zeros(years.shape)

# Initialisation:
lamda = np.random.gamma(a, b)
mu  = np.random.gamma(a, b)
theta = np.random.choice(years)
old_theta = theta + 4


while abs(theta - old_theta) > 3:
    old_theta = theta
    # Sample lamda:
    a_lamda = sum(y[first_year:theta]) + 111*a_lamda -110
    b_lamda = 1/(b_lamda + theta - first_year)
    print(a_lamda)
    #print(b_lamda)
    lamda = np.random.gamma(a_lamda, b_lamda)
    
    # Sample mu:
    a_mu= sum(y[theta+1:last_year]) + a_mu
    b_mu = 1/(b_mu + last_year-theta)
    print(a_mu)
    #print(b_mu)
    mu = np.random.gamma(a_mu, b_mu)
    
    # Sample theta:
    for j in years:
        if j == 0:
            p_theta[j] = lamda**(sum(y))*np.exp((mu+lamda))
        else:
            p_theta[j] = mu**(sum(y[first_year:j-1]))*lamda**(sum(y[j:last_year]))*np.exp((mu+lamda))
    p_theta = p_theta/sum(p_theta)
    theta = np.random.multinomial(1, p_theta, size=1)
    for i, element in enumerate(theta[0]):
        if element != 0:
            theta = i
            break
                
    



```


```python
a = np.ones(4)
print(sum(a[0:1]))
```

    1.0


<div style="height: 20em"></div>

# Appendix

### Data import

Here is a crude ranking of the players, based on a simple ratio
$$
\frac{\text{number of wins for player }i}{\text{number of games played by player }i}.
$$


```python
wins = np.zeros(M, dtype=int)
losses = np.zeros(M, dtype=int)
w, x = np.unique(G[:,0], return_counts=True)
wins[w] = x
w, x = np.unique(G[:,1], return_counts=True)
losses[w] = x
score = wins / (wins + losses)
rank_order = np.argsort(score)[::-1]

with matplotlib.rc_context({'figure.figsize': [20,5], 'xtick.labelsize': 8}):
    x = np.arange(M)
    plt.bar(x, score[rank_order], align='center', width=.8)
    plt.xticks(x, W[rank_order], rotation=-90, ha='right')
    plt.ylabel('%win')
plt.show()
```

### Question (a)

Some notes about the code in Question (a), and how it uses normally distributed random variables:

* If we just wanted to sample `t` given skill differences `s`, we could just use `t=np.random.normal(loc=s, scale=σ)`. But we actually want to generate `t` conditional on `t>0`. Using the [inversion transformation method](https://en.wikipedia.org/wiki/Inverse_transform_sampling) as mentioned in the notes, we end up with the code given here.
* We could also compute `μ` by defining `Rinv = scipy.linalg.cholesky(Σinv)` and then `μ=solve(Rinv, solve(Rinv.T, μtilde))`.
* We could also generate `w` by `w = μ + solve(Rinv, np.random.normal(size=M))`
* To generate a simple normal random variable with mean $\mu$ and variance $\sigma^2$, we can first generate $X\sim\operatorname{Normal}(0,1)$ and then return $\mu+\sigma X$. To generate a multivariate normal random variable with mean vector $\mu$ and covariance matrix $\Sigma$, we first generate $X\sim\operatorname{MVN}(0,I)$ and then return $\mu+R X$ where $R$ is the matrix square root of $\Sigma$. It can be found from the Cholesky decomposition $R R^\top=\Sigma$.

To plot many panels in a single figure,


```python
with matplotlib.rc_context({'figure.figsize': [6,4], 'axes.titlesize': 10}):
    fig = plt.figure()
    for i in range(5):
        ax = fig.add_subplot(2, 3, i+1)  # 2 rows of panels, 3 columns
        ax.plot(np.random.normal(size=10), np.random.normal(size=10))
        ax.set_title(f'Panel {i}')
```

To plot autocorrelations,


```python
# Generate synthetic data with autocorrelations
x = [5, 5]
for _ in range(1000):
    lastx = x[-2:]
    nextx = np.mean(lastx) + 0.1 * np.random.normal() - 0.3 * (lastx[0] - 3)
    x.append(nextx)
x = np.array(x)
    
# Subtract the mean, then plot autocorrelation
plt.acorr(x - np.mean(x), maxlags=40)
plt.xlim(-1,40)
plt.show()
```

### Question (b)

The `gaussian_ep` function includes a `yield` statement, which means it is a lazy list generator, and you can use it with e.g.
```
g = gaussian_ep(G, M)
for _ in range(10):
    print(next(g))
```

### Question (d)

To evaluate $\mathbb{P}(\operatorname{Normal}(\mu,\sigma^2)>x)$, you can use the appropriate scientific function directly,
```
scipy.stats.norm.sf(x, loc=μ, scale=σ)
```
or you can compute it using Monte Carlo integration,
```
np.mean(np.random.normal(size=10000, loc=μ, scale=σ) > x)
```
For this simple calculation there's no point in Monte Carlo integration, but the technique
is useful in situations where there is no simple built-in scientific function.

fe63e8ff8cc39eab40bdb5303d5fda543f23751d
ipynb
Jupyter Notebook
coursework2.ipynb
jt2718/ProbMLCambridge
aae8c4f38dcf2b619f36bad7c1abd41177710f9e
MIT
coursework2.ipynb
jt2718/ProbMLCambridge
aae8c4f38dcf2b619f36bad7c1abd41177710f9e
MIT
coursework2.ipynb
jt2718/ProbMLCambridge
aae8c4f38dcf2b619f36bad7c1abd41177710f9e
MIT
Qwen/Qwen-72B
1. YES
2. YES