@@ -35,7 +35,7 @@ In the lecture we'll learn about

3535* using a Monte Carlo simulation of a multivariate normal distribution to evaluate the quality of a normal approximation3636* the administrator's problem and why the multivariate hypergeometric distribution is the right tool373738-## The Administrator's Problem38+## The administrator's problem39394040An administrator in charge of allocating research grants is in the following situation.4141

@@ -62,7 +62,7 @@ The $n$ balls drawn represent successful proposals and are awarded research fu

62626363The remaining $N-n$ balls receive no research funds.646465-### Details of the Awards Procedure Under Study65+### Details of the awards procedure under study66666767Let $k_i$ be the number of balls of color $i$ that are drawn.6868

@@ -106,11 +106,11 @@ the population of $N$ balls.

106106107107The right tool for the administrator's job is the **multivariate hypergeometric distribution**.108108109-### Multivariate Hypergeometric Distribution109+### Multivariate hypergeometric distribution110110111111Let's start with some imports.112112113-```{code-cell} ipython113+```{code-cell} ipython3114114import matplotlib.pyplot as plt115115import numpy as np116116from scipy.special import comb

@@ -159,7 +159,7 @@ $$

159159160160To do our work for us, we'll write an `Urn` class.161161162-```{code-cell} python3162+```{code-cell} ipython3163163class Urn:164164165165 def __init__(self, K_arr):

@@ -209,20 +209,14 @@ class Urn:

209209 number of draws.210210 """211211212- K_arr, N, c = self.K_arr, self.N, self.c212+ K_arr, N = self.K_arr, self.N213213214214 # mean215215 μ = n * K_arr / N216216217217 # variance-covariance matrix218- Σ = np.full((c, c), n * (N - n) / (N - 1) / N ** 2)219- for i in range(c-1):220- Σ[i, i] *= K_arr[i] * (N - K_arr[i])221- for j in range(i+1, c):222- Σ[i, j] *= - K_arr[i] * K_arr[j]223- Σ[j, i] = Σ[i, j]224-225- Σ[-1, -1] *= K_arr[-1] * (N - K_arr[-1])218+ p = K_arr / N219+ Σ = n * (N - n) / (N - 1) * (np.diag(p) - np.outer(p, p))226220227221 return μ, Σ228222

@@ -265,15 +259,15 @@ $$

265259

P(2{\text{ black}},2{\text{ white}},2{\text{ red}})={{{5 \choose 2}{10 \choose 2}{15 \choose 2}} \over {30 \choose 6}}=0.079575596816976

266260$$267261268-```{code-cell} python3262+```{code-cell} ipython3269263# construct the urn270264K_arr = [5, 10, 15]271265urn = Urn(K_arr)272266```273267274268

Now use the Urn Class method `pmf` to compute the probability of the outcome $X = \begin{bmatrix} 2 & 2 & 2 \end{bmatrix}$

275269276-```{code-cell} python3270+```{code-cell} ipython3277271k_arr = [2, 2, 2] # array of number of observed successes278272urn.pmf(k_arr)279273```

@@ -283,99 +277,99 @@ constructing a 2-dimensional

283277array `k_arr` and `pmf` will return an array of probabilities for284278observing each case.285279286-```{code-cell} python3280+```{code-cell} ipython3287281k_arr = [[2, 2, 2], [1, 3, 2]]288282urn.pmf(k_arr)289283```290284291285Now let's compute the mean vector and variance-covariance matrix.292286293-```{code-cell} python3287+```{code-cell} ipython3294288n = 6295289μ, Σ = urn.moments(n)296290```297291298-```{code-cell} python3292+```{code-cell} ipython3299293μ300294```301295302-```{code-cell} python3296+```{code-cell} ipython3303297Σ304298```305299306-### Back to The Administrator's Problem300+### Back to the administrator's problem307301308302Now let's turn to the grant administrator's problem.309303310304Here the array of311305numbers of $i$ objects in the urn is312306$\left(157, 11, 46, 24\right)$.313307314-```{code-cell} python3308+```{code-cell} ipython3315309K_arr = [157, 11, 46, 24]316310urn = Urn(K_arr)317311```318312319313Let's compute the probability of the outcome $\left(10, 1, 4, 0 \right)$.320314321-```{code-cell} python3315+```{code-cell} ipython3322316k_arr = [10, 1, 4, 0]323317urn.pmf(k_arr)324318```325319326320We can compute probabilities of three possible outcomes by constructing a 3-dimensional327321arrays `k_arr` and utilizing the method `pmf` of the `Urn` class.328322329-```{code-cell} python3323+```{code-cell} ipython3330324k_arr = [[5, 5, 4 ,1], [10, 1, 2, 2], [13, 0, 2, 0]]331325urn.pmf(k_arr)332326```333327334328Now let's compute the mean and variance-covariance matrix of $X$ when $n=6$.335329336-```{code-cell} python3330+```{code-cell} ipython3337331n = 6 # number of draws338332μ, Σ = urn.moments(n)339333```340334341-```{code-cell} python3335+```{code-cell} ipython3342336# mean343337μ344338```345339346-```{code-cell} python3340+```{code-cell} ipython3347341# variance-covariance matrix348342Σ349343```350344351345

We can simulate a large sample and verify that sample means and covariances closely approximate the population means and covariances.

352346353-```{code-cell} python3347+```{code-cell} ipython3354348size = 10_000_000355349sample = urn.simulate(n, size=size)356350```357351358-```{code-cell} python3352+```{code-cell} ipython3359353# mean360354np.mean(sample, 0)361355```362356363-```{code-cell} python3357+```{code-cell} ipython3364358# variance covariance matrix365359np.cov(sample.T)366360```367361368362Evidently, the sample means and covariances approximate their population counterparts well.369363370-### Quality of Normal Approximation364+### Quality of normal approximation371365372366

To judge the quality of a multivariate normal approximation to the multivariate hypergeometric distribution, we draw a large sample from a multivariate normal distribution with the mean vector  and covariance matrix for the corresponding multivariate hypergeometric distribution and compare the simulated distribution with the population multivariate hypergeometric distribution.

373367374-```{code-cell} python3368+```{code-cell} ipython3375369sample_normal = np.random.multivariate_normal(μ, Σ, size=size)376370```377371378-```{code-cell} python3372+```{code-cell} ipython3379373def bivariate_normal(x, y, μ, Σ, i, j):380374381375 μ_x, μ_y = μ[i], μ[j]

@@ -392,7 +386,7 @@ def bivariate_normal(x, y, μ, Σ, i, j):

392386 return np.exp(-z / (2 * (1 - ρ**2))) / denom393387```394388395-```{code-cell} python3389+```{code-cell} ipython3396390@jit397391def count(vec1, vec2, n):398392 size = sample.shape[0]

@@ -404,7 +398,7 @@ def count(vec1, vec2, n):

404398 return count_mat405399```406400407-```{code-cell} python3401+```{code-cell} ipython3408402c = urn.c409403fig, axs = plt.subplots(c, c, figsize=(14, 14))410404

@@ -454,7 +448,7 @@ The null hypothesis is that the sample follows normal distribution.

454448455449> `normaltest` returns an array of p-values associated with tests for each $k_i$ sample.456450457-```{code-cell} python3451+```{code-cell} ipython3458452test_multihyper = normaltest(sample)459453test_multihyper.pvalue460454```

@@ -463,7 +457,7 @@ As we can see, all the p-values are almost $0$ and the null hypothesis is soundl

463457464458By contrast, the sample from normal distribution does not reject the null hypothesis.465459466-```{code-cell} python3460+```{code-cell} ipython3467461test_normal = normaltest(sample_normal)468462test_normal.pvalue469463```