## 1.

By linearity of expectation,

$\mathbb{E}\left[a^{\intercal}X\right]=\sum_{i}a_{i}\mathbb{E}X_{i}=\sum_{i}a_{i}\mu_{i}=a^{\intercal}\mu.$

Since

$AX=\begin{pmatrix}A_{1*}\\ \vdots\\ A_{n*} \end{pmatrix}X=\begin{pmatrix}A_{1*}X\\ \vdots\\ A_{n*}X \end{pmatrix},$

it follows that

$\mathbb{E}\left[AX\right]=\begin{pmatrix}A_{1*}\mu\\ \vdots\\ A_{n*}\mu \end{pmatrix}=A\mu.$

Similarly,

$\mathbb{V}(a^{\intercal}X)=\sum_{ij}a_{i}a_{j}\operatorname{Cov}(X_{i},X_{j})=a^{\intercal}\Sigma a.$

It follows that

$[\mathbb{V}(AX)]_{ij}=(A_{i*}X,A_{j*}X)=\sum_{uv}A_{iu}A_{jv}\operatorname{Cov}(X_{u},X_{v})=[A\Sigma A^{\intercal}]_{ij}.$

## 2.

The log-probability of a multinomial is

$\log f(X;p)=\log(n!)+\sum_{j}-\log(X_{j}!)+X_{j}\log(p_{j}).$

Therefore,

$\mathbb{E}\left[\frac{\partial^{2}\log f(X;p)}{\partial p_{j}^{2}}\right]=\mathbb{E}\left[-\frac{X_{j}}{p_{j}^{2}}\right]=-\frac{n}{p_{j}}$

so that the Fisher information matrix is

$I(p)=n\operatorname{diag}(1/p_{1},\ldots,1/p_{k}).$

## 3.

import numpy as np

def sample_multinomial(n, p):
"""Samples a multinomial distribution.

Parameters
----------
n : int
Number of bins.
p : array_like, shape=[k]
Bin probabilities (summing to one).

Returns
-------
counts : array_like, shape=[k]
Realized bin counts.
"""
k, = p.shape
cdf = np.cumsum(p)
samples = np.random.uniform(low=0.0, high=1.0, size=n)
indices = np.searchsorted(cdf, samples)
indices = np.concatenate([np.arange(k), indices])
_, counts = np.unique(indices, return_counts=True)
counts -= 1
return counts


## 4.

def sample_multivariate_normal(mean, cov):
"""Samples a multivariate normal distribution.

Parameters
----------
mean : array_like, shape=[k]
Mean.
cov : array_like, shape=[k, k]
Covariance matrix.

Returns
-------
z : array_like, shape=[k]
Realized samples.
"""
k, = mean.shape
sqrt_cov = sqrtm(cov)
samples = np.random.randn(k)
return mean + sqrt_cov @ samples


## 5.

The true correlation between $$X1$$ and $$X2$$ is $$1 / \sqrt{2}$$.

100 points were sampled:

Various quantities were computed from these samples: (see Appendix for code)

Sample mean:
[3.19678877 8.08354934]

Sample covariance matrix:
[[0.68362962 0.67545853]
[0.67545853 1.62768145]]

Sample correlation:
0.6403295531123555

Fisher's 95% confidence interval:
(0.5047523958530407, 0.7450792923252504)

Bootstrap 95% confidence interval:
(0.5124584072915059, 0.7682006989332051)


Fisherâ€™s method and the (nonparametric) bootstrap yield comparable results.

## 6.

1000 points were sampled:

Various quantities were computed from these samples: (see Appendix for code)

Sample mean:
[3.04724353 8.05340182]

Sample covariance matrix:
[[0.98527055 0.95830459]
[0.95830459 1.97142287]]

Sample correlation:
0.6876000872761875

Fisher's 95% confidence interval:
(0.6527346353163089, 0.7195590209752238)

Bootstrap 95% confidence interval:
(0.6532707565958478, 0.7219294179565272)


Fisherâ€™s method and the (nonparametric) bootstrap yield comparable results.

## Appendix: Code for Exercises 5 and 6

n_points = 100
mean = np.array([3., 8.])
cov = np.array([[1., 1.], [1., 2.]])

np.random.seed(1)

def sample_cov_and_corr(points):
sample_cov = np.cov(points, ddof=1, rowvar=False)
sample_std = np.sqrt(np.diag(sample_cov))
sample_corr = ((sample_cov / sample_std).T / sample_std).T
return sample_cov, sample_corr

points = np.array(
[sample_multivariate_normal(mean, cov) for _ in range(n_points)])
sample_mean = np.mean(points, axis=0)
sample_cov, sample_corr = sample_cov_and_corr(points)
sample_corr = sample_corr[0, 1]

def fisher_func(r):
return 0.5 * (np.log(1 + r) - np.log(1 - r))

def fisher_func_inv(z):
return (np.exp(2. * z) - 1) / (np.exp(2. * z) + 1)

theta = fisher_func(sample_corr)
a = theta - 2. / np.sqrt(n_points - 3)
b = theta + 2. / np.sqrt(n_points - 3)
fisher_lb = fisher_func_inv(a)
fisher_ub = fisher_func_inv(b)

bstrap_sample_corrs = []
for _ in range(1000):
indices = np.random.randint(0, high=points.shape[0], size=points.shape[0])
new_points = points[indices]
_, bstrap_sample_corr = sample_cov_and_corr(new_points)
bstrap_sample_corr = bstrap_sample_corr[0, 1]
bstrap_sample_corrs.append(bstrap_sample_corr)
se = np.std(bstrap_sample_corrs)
bstrap_lb = sample_corr - 2. * se
bstrap_ub = sample_corr + 2. * se

print("""
Sample mean:
{}

Sample covariance matrix:
{}

Sample correlation:
{}

Fisher's 95% confidence interval:
({}, {})

Bootstrap 95% confidence interval:
({}, {})
""".format(sample_mean, sample_cov, sample_corr, fisher_lb, fisher_ub,
bstrap_lb, bstrap_ub))