Method of moments

Peter Ralph

https://uodsci.github.io/dsci345

import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['figure.figsize'] = (15, 8)
import numpy as np
import pandas as pd
import dsci345 as dsci

rng = np.random.default_rng(seed=123)

$$\renewcommand{\P}{\mathbb{P}} \newcommand{\E}{\mathbb{E}} \newcommand{\var}{\text{var}} \newcommand{\sd}{\text{sd}}$$ This is here so we can use \P and \E and \var and \sd in LaTeX below.

Fitting distributions to data

lam = 1
npanels = 1000
nd = rng.poisson(lam * rng.exponential(size=npanels), size=npanels)
defects = pd.DataFrame({ "value" : np.arange(11, dtype='int') })
defects['observed'] = [np.sum(nd == k) for k in defects['value']]

Back to the factory

In our Poisson example, we had data like this for the distribution of numbers of defects per solar panel. Our goal was to get a good model for this number. How did we do that?

defects

	value	observed
0	0	514
1	1	252
2	2	127
3	3	53
4	4	26
5	5	17
6	6	4
7	7	2
8	8	2
9	9	3
10	10	0

Well, here's what we did:

The mean number of defects per panel is 1, and they're "counts", so... maybe Poisson? With mean 1?

In other words, we used the Method of Moments:

1. Pick a particular form of the distribution.
2. Choose parameter values for the distribution so that the "moments" (here, the mean) match.

Our next step was to say

Gee, it's more spread out ("overdispersed") than that... maybe an Exponential mixture of Poissons? (i.e., Poisson but where the mean is random, from the Exponential distribution)

Let's have a look at applying the method of moments here, too.

Goodness-of-fit

Do these look similar? How or how not?

obs_mean = np.sum(defects['value'] * defects['observed']) / np.sum(defects['observed'])
obs_sd = np.sqrt(
    np.sum((defects['value'] ** 2) * defects['observed']) / np.sum(defects['observed'])
    - obs_mean ** 2
)
x = rng.poisson(obs_mean, size=np.sum(defects['observed']))
defects['expected'] = [np.sum(x==u) for u in defects['value']]
defects

	value	observed	expected
0	0	514	382
1	1	252	368
2	2	127	185
3	3	53	55
4	4	26	9
5	5	17	1
6	6	4	0
7	7	2	0
8	8	2	0
9	9	3	0
10	10	0	0

A plot of the numbers in the previous slide: observed counts and the expected counts under the Poisson model, as lines. The "observed" line starts higher than the "expected, but then drops below it.

plt.plot(defects['value'], defects['observed'], label=f'observed, mean={obs_mean}, sd={obs_sd:.2f}')
plt.xlabel("count"); plt.ylabel("frequency")
plt.plot(defects['value'], defects['expected'], label=f'poisson, mean={obs_mean}, sd={np.sqrt(obs_mean):.2f}')
plt.legend();

Variance of a mixture

We'd like to fit two moments (mean and SD), so let's use a model with two parameters: $$\begin{aligned} \text{error rate: } R &\sim \text{Gamma}(\text{scale}=\theta, \text{shape}=k) \\ \text{number of defects: } X &\sim \text{Poisson}(R) . \end{aligned}$$ We'd like to see if we can choose $\theta$ and $k$ so the mean and SD match the observed values.

Law of Total Variance: $$ \var[X] = \E[\var[X|R]] + \var[\E[X|R]] $$

Here we have $$\begin{aligned} \var[X|R] &= R \\ \E[X|R] &= R \\ \E[R] &= \theta k \\ \var[R] &= \theta^2 k \end{aligned}$$ and so $$ \var[X] = \theta k (1 + \theta ) , $$ i.e., the variance is inflated by $(1 + \theta)$.

Okay, so we want to find $\theta$ and $k$ so that $$\begin{aligned} \E[X] = \E[R] = \theta k = \text{observed mean} \\ \sd[X] = \sqrt{\theta k (1 + \theta)} = \text{observed SD} . \end{aligned}$$

obs_mean, obs_sd

(np.float64(0.935), np.float64(1.3478779618348242))

theta = obs_sd**2  / obs_mean - 1
k = obs_mean / theta
print(f"Estimates: theta = {theta:.2f}, k = {k:.2f}.")

Estimates: theta = 0.94, k = 0.99.

Goodness-of-fit, take 2

Here's our fit, having used the method of moments to fit a Gamma-Poisson mixture to the data:

R = rng.gamma(shape=k, scale=theta, size=np.sum(defects['observed']))
X = rng.poisson(R, size=np.sum(defects['observed']))
defects['expected'] = [np.sum(X==u) for u in defects['value']]
defects

	value	observed	expected
0	0	514	539
1	1	252	234
2	2	127	112
3	3	53	56
4	4	26	30
5	5	17	11
6	6	4	11
7	7	2	3
8	8	2	1
9	9	3	0
10	10	0	1

A plot of the numbers in the previous slide, as before, except this time, the 'observed' and 'expected' lines match closely.

plt.plot(defects['value'], defects['observed'], label=f'observed, mean={obs_mean:.2f}, sd={obs_sd:.2f}')
plt.xlabel("count"); plt.ylabel("frequency")
plt.plot(defects['value'], defects['expected'], label=f'poisson, mean={np.mean(X):.2f}, sd={np.std(X):.2f}')
plt.legend();