---
title: Modeling Income and Wealth Distributions II --- Fitting and Checking
author: 36-313, Fall 2021
date: 7 September 2021 (Lecture 4)
output:
  pdf_document:
    toc: true
    number_sections: true	
    toc_depth: 4
bibliography: locusts.bib
header-includes:
  - \DeclareMathOperator*{\argmin}{argmin}
  - \DeclareMathOperator*{\argmax}{argmax}
---

\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right) }
\newcommand{\Expect}[1]{\mathbb{E}\left( #1 \right) }
\newcommand{\Var}[1]{\mathrm{Var}\left( #1 \right) }
\newcommand{\Cov}[1]{\mathrm{Cov}\left( #1 \right) }
\newcommand{\myexp}[1]{\exp{\left\{ #1 \right\}}}
\newcommand{\xmin}{x_{\min}}
\newcommand{\TrueParam}{\theta^*}
\newcommand{\Jacobian}{\mathbf{j}}

```{r, include=FALSE}
# Set up global options for displaying code and so forth
library(knitr)
opts_chunk$set(size="small", 
               tidy=TRUE, tidy.opts=list(comment=FALSE),
               cache=TRUE, cache.comments=FALSE, autodep=TRUE,
               warning=FALSE, message=FALSE, error=FALSE, echo=FALSE)
```


```{r}
# Data from HW on the US income distribution in 2019, from the Current
# Popilation Survey (in dollars per year)
cps.median <- 68703
cps.median.se <- 550
cps.mean <- 98088
cps.mean.se <- 633
```

In this lecture, we'll look at how we actually fit models, like the log-normal
or Pareto, to data, at some ways of checking the fit.  We'll hold off on
talking (much) about processes which produce distributions like
these until a later date.

\clearpage

# Fitting

There are two big ways of fitting probability distributions to data.  One is to
try to match some features or aspects of the data: pick some summary
statistics, calculate them on the data, and then tune the parameters until they
match what we see.  The other leading type of estimator is the method of
maximum likelihood, where we use the model to calculate the probability of the
data set we got at various possible parameters, and adjust the parameters until
the data is as "likely" as possible.  Matching summary features is older,
easier to understand, and widely applicable, but maximum likelihood (usually)
gives better estimates, when we can use it.   We'll start with
feature matching.

## Common Notation

We have a probability model which specifies the probability density
$f(x;\theta)$ as a function of the parameter $\theta$.  $\theta$ will usually
be a vector, e.g., for the log-normal distribution $\theta = (\mu, \sigma^2)$
with $\mu = \Expect{log{X}}$ and $\sigma^2 = \Var{\log{X}}$.  We'll write
$d$ for the number of dimensions or coordinates in $\theta$.  ($d=2$ for the
log-normal.)  We'll need a symbol to distinguish the correct value
of $\theta$ from the others; I'll write $\TrueParam$.

Our data set consists of measurements of $X$ for $n$ members of the population,
which we'll imaginative write as $X_1, X_2, \ldots X_n$.  (The capital letters
here are reminder that the data are random.)  If we need to abbreviate we'll
write $X_{1:n}$ for the whole data, $X_{1:3}$ for the first 3 data points, etc.

We'll assume
that every individual in the data is statistically independent of every
other, so
\begin{equation}
f(x_{1:n}; \theta) = \prod_{i=1}^{n}{f(xi;\theta)}
\end{equation}
(Independence isn't strictly necessary but it simplifies our math in some
places.)

Our _main_ goal is to go from $X_1, \ldots X_n$ to an estimate of the true
$\theta$, say $\hat{\theta}$.  We should really write this as
$\hat{\theta}(X_{1:n}$ because it's a function of the data, but we'll suppress
that to keep the notation under control.  Because the data are random,
the estimate $\hat{\theta}$ is also random.  But we'd like it to be _close_
to $\TrueParam$, and we'd like it to get closer as we get more
data --- we want $\hat{\theta}$ to **converge on** $\TrueParam$, $\hat{\theta}(X_{1:n}) \rightarrow \TrueParam$ as $n\rightarrow\infty$.
Estimators which converge on the truth are called
**consistent**[^consistency].

[^consistency]: The name makes more sense historically; see [http://bactra.org/notebooks/consistency-pac.html].




\clearpage

# Fitting by Matching Features or Summary Statistics

In this approach, we don't work directly with the individual observations, but
only with selected **features** or **statistics** which we have calculated from
the data (or have had handed to us).  These features or statistics are
_functions_ of the data, and they're supposed to summarize important and
informative aspects of it, while being smaller and easier to understand than
the full data.

In algebra, we have $m$ statistics $G_1, G_2, \ldots G_m$, each defined
as a function of the data:
\begin{equation}
G_k = g_k(X_{1:n}), ~ k \in 1:m
\end{equation}
(Notice that because the data are random, the values of the summary statistics
are also random variables, hence the capital letters $G_k$, while the
_functions_ used to calculate them are fixed, hence the lower-case $g_k$.)
We can also think of a single $m$-dimensional vector-valued statistic, $G = g(X_{1:n})$.

Examples of summary statistics include, but are not limited to:

- **Moments**, such as the sample mean $\overline{x}$, the sample mean square $\overline{x^2}$, up to the kth moment $\overline{x^k}$ and beyond
- **Central moments** such as the sample variance $\overline{(x - \overline{x})^2}$ (or the square root of this, the sample standard deviation) or the kth central moment $\overline{(x-\overline{x})^k}$
- **Percentiles** or other quantiles, including the median
- Ratios of or differences of any of these
- etc., etc., etc.: anything we can calculate from the data can, in principle, be used as a summary statistic, e.g., the Gini index.

We have seen all of these as ways of just describing the data.

So far, these summary statistics are _just_ functions of the data, and don't
involve any probability model at all.  But, once we have a probability model,
we can work out what these quantities _ought_ to be, for the whole
population or the distribution.  At its most basic, if we consider
the sample mean
\begin{equation}
\overline{x} \equiv \frac{1}{n}\sum_{i=1}^{n}{x_i}
\end{equation}
the corresponding theoretical quantity is the expected value,
\begin{equation}
\Expect{X} = \int_{-\infty}^{\infty}{x f(x;\theta) dx}
\end{equation}
Notice that the right-hand side of that last equation, once we've gone and done
the integral, is going to be a function of the parameter $\theta$.  For
instance, for the log-normal distribution, $\Expect{X} = e^{\mu + \sigma^2/2}$.

_All_ of the kinds of summary statistics I mentioned have theoretical,
model-based counter-parts for the whole distribution.  Once we have (say) a
model for income distribution, we can use it to calculate what (say) the ratio
between the median income and the 99th percentile of income _should_ be, in the
whole population, as a function of $\theta$. 

> _For every choice of summary statistics $g(X_1, \ldots X_n)$ and probability model $f(x;\theta)$, there is a corresponding theoretical value $\gamma(\theta)$_

(Notice I am using the Greek letter $\gamma$, which is the counter-part to the
Roman letter $g$, to stand for the theoretical quantity.)

Here is the key idea of the feature-matching approach to estimating parameters:

> Adjust the parameter $\theta$ so that the theoretical values $\gamma_1(\theta), \ldots \gamma_m(\theta)$ all match the observed summary statistics $g_1(X_{1:n}), \ldots g_m(X_{1:n})$

One way to think of this is to set up a system of equations:
\begin{eqnarray}
\gamma_1(\theta) & = & G_1 = g_1(X_1, \ldots X_n)\\
\gamma_2(\theta) & = & G_2 = g_2(X_1, \ldots X_n)\\
& \vdots & \\
\gamma_m(\theta) & = & G_m = g_m(X_1, \ldots X_n)
\end{eqnarray}
All the right-hand sides are known (or at least stuff we can calculate
from the data), and all the left-hand sides are functions of the unknown
parameter $\theta$.  So we try to solve for $\theta$, and call the result
our estimate $\hat{\theta}$.

Remember that $\theta$ is a vector, with dimension $d$, so we've got $m$
equations and $d$ unknowns.  And we remember some rules about counting equations and unknowns from algebra:

- If $m < d$, there are fewer equations than unknowns, and there won't be a
  unique solution.  (The system of equations is "under-determined".)  That is, we can't possibly hope to get _the_ unique
  $\theta$ this way (though we might put some useful constraints on it).  We
  need at least one feature/statistic per parameter.
- If $m=d$,  there will (generally[^weasle]) be a unique solution.  (I will explain that
  weasily "hope" in a moment.)  --- Again, we want at least one feature per
  parameter.
- If $m > d$, there will generally not be _one_ value of $\theta$ which can
  solve _all_ the equations.  (The system of equations is "over-determined".)  We therefore seem to be in trouble if we use _more_ information to estimate
  the parameters rather than less, but that sounds stupid.
  
[^weasle]: I put in the weasle-word "generally" because there can be some  exceptions, say where random sampling noise gives a value for one of the $G_k$ which just can't _exactly_ match $\gamma_k(\theta)$ for any $\theta$.  (This can happen when we have continuous parameters for discrete observations and a small sample size, for example.)  This sort of annoying edge case will be handled by the same least-squares trick I'm about to explain for over-determined systems.

Fortunately, about two hundred years ago the Ancestors realized that there was
a sensible way to get an estimate of $d$ parameters from $m>d$ observations,
even when the equations can't be solved uniquely[^premoderndata]: this is the method of
least squares [@Farebrother-on-linear-relationships].  For each quantity, we calculate the squared error
$(G_k - \gamma_k(\theta))^2$, we add up, and we adjust $\theta$ to minimize
the sum of squared errors.  In modern symbols:
\begin{equation}
\hat{\theta} = \argmin_{\theta}{\sum_{k=1}^{m}{(G_k - \gamma_k(\theta))^2}}
\end{equation}
Since the sum of squares is necessarily $\geq 0$, if there's a $\theta$ which
_does_ exactly match the observed values of the statistics, that's automatically
what this method will pick out, but this _also_ works even when $m > d$.

[^premoderndata]: Before that, for thousands of years, the best scientific practice was to try to determine which $d$ observations were the most accurate, use them to solve for the parameters, and discard the rest of the data, or at best use it to check later [@Farebrother-on-linear-relationships].  This was wasteful of data!

(Of course, there is nothing magical about treating all the summary statistics
equally.  We could also do weighted least squared
\begin{equation}
\hat{\theta} = \argmin_{\theta}{\sum_{k=1}^{m}{(G_k - \gamma_k(\theta))^2 w_k}}
\end{equation}
with weights $w_k > 0$, and this is often useful in practice, say if one
statistic is more important for applications, or better-measured.  See Complementary Problem 2.)

## A little worked example with the log-normal distribution

This has all been very abstract so it's good to give a concrete worked example
with a distribution we're using a lot, like the log-normal.  In Lecture 3,
we saw that the expected value of the log-normal is $e^{\mu + \sigma^2/2}$.
We also saw there that the theoretical median is just $e^{\mu}$.  The sample  mean and median are of course among our very first and simplest summary statistics, let's say $G_1$ and $G_2$.
So _one_ way of estimating the parameters of the log-normal distribution would
be to adjust $\mu$ and $\sigma^2$ until it matched the observed values
of the mean and the median:
\begin{eqnarray}
\gamma_1(\mu, \sigma^2) = e^{\mu+\sigma^2/2} & = & G_1 = \text{sample mean}\\
\gamma_2(\mu, \sigma^2 = e^{\mu} & = & G_2 = \text{sample median}
\end{eqnarray}
In this case, $d=m$ so there'll be a unique solution, and you can, in fact,
work it out explicitly (because the functions involved are very well-behaved):
\begin{eqnarray}
\hat{\mu} & = & \log{G_2}\\
\widehat{\sigma^2} & = & 2\log{(G_1/G_2)}
\end{eqnarray}

It is _also_ true that, for a log-normal distribution, the Gini index is
given by $2\Phi(\sigma/\sqrt{2})-1$, where $\Phi$ is the CDF of the standard
Gaussian distribution.  If we had $m=3$ summary statistics available
to us, which were the mean, the median, and the Gini index, we'd need to
set up the least-squares problem I mentioned above[^discardinginfo].

[^discardinginfo]: Alternatively, we could pick any two of the three summary statistics, use them to estimate $\mu$ and $\sigma^2$, and then check that against the predicted value for the third statistic...

## When will feature matching work, and how well?

We pick our features or summary statistics $g_1, \ldots g_m$; we calculate them
on data $X_1, \ldots X_n$ to get numbers $G_1, \ldots G_m$; we try to match
them to the corresponding theoretical values $\gamma_1(\theta), \ldots
\gamma_m(\theta)$.  We get a value $\hat{\theta}$ for the parameter, either by
exactly solving the system of equations $\gamma_k(\theta) = G_k$, or by least
squares.  How well does this $\hat{\theta}$ approximate the true $\theta$?
Better said, what makes for better or worse estimates from this approach?

There are two _very_ important requirements which have to hold if feature-matching is to have any chance of working:

1. **Convergence**: As we get more and more data, our summary
statistics or features have to approach their theoretical limits:
\begin{equation}
\lim_{n\rightarrow\infty}{g_k(X_{1:n})} \rightarrow \gamma_k(\TrueParam)
\end{equation}
This is where probability results like the law of large numbers come in: it
needs to be the case that increasingly big samples become closer and closer
approximations to the whole data-generating distribution.  Most of the common summary statistics have this property!
2. **Identifiability**: Differences in the parameters _must_ make a difference to the theoretical values of the features.  That is, if we have two different parameter vectors $\theta \neq \theta^{\prime}$, there has to be _some_ $\gamma_k$ where $\gamma_k(\theta) \neq \gamma_k(\theta^{\prime})$.  Said differently the vector-valued function $\gamma(\theta) = (\gamma_1(\theta), \ldots \gamma_m(\theta))$ must be invertible, so $\gamma^{-1}$ must be well-defined[^continuity].  Checking this is, in
principle, just a matter of doing math with our favorite probability model.

[^continuity]: We sometimes demand a stronger condition, that $\gamma^{-1}$ isn't just well-defined, but also "smooth", say continuous, or even differentiable.  This is so that small errors in the summaries will "track back" through $\gamma^{-1}$ to small changes in the parameters.  It is genuinely hard to come up with examples of real models and summaries where $\gamma^{-1}$ exists but _isn't_ smooth, however.

If convergence fails, matching summary statistics to theoretical values is
senseless.  If identifiability or invertibility fails, even perfect matching
can't pick out ("identify") the correct parameter value.

So, let's assume that we're using summary statistics which converge, and that
the model parameters are identifiable from them.  There usually going to be
_lots_ of different combinations of features with these properties.  Some of
them will work better than others, i.e., give more accurate and precise
estimates.  Two aspects of the features will improve the quality of the
estimates:

3. **Sensitivity** to the parameters: small changes to the parameter vector should lead to big changes to the theoretical values of the features.  More precisely, if we think of $\theta$ as a vector $\theta=(\theta_1, \theta_2, \ldots \theta_d)$, we'd like all the partial derivatives $\frac{\partial \gamma_k}{\partial \theta_i}$ to be large.  In fact, all else being equal, $\Var{\hat{\theta}}$ will be proportional to $\left(\frac{\partial \gamma_k}{\partial \partial \theta_i}\right)^{-2}$.  In words: we'll estimate $\theta$ better if even small changes in $\theta$ are easily to detect in the features.
4. **Precision** of the statistics: the smaller the variance of each summary statistic, $\Var{G_k}$, the more precise the estimates will be.  In fact,
all else being equal, $\Var{\hat{\theta}}$ will be proportional to $\Var{G_k}$.
In words: we'll estimate $\theta$ better if we're matching features which aren't themselves very noisy.

These last two points are, I hope, ones which make sense intuitively, but the
next sub-section will explain them in some mathematical detail --- and cash out
my vague "all else being equal" in actual formulas.

```{r, fig.cap="Dotted curve: theoretical median income as a function of the $\\mu$ parameter in a log-normal model (namely $e^{\\mu}$).  Black horizontal line: measurement of median income by the Current Population Survey in 2019.  Black vertical line: implied value of $\\mu$ from matching.  Grey lines horizontal lines: a random scatter around the point estimate of the median, with the same uncertainty as the reported standard error for the median. Grey vertical lines: scatter of implied $\\mu$ values, showing the induced uncertainty in this estimate of the parameter.  Notice the slope ${\\frac{\\partial \\gamma}{\\partial \\mu}}$ (dashed blue line) helps propagate uncertainty in the observed statistic into uncertainty in the estimated parameter."}
curve(exp(x), from=10, to=12, xlab=expression(mu), ylab="Median income ($/yr)",
      lty="dotted")
segments(x0=10, x1=log(cps.median), y0=cps.median, y1=cps.median)
segments(x0=log(cps.median), x1=log(cps.median), y0=0, y1=cps.median)
# Variance of a uniform distribution over [a,b] is (1/12)*(b-a)^2
# so if we know the desired variance v then the range is (b-a)^2=12*v
# or (b-a)=sqrt(12*v)
cps.median.noise.range <- cps.median.se*sqrt(12)
alternative.medians <- runif(n=5, min=cps.median-cps.median.noise.range/2,
                                  max=cps.median+cps.median.noise.range/2)
alternative.muhats <- log(alternative.medians)
for (i in seq_along(alternative.medians)) {
  segments(x0=10, x1=alternative.muhats[i], y0=alternative.medians[i],
           y1=alternative.medians[i], col="grey")
  segments(x0=alternative.muhats[i], x1=alternative.muhats[i],
          y0=0, y1=alternative.medians[i], col="grey")
}
# Tangent line
  # The curve is e^{\mu} whose derivative is e^{\mu} so slope = b = cps.median
  # Intercept a needs to give 
```



## Asymptotics for feature matching

> _This sub-section is a bit more mathematically involved than the rest of this section; it's useful stuff but not quite so crucial to what we'll do in the rest of the course, so you might want to skim it first, then come back to it later._

I said above that we can think of fitting a model by matching features as
solving a least squares problem:
\begin{equation}
\hat{\theta} = \argmin_{\theta}{\sum_{k=1}^{m}{(G_k - \gamma_k(\theta))^2}}
\end{equation}
It'll be convenient to have a name for the function we're minimizing, so
let's call it $M$:
\begin{equation}
M(G, \theta) \equiv \sum_{k=1}^{m}{(G_k - \gamma_k(\theta))^2}
\end{equation}
$M$ is a sum of squares, so it must be $\geq 0$.  The only way we can have
$M=0$ is if $G_k = \gamma_k(\theta)$ for all $k \in 1:m$.  That is, the
smallest conceivable value of $M$ is attained when we match all the observed
features exactly.

### Consistency

Let's use the assumption I called "convergence" and numbered (1) above: as $n
\rightarrow\infty$, $G_k \rightarrow \gamma_k(\TrueParam)$.  (Remember
$\TrueParam$ is the true value of $\theta$, which we're trying to learn.)
So
\begin{equation}
M(G, \TrueParam) \rightarrow 0
\end{equation}
For any other $\theta \neq \TrueParam$,
\begin{equation}
M(G, \theta) \rightarrow \sum_{k=1}^{m}{(\gamma_k(\theta) - \gamma_k(\TrueParam))^2}
\end{equation}
Now use the assumption I called "identifiability" or "invertibility"  and numbered (2) above: if $\theta \neq \TrueParam$, then $\gamma_k(\theta) \neq \gamma_k(\TrueParam)$ for (at least) some $i$.  Thus, with this extra assumption, if $\theta \neq \TrueParam$,
\begin{equation}
M(G, \theta) \rightarrow \sum_{k=1}^{m}{(\gamma_k(\theta) - \gamma_k(\TrueParam))^2} > 0
\end{equation}

Hence, under these assumptions, as $n\rightarrow\infty$, $\TrueParam$ will come
closer and closer to being the minimizer of $M(G,\theta)$, and any _other_
parameter value will look worse and worse by comparison.  Conclusion[^sloppy]:

> _Assuming convergence and identifiability, then $\hat{\theta} \rightarrow \TrueParam$ as $n\rightarrow \infty$.

[^sloppy]: If you worry that this isn't a totally air-tight proof, you're not wrong.  The missing extra assumptions ("regularity conditions") to make it rigorous aren't too onerous, but they do take us too far afield.  If you really want to know more, see @van-der-Vaart-asymptotic-stats, chapter 5.


### Uncertainties (variances) and the "sandwich variance" formula

We can say more about _how_ $\hat{\theta} \rightarrow \TrueParam$ if we use a
little calculus.  The first point is that $\hat{\theta}$ minimizes a function,
and we know from calculus that the derivative is zero at a
minimum[^smoothinteriorminimum].  So
\begin{equation}
\left. \frac{\partial M}{\partial \theta_i}\right|_{\theta=\hat{\theta}} = 0
\end{equation}
Remember that $M(G,\theta) = \sum_{k=1}^{m}{(G_k - \gamma_k(\theta))^2}$, so
this becomes
\begin{equation}
2\sum_{k=1}^{m}{(G_k-\gamma_k(\hat{\theta}))\left. \frac{\partial \gamma_k}{\partial \theta_i}\right|_{\theta=\hat{\theta}}} = 0
\end{equation}
(I have my reasons for not dropping the factor of 2 here.)

[^smoothinteriorminimum]: Assuming the function is smooth, the minimum is located in the interior of the domain, etc., etc.


It's annoying to keep track of all the derivatives separately, so we
bundle them into a vector, the **gradient** of the function (with respect to $M$):
\begin{equation}
\nabla M(G, \hat{\theta}) = 0
\end{equation}
I will leaving writing this out in terms of the vector $G-\gamma(\hat{\theta})$
and a matrix of partial derivatives as an exercise.

Now we use another idea from basic calculus, that of a Taylor expansion or
Taylor approximation: $h(u) \approx h(u_0) + (u-u_0)\left.\frac{dh}{du}\right|_{u=u_0}$.  Applied here, and expanding around the true parameter value $\TrueParam$,
\begin{equation}
\nabla M(G, \hat{\theta}) \approx \nabla M(G, \TrueParam) + (\nabla \nabla M(G, \TrueParam)) (\hat{\theta} - \TrueParam)
\end{equation}
Here, $\nabla\nabla M(G, \TrueParam)$ is the matrix of second partial
derivatives (with respect to the coordinates of $\theta$), also called the
**Hessian** of $M$.

Putting these last two equations together,
\begin{eqnarray}
0 & \approx & \nabla M(G, \TrueParam) + (\nabla \nabla M(G, \TrueParam)) (\hat{\theta} - \TrueParam)\\
0 & \approx & (\nabla \nabla M(G, \TrueParam))^{-1}\nabla M(G, \TrueParam) + (\hat{\theta} - \TrueParam)\\
\hat{\theta} & \approx & \TrueParam - (\nabla \nabla M(G, \TrueParam))^{-1}\nabla M(G, \TrueParam)
\end{eqnarray}

As $n\rightarrow\infty$, by convergence, $\nabla M(G, \TrueParam) \rightarrow
0$ (why?), so this again tells us that $\hat{\theta} \rightarrow \TrueParam$.
It tells us a little more, too.  As $n$ grows, $\nabla \nabla M(G, \TrueParam)$ approaches a limiting Hessian matrix, say $\mathbf{h}$, and then
\begin{equation}
\hat{\theta} \approx \TrueParam - \mathbf{h}^{-1} \nabla M(G, \TrueParam)
\end{equation}
so all the randomness in $\hat{\theta}$ comes from the randomness in $\nabla M(G, \TrueParam)$.  In fact, using the rules for algebra with variances,
\begin{equation}
\Var{\hat{\theta}} \approx \mathbf{h}^{-1}\Var{\nabla M(G,\TrueParam)} \mathbf{h}^{-1}
\end{equation}
This is called the **sandwich variance formula**.

Let's explore the two parts of this, the (limiting) Hessian $\mathbf{h}$ and
the variance of the gradient $\Var{\nabla M(G,\TrueParam)}$.

As $n\rightarrow\infty$, $G(X_{1:n}) \rightarrow \gamma(\TrueParam)$, so
\begin{equation}
M(G, \theta) \approx \sum_{k=1}^{m}{(\gamma_k(\theta) - \gamma_k(\TrueParam))^2}
\end{equation}
We can just take the second partial derivatives:
\begin{eqnarray}
\frac{\partial^2 M}{\partial \theta_i \partial \theta_k} & \approx &\frac{\partial}{\partial \theta_k} 2\sum_{k=1}^{m}{(\gamma_k(\theta) - \gamma_k(\TrueParam))\frac{\partial \gamma_k}{\partial \theta_i}}\\
& \approx & 2\sum_{k=1}^{m}{(\gamma_k(\theta) - \gamma_k(\TrueParam))\frac{\partial^2 \gamma_k}{\partial \theta_i \partial \theta_k} + \frac{\partial\gamma_k}{\partial \theta_i}\frac{\partial\gamma_k}{\partial \theta_k}}
\end{eqnarray}
Evaluated at $\theta=\TrueParam$, we get
\begin{equation}
\mathbf{h}_{jk} = \left.\frac{\partial^2 M}{\partial \theta_i \partial \theta_k}\right|_{\theta=\TrueParam} \approx 2\sum_{k=1}^{m}{\frac{\partial \gamma_k}{\partial \theta_i}\frac{\partial \gamma_k}{\partial\theta_k}}
\end{equation}
The partial derivatives $\left.\frac{\partial \gamma_k}{\partial \theta_i}\right|_{\theta}$ form an $m\times d$ matrix; in vector calculus we call
this the **Jacobian matrix** of $\gamma$, so let's write this matrix $\Jacobian(\theta)$. Then we can say
\begin{equation}
\mathbf{h}_{jk} \approx 2 \Jacobian(\TrueParam)^T \Jacobian(\TrueParam)
\end{equation}
(You should check that this gives a $d\times d$ matrix, as required.)

Turning to the variance of the gradient, we have from above that the jth component, $\nabla M(\TrueParam)_i$, is
\begin{equation}
\nabla M(G,\TrueParam)_i = 2\sum_{k=1}^{m}{(G_k-\gamma_k(\TrueParam))\left. \frac{\partial \gamma_k}{\partial \theta_i}\right|_{\theta=\TrueParam}}
\end{equation}
In matrix-vector form, this is
\begin{equation}
\nabla M(G, \TrueParam) = 2 (G - \gamma(\TrueParam)) \Jacobian(\TrueParam) = 2G \Jacobian(\TrueParam) - 2\gamma(\TrueParam)\Jacobian(\TrueParam)
\end{equation}
The non-random term doesn't matter for the variance, so
\begin{eqnarray}
\Var{\nabla M(G, \TrueParam)} & = & 4 \Var{G \Jacobian(\TrueParam)}\\
& = & 4 \Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam)
\end{eqnarray}
In reading these formulas, remember that $G = g(X_{1:n})$, so the variance on
the right-hand side changes --- we hope it shrinks! --- as $n$ grows, even though
the notation doesn't make that explicit.

Putting everything together,
\begin{eqnarray}
\Var{\hat{\theta}} & \approx &  (2 \Jacobian(\TrueParam)^T \Jacobian(\TrueParam))^{-1} 4\Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam) (2\Jacobian(\TrueParam)^T \Jacobian(\TrueParam))^{-1}\\
& = & (\Jacobian(\TrueParam)^T \Jacobian(\TrueParam))^{-1} \Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam) (\Jacobian(\TrueParam)^T \Jacobian(\TrueParam))^{-1}
\end{eqnarray}
Notice that the smaller $\Var{G}$ is, the smaller $\Var{\hat{\theta}}$ is --- less noise in the features we're matching on translates directly into more
precise estimates.

The role of the derivative matrices is harder to grasp, but we can make a _little_ more sense of this by considering the special case where
if $m=p$, so that $\Jacobian$ is a square matrix which could be inverted.  Then we have some cancellations on the right-hand side:
\begin{eqnarray}
Var{\hat{\theta}} & \approx & (\Jacobian(\TrueParam)^T \Jacobian(\TrueParam))^{-1} \Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam) (\Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam))^{-1}\\
& = & (\Jacobian(\TrueParam))^{-1} ((\Jacobian(\TrueParam)^T)^{-1}  \Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam) (\Jacobian(\TrueParam))^{-1} ((\Jacobian(\TrueParam)^T)^{-1}\\
& = & (\Jacobian(\TrueParam))^{-1} \Var{G} ((\Jacobian(\TrueParam)^T)^{-1}
\end{eqnarray}

In this situation, we can say that increasing the derivatives $\frac{\partial \gamma_k}{\partial \theta_i}$ makes $\Var{\hat{\theta}}$ smaller.  That is,
making the features more "sensitive" to the parameters we're trying to estimate
leads to more precise estimates.

(The last claim might suggest a "cheat" to get very precise estimates: just
multiply all the features through by a big number $c$.  This will, after all,
increase the derivatives by the same factor $c$, and so reduce the variance of
$\hat{\theta}$ by a factor of $c^2$.  The catch, though, is that this will also
increase $\Var{G}$ by a factor of $c^2$, so in the end it does nothing.  Can
you show that scaling all the features doesn't change $\Var{\hat{\theta}}$ in
the general case, where $\Jacobian$ isn't invertible?)

Now that you've seen how things work when there are some cancellations, it
should be easier to convince yourself that the same idea also applies when
$\Jacobian$ is non-invertible: the larger the derivatives $\frac{\partial
\gamma_k}{\partial \theta_i}$ get, the smaller $\Var{\hat{\theta}}$ gets.

#### How big is $\Var{G}$?

$\Var{G}$ is really $\Var{g(X_{1:n})}$, i.e., we're really looking at $n$
samples.  If those samples are independent of each other, then for a _lot_ of
summary statistics, $\Var{G} \propto 1/n$.  This is true of sample means, but
also (with some caveats) of medians, correlations, and many other summary
statistics.  It can even remain true if there is _some_ dependence among the
$X$'s, but not too much.  (We could make that precise but it's a different
class, namely 36-467.)

#### Practical application of the sandwich variance formula

We've written the sandwich variance formula in two ways:
\begin{equation}
\Var{\hat{\theta}} \approx \mathbf{h}^{-1}\Var{\nabla M(G,\TrueParam)} \mathbf{h}^{-1}
\end{equation}
and
\begin{equation}
\Var{\hat{\theta}} \approx (\Jacobian(\TrueParam)^T \Jacobian(\TrueParam))^{-1} \Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam) (\Jacobian(\TrueParam)^T \Var{G} \Jacobian(\TrueParam))^{-1}
\end{equation}
Both involve knowing something about derivatives at the true parameter value
$\TrueParam$, which of course is what we're trying to figure out!

To actually do calculations, what we usually do is substitute in $\hat{\theta}$
for $\TrueParam$ throughout.  This is justified on the grounds that, after all,
$\hat{\theta} \rightarrow \TrueParam$, and the hope that all the functions are
continuous in the parameter.  That still leaves the problem of getting
$\Var{G}$, the variance of the summary statistics, but there are lots of ways
of tackling this.  (You know how to handle it if the summary statistics are
means, for instance.)


#### A little worked example with the log-normal

Let's go back to our example with estimating the log-normal by matching
the sample mean and sample median.  Remember that
\begin{eqnarray}
\gamma_1(\mu, \sigma^2) & = & e^{\mu+\sigma^2/2} ~  (\text{theoretical mean})\\
\gamma_2(\mu, \sigma^2 & = & e^{\mu}  ~ (\text{theoretical median})
\end{eqnarray}
So
\begin{eqnarray}
\Jacobian(\mu, \sigma^2) & =  & \left[\begin{array}{cc} \frac{\partial \gamma_1}{d\mu} & \frac{\partial \gamma_1}{\partial \sigma^2}\\
\frac{\partial \gamma_2}{\partial \mu} & \frac{\partial \gamma_2}{\partial \sigma^2} \end{array}\right]\\
& = & \left[\begin{array}{cc} e^{\sigma^2/2} \frac{\partial e^{\mu}}{\partial \mu} & e^{\mu}\frac{\partial e^{\sigma^2/2}}{\partial \sigma^2}\\
\frac{\partial e^{\mu}}{d\mu} & \frac{\partial e^{\mu}}{\partial \sigma^2}\end{array}\right]\\
& = & \left[\begin{array}{cc} e^{\mu + \sigma^2/2} & \frac{1}{2}e^{\mu + \sigma^2/2}\\
e^{\mu} & 0\end{array}\right]
\end{eqnarray}
We need to invert this, but it's just a $2\times 2$ matrix:
\begin{eqnarray}
\Jacobian(\mu,\sigma^2)^{-1} & = & \frac{1}{-e^{\mu+\sigma^2/2}e^{\mu}/2}\left[ \begin{array}{cc} 0 & -e^{\mu+\sigma^2/2}/2 \\ - e^{\mu} & e^{\mu+\sigma^2/2}\end{array}\right]\\
& = & \left[ \begin{array}{cc} 0 & 1/e^{\mu} \\ 2/e^{\mu+\sigma^2/2} & -2/e^{\mu} \end{array}\right]
\end{eqnarray}

When we've matched on $G_1$ and $G_2$, we get
\begin{equation}
(\Jacobian)^{-1} = \left[\begin{array}{cc} 0 1/G_2 \\ 2/G_1 & 2/G_2 \end{array}\right]
\end{equation}
So we have
\begin{equation}
\Var{\hat{\theta}} \approx  \left[\begin{array}{cc} 0 1/G_2 \\ 2/G_1 & 2/G_2 \end{array}\right] \left[ \begin{array}{cc} \Var{G_1} & \Cov{G_1, G_2} \\ \Cov{G_1, G_2} & \Var{G_2} \end{array}\right] \left[\begin{array}{cc} 0 1/G_2 \\ 2/G_1 & 2/G_2 \end{array}\right]
\end{equation}

\clearpage

# Fitting by Maximum Likelihood

Let's recall how the **method of maximum likelihood** works in general.

We have a sample $X_{1:n}$.  We also have a model which gives us a
probability density for such data sets, and this model contains one or more
adjustable parameters, which generically we'll write as $\theta$.  The **likelihood** of a particular data set $x_{1:n} = (x_1, x_2, \ldots x_n)$ at a particular parameter value $\theta$
is the joint probability density of that data under that parameter value,
\begin{equation}
f(x_{1:n}; \theta)
\end{equation}
If the data points are independent and identically distributed (which we usually
try very hard to arrange!), then
\begin{equation}
f(x_{1:n}; \theta) = \prod_{i=1}^{n}{f(x_i; \theta)}
\end{equation}
The **method of maximum likelihood** is to select the parameter value which
gives the data the highest possible likelihood.

People usually work with the log-likelihood rather than the likelihood, for
several reasons.  The first reason is that the $\theta$ which maximizes one is
always the same as the $\theta$ which maximizes the other (because $\log$
itself is a monotonically-increasing function); the second reason is that it
leads to easier calculations (as we'll see in a moment or two); and the third
is that it leads to easier proofs (as in the sketch in the optional section
below).  So
\begin{equation}
L(\theta) \equiv \sum_{i=1}^{n}{\log{f(x_i; \theta)}}
\end{equation}
The **maximum likelihood estimator** is the value of $\theta$ which
maximizes $L(theta)$,
\begin{equation}
\hat{\theta} \equiv \argmax_{\theta}{L(\theta)}
\end{equation}
(Notice that subscript here indicates what variable or variables we're
maximizing over: here it's just the parameter(s); we treat the data as given,
fixed.)

There are, broadly speaking, two ways to actually find $\hat{\theta}$.  One is
to say that $L(\theta)$ is just a function, which we can calculate at any
particular value of $\theta$, and hand the problem over to some general-purpose
numerical optimization algorithm.  This is _often_ what we end up doing in
practice, particularly when the models become more complicated.

The other approach is to use our calculus, and remember that generally a
maximum is a point where the first derivative is zero and the second derivative
is negative.  (The exceptions are typically when the maximum is at the boundary
of the allowed parameter region.)  Taking the first derivatives and setting
them to zero would give a system of equations
\begin{equation}
\left.\frac{partial L}{\partial \theta_i}\right|_{\theta=\hat{\theta}} = 0
\end{equation}
for each $j \in 1:d$.  This becomes
\begin{equation}
\sum_{i=1}^{n}{\frac{\partial}{\partial \theta_i}\log{f(x_i;\hat{\theta})}} = 0
\end{equation}
which, using the chain rule from calculus, can also be written
\begin{equation}
\sum_{i=1}^{n}{\frac{1}{f(x_i;\hat{\theta})}\frac{\partial}{\partial \theta_i}f(x_i;\hat{\theta})} = 0
\end{equation}
This system of equations are called the **estimating equations** for the model,
or sometimes the **first-order conditions** (because they involve only first
derivatives).  The system of estimating equations will have one equation for
each unknown parameter $\theta_i$.  Solving these equations gives us a
candidate[^checkformax] value for $\hat{\theta}$.  Sometimes, especially for
very well-behaved models, the equations can be solved explicitly, and we end up
writing down nice formulas for $\hat{\theta}$ in terms of the data.  Even when
there aren't closed-form solutions, however, there are techniques to solve the
equations numerically.


[^checkformax]: In principle, we should still check that we've found a maximum rather than a minimum or an inflection point, but it's honestly a little unusual for that to be a serious problem.

(Notice, by the way, that it'd be much harder to do the equivalent calculus if
we worked with the likelihood rather than the log-likelihood, because
differentiating a product is more complicated than differentiating a sum.  This
is one of the promised reasons for preferring _log_-likelihood.)


## Why should we use maximum likelihood estimates?

There are three big reasons.

1. Generally, maximum likelihood estimates are consistent, which remember means that as we get more and more data ($n\rightarrow\infty$), they will converge on the true value of the parameter[^pac].  There are some exceptional cases where maximum likelihood estimates don't converge on the truth, but they're genuinely odd.
2. Generally, maximum likelihood estimates are **efficient**.  We want our estimates to converge on the true value, and we can measure how far off they are by, say, the expected squared error, $\Expect{|\hat{\theta} - \TrueParam|^2}$.  If we consider any other consistent (probably-approximately-correct) estimator, say $\tilde{\theta}$, it's generally the case that $\Expect{|\tilde{\theta} - \TrueParam|^2} \geq \Expect{|\hat{\theta} - \TrueParam|^2}$.  Other estimators _usually_ make bigger errors[^cramerrao].
3. We can be precise about the standard errors and confidence sets of the MLE, at least asymptotically.

The next sub-section, which is optional, sketches the argument for
consistency.  The one after that, which is less optional, states the
results about standard errors and confidence sets.

[^pac]: The formal definition is this: we can pick any margin of error $\epsilon$, and any error probability $\delta$, and if $n$ is big enough, $\Prob{|\hat{\theta} - \theta| \epsilon} \leq \delta$.  We also say that $\hat{\theta}$ **converges in probability** to $\TrueParam$.  In the 1980s, computer scientists re-invented this concept, which statisticians had been using explicitly since the 1920s, under the name of "probably approximately correct".  This was a clear case of re-inventing the wheel, but the re-invention did have a much more descriptive and transparent name, which statisticians should probably adopt.

[^cramerrao]: If this doesn't sound like a very crisply stated mathematical theorem, that's because it isn't.  The actual result, the Cramer-Rao theorem, is fairly complicated even to state, but what I say in the main text does get the main idea across.  If you want the gory details, see, e.g., @Pitman-basic-theory.

### Why the MLE is consistent (optional and sketchy)

Let's define
\begin{equation}
\overline{L}(\theta) = \frac{1}{n}L(\theta) = \frac{1}{n}\sum_{i=1}^{n}{\log{f(x_i;\theta)}}
\end{equation}
The factor of $1/n$ doesn't change the location of the maximum, so
\begin{equation}
\hat{\theta} = \argmax_{\theta}{\overline{L}(\theta)}
\end{equation}
The point of introducing $\overline{L}$ is that it's a sample average of a function of the data, and we know something about the behavior of sample
averages.  Since we're supposing the data points $X_i$ are independent and identically distributed[^iid], the terms $\log{f(X_i;\theta)}$ are  also
IID.  Now the law of large numbers tells us that, with IID samples, for any function $h$,
\begin{equation}
\frac{1}{n}\sum_{i=1}^{n}{h(X_i)} \rightarrow \Expect{h(X)} = \int{h(x) f(x; \TrueParam) dx}
\end{equation}
That is, sample averages converge on expected values.  So
\begin{equation}
\overline{L}(\theta) \rightarrow \int{\log{(f(x;\theta))} f(x;\TrueParam) dx} \equiv \lambda(\theta, \TrueParam)
\end{equation}
The function $\lambda(\theta, \TrueParam)$ says what _expected_ log-likelihood
the parameter value $\theta$ will get when the true parameter value is
$\TrueParam$.

[^iid]: Maximum likelihood also works if the data points _aren't_ IID, but we need to find some replacement or generalization for the law of large numbers that works with whatever kind of dependence or heterogeneity we're dealing with.

Since the random function $\overline{L}(\theta)$ is approaching the deterministic
function $\lambda(\theta, \TrueParam)$, the MLE $\hat{\theta}$ will approach
\begin{equation}
\argmax_{\theta}{\lambda(\theta, \TrueParam)}
\end{equation}
We can say something about what that maximizer is by using a fundamental (if hardly obvious) result from probability theory, sometimes called **Gibbs's inequality**, which says that for any two pdfs $f$ and $f^{\prime}$,
\begin{equation}
\int{\log{f^{\prime}(x)} f(x)} \leq \int{\log{f(x)} f(x) dx}
\end{equation}
and that the only way to get equality is if $f^{\prime}(x) = f(x)$ everywhere[^ae] (Complementary Problem 1).  Applied here, we see that the $\theta$ which maximizes $\lambda(\theta, \TrueParam)$ is in fact  nothing other than $\TrueParam$.

[^ae]: More exactly, $f^{\prime}(x) = f(x)$ "almost everywhere", i.e., except on a set of points of total length ("measure") zero.  If you know enough to demand this qualifier, though, you also know enough to insert it for yourself.

Conclusion: the MLE $\hat{\theta}$ should converge on the true parameter value
$\TrueParam$ as $n\rightarrow\infty$.

#### Going deeper into the weeds: what if the model is wrong?

In the previous argument, there was only one place where I assumed that the
probability model was correct.  This is when I asserted that the true pdf
from which the data are drawn is $f(x;\TrueParam)$ for the true-but-unknown
$\TrueParam$.  Suppose the _actual_ pdf was instead some $p(x)$, and this
pdf is not equal to $f(x;\theta)$ no matter what value of $\theta$ we try.  In
the jargon, we say that the model is then **mis-specified**.  What
happens to the MLE when the model is mis-specified?

Well, the law of large numbers part of the argument still holds, so
\begin{equation}
\overline{L}(\theta) \rightarrow \Expect{\log{f(X;\theta)}}
\end{equation}
but now
\begin{equation}
\Expect{\log{f(X;\theta)}} =  = \int{\log{f(x;\theta)} p(x) dx}
\end{equation}
$\hat{\theta}$ will converge to the $\theta^*$ which maximizes _this_.  What
can we say about it?

\begin{eqnarray}
\int{\log{f(x;\theta)} p(x) dx} & = & \int{\log{p(x)} p(x) dx} + \int{\log{\left( \frac{f(x;\theta)}{p(x)}\right)} p(x) dx} \leq \int{\log{p(x)} p(x) dx}
\end{eqnarray}
This tells us that $\int{\log{\left( \frac{f(x;\theta)}{p(x)}\right)} p(x) dx} \leq 0$.  It's more common to work with the negative of this quantity,
\begin{equation}
D(p \| f(\cdot;\theta)) = \int{\log{\left( \frac{p(x)}{f(x;\theta)}\right)} p(x) dx}
\end{equation}
is called the **divergence**[^divergence] of the density $f(\cdot;\theta)$ from
$p$.  It's the expected log likelihood ratio between the true distribution
($p$) and the one our model hypothesizes ($f(\cdot; \theta)$).  The divergence
is $\geq 0$, and it's only $=0$ if the two distributions coincide completely.

[^divergence]: Other names for this quantity include the **Kullback-Leibler divergence** (after @Kullback-Leibler), the KL divergence, the KL **information** (or "KL information number" or "KL information criterion"), and the **relative entropy**.

In general, then, when the model is not right, the MLE will converge on the
parameter value $\theta^*$ which minimizes the divergence from the true
distribution_.  It will, in this very particular sense, find the best
approximation to the true distribution that the model can provide.
Econometricians have a wonderful name for $\theta^*$; they call it the
"pseudo-truth", or the "pseudo-true parameter value".  So when the model
is correctly specified, the MLE converges on the truth, and when the model
is mis-specified, the MLE converges on the pseudo-truth.


### Uncertainty around the MLE (less optional)

Since the MLE is the maximum of the log-likelihood, the first derivative
there is zero, and the second derivative tells us how sharply peaked the
log-likelihood function is.  It should make some sense that the more sharply
peaked the the log-likelihood function, the more easily we can locate
the maximum, and the less uncertainty there is.  This not only makes
sense, but is true.  The **observed information** is
\begin{equation}
I(\hat{\theta}) \equiv -\frac{d^2 L}{d\theta^2}(\hat{\theta})
\end{equation}
and, quite generally,
\begin{equation}
\Var{\hat{\theta}} \approx I(\hat{\theta})^{-1}
\end{equation}
when the model is well-specified.

Notice that $L$ is a sum of $n$ terms, so, with independent samples, we can
expect it to be $\propto n$.  Hence $I(\hat{\theta})^{-1} \propto n^{-1}$.  The
standard error of the MLE should therefore be $\propto n^{-1/2}$.

If we're not sure that the model is well-specified, that the MLE is well-behaved, etc., a robust, if computationally more expensive, alternative, is to
use the bootstrap (lecture 6).

#### What if there's more than one parameter?

We need the $d\times d$ matrix of all partial derivatives,
\begin{equation}
I_{ij}(\hat{\theta}) \equiv - \frac{\partial^2 L}{\partial \theta_i \theta_j}
\end{equation}
and then
\begin{equation}
\Var{\hat{\theta}} \approx \mathbf{I}(\hat{\theta})^{-1}
\end{equation}

The matrix of second partial derivatives of a function is called that
function's **Hessian**, so you will often encounter that word in statistical
contexts, R documentation, etc.

The argument that this should be $\propto n^{-1}$, and standard errors $\propto
n^{-1/2}$, works just as though there was one parameter.

#### Expected information, Fisher information

The **Fisher information** in one observation at parameter value $\theta$ is
a matrix, with entries
\begin{equation}
F_{ij}(\theta) \equiv \Expect{\frac{\partial^2 \log{f(X; \theta)}}{\partial \theta_i \partial \theta_j}}
\end{equation}
The **expected information** in $n$ observations is
\begin{equation}
J_{ij}(\theta) \equiv \Expect{\frac{\partial^2 f(X_{1:n}; \theta)}{\partial \theta_i \theta_j}}
\end{equation}
With IID observations, $J_{ij}(\theta) = n F_{ij}(\theta)$.  That is,
information adds up across independent observations.  This is relevant
because when the true parameter value is $\TrueParam$,
\begin{equation}
\Var{\hat{\theta}} \rightarrow \mathbf{J}(\TrueParam) = n\mathbf{F}(\TrueParam)
\end{equation}
More exactly,
\begin{equation}
\frac{1}{n}\Var{\hat{\theta}} \rightarrow \mathbf{F}(\TrueParam)
\end{equation}

We don't know $\TrueParam$, but we do have an estimate of it, $\hat{\theta}$,
and that's generally consistent, so we can say
\begin{equation}
\frac{1}{n}\Var{\hat{\theta}} \rightarrow \mathbf{F}(\hat{\theta})
\end{equation}

Taking the expectation buried in the definition of $\mathbf{F}$ can be hard,
but we usually don't have to.  Again assuming IID observations, the law of
large numbers applies: at any fixed $\theta$,
\begin{equation}
\frac{1}{n}\mathbf{I}(\theta) \rightarrow \mathbf{F}(\theta)
\end{equation}
since $F\mathbf{F}(\theta)$ is the expected value of the observed information
$\mathbf{I}$ in a single data point.  Since the MLE is consistent, and
$\mathbf{F}$ is (usually) continuous in $\theta$,
\begin{equation}
\frac{1}{n}\mathbf{I}(\hat{\theta}) \rightarrow \mathbf{F}(\hat{\theta}) \rightarrow \mathbf{F}(\TrueParam)
\end{equation}
So we're back to saying that
\begin{equation}
\Var{\hat{\theta}} \approx \mathbf{I}(\hat{\theta})^{-1}
\end{equation}


### Situations where maximum likelihood gives unreasonable or inconsistent answers

Maximum likelihood _usually_ works, but it doesn't always.  Some of the places
where it fails are carefully-arranged mathematical horror shows.  Others are
more logical problems.  For instance, if we have a collection of $n$ data
points, $x_1, x_2, \ldots x_n$, and ask "what pdf, out of all possible pdfs,
gives these points the highest likelihood?", we'll find that the answer is
always a set of infinitely high, infinitely narrow spikes centered at the data
points.  (In symbols, $f(x) = \frac{1}{n}\sum_{i=1}^{n}{\delta(x-x_i)}$.)  (If
we try using merely very tall, very narrow spikes centered at the $x_i$, we can
always increase the likelihood by making them taller and narrower.)  This isn't
_wrong_, exactly --- any probability we put on a point $x$ which isn't one of
the $x_i$ has to be taken _from_ the $x_i$ --- but it seems unreasonable.
Maximum likelihood generally does poorly when it faces such
"infinite-dimensional" or "nonparametric" problems.

Finally, I should mention that the MLE also does poorly for finding a parameter
like $\xmin$ in the Pareto distribution --- it always suggests using the
smallest observed $x$.  This is actually consistent _if_ the data come purely
from a Pareto distribution.  But usually we just want a Pareto _tail_, and are
being tactically (or tactfully) vague about the distribution of the body.
Since we're not including the distribution of the body in our calculation of
the likelihood, we shouldn't be surprised that the MLE does something
unreasonable.  (Once we have $\xmin$, though, using the MLE to get $\alpha$
is definitely the right way to go.)  For a more reasonable approach to
estimating $\xmin$, see @power-law-or-not.

\clearpage

# Matching Features vs. Maximizing Likelihood

I have gone on at some length about how to use maximum likelihood because it's
a fundamental technique in statistics, and usually the best thing to do _when
we can_.  The problem, for a lot of social analyses, is with the italicized
phrase.

Maximum likelihood presumes not just that _someone_ collected a sample of
individual data points, but that we have access to those data points, so we can
evaluate $\log{f(x_i;\theta)}$ and add it up across the data set.  Often,
however, what _we_ have access to aren't the individual $x_i$, but _only_
summary statistics  $g_1(x_1, \ldots x_n)$, $g_2(x_1, \ldots
x_n)$, down to $g_m(x_1, \ldots x_n)$.  These are functions of the data, as the
notation indicates, but they are _not_ the data themselves.  We then
need to fall back on the technique of matching these summaries[^mleforsummaries].

[^mleforsummaries]: We could in principle calculate the joint distribution of all the summary statistics as functions of the parameter $\theta$, and so maximize the likelihood of the summaries.  The key phrase in that sentence is "in principle" --- working out exact sampling distributions of summary statistics is usually extremely hard.  (For an example of the sort of heroic calculations required, see, e.g., @Virkar-Clauset-binned-power-laws, where the model is a power law and the summary statistics are the number of observations in fixed bins or brackets.)  For very large $n$, many sampling distributions become approximately Gaussian, which simplifies this approach [@Wood-inference-for-nonlinear-dynamic], but also makes it more similar to a least-squares matching estimate...

There are various reasons why we might need to do this.

1. _The original data is too large to release._  This can still be true, but is less of an issue today than in, say, 1900, or even 1980, and will presumably be less and less of an issue in the future.
2. _The original data is too sensitive to release._  A lot of our knowledge about income and wealth distributions comes from official records.  Part of the reason that people are willing to be (fairly) honest with the authorities is that modern governments are really powerful and intrusive and lying to them too blatantly about money has reliably bad consequences.  But part of the way modern governments soften this pill is by promising not to disclose everyone's personal information to anyone with an Internet connection.  If you are a trusted employee of the IRS or the Census Bureau, you could probably calculate MLEs on the full tax records of the US population.  But the rest of us are never going to see _that_ data.
3. _The original data was not preserved_.  Especially if we're interested in
historical comparisons, it becomes harder and harder to retrieve the original data, but summaries may be easily obtained.

Of these three reasons, the second is probably the one which will be the most
important to you in practice.



## Summary matching estimates are usually less efficient than maximum-likelihood estimates

I said above that maximum likelihood estimates are generally _efficient_,
meaning that nothing else converges on the truth faster.  Summary-matching
estimates are "something else", so they will generally converge slower than the
MLE.  This is disappointing, but we're usually using them because we don't have
access to the data which we'd need to calculate the MLE.

### Sufficient Statistics

There is a theoretically-important case where summary-matching estimates _are_
efficient.  This is when we're matching what are called **sufficient
statistics**.  This only happens when the likelihood can be written in the
following special form:
\begin{equation}
f(x_{1:n}; \theta) = t(s(x_{1:n}); \theta) \phi(x_{1:n})
\end{equation}
That is, we can factor the whole likelihood into two terms, one of which,
$\phi(x_{1:n})$, doesn't involve the parameter at all, and the other, which does involve
the parameter, doesn't use the data directly, but only a function of the data
$s(x_{1:n})$.  That function $s$ is the "sufficient statistic"[^suff].

[^suff]: That definition is based on what's called the "Neyman factorization criterion".  There are actually multiple, equivalent ways of defining "sufficient statistic".  Another way to do so is to say that a statistic is sufficient when it keeps all the information the data had about the parameter $\theta$.  This can be made precise using information theory [@Kullback-Leibler].

Most models do _not_ have sufficient statistics, but many familiar ones do.
For the Gaussian distribution, for instance, the sample mean and the sample
variance are sufficient.  For the log-normal, the sufficient statistics are the
sample mean logarithm and the sample variance of the logarithm.  For a Pareto
distribution with known minimum, the mean logarithm is a sufficient statistic.

When we can identify a model's sufficient statistics, and they're available to
us, we should certainly prefer matching them when we can.  In fact, a lot of
the time[^exponentialfamily], matching sufficient statistics gives us the MLE.

[^exponentialfamily]: The distribution needs to be (or be put into) what's called **exponential family** form, so $f(x_1, x_2, \ldots x_n; \theta) = h(x) \frac{e^{ s(x_{1:n}) \cdot \theta}}{z(\theta)}$.  (The denominator here, $z(\theta)$, makes sure that this is a valid probability distribution that adds up to 1 across all possible values of $x_{1:n}$.  It's called the **normalizing factor** or sometimes the **partition function**.)}  If a distribution has statistics and satisfies some additional technical requirements, it _must_ be an exponential family, but there are many common distributions which aren't exponential families, like the $t$ distribution, or some of the natural generalizations of the Pareto distribution.

\clearpage

# Summing Up on Fitting

- If we have individual-level data, the method of maximum likelihood is usually
better than matching summary statistics.
    + The MLE will generally converge on the
true parameter values ("consistency").
    + The MLE will converge at least as fast as any other
consistent estimator ("efficiency").
    + For large samples, the variance of the MLE can be found from the matrix of 2nd derivatives of the log-likelihood.   (Parameters which make a big difference to the log-likelihood will be more precisely estimated than the others.)
    + There's usually no exact formula for the MLE, but that's why we have numerical optimization algorithms.
- Matching summary statistics or features is _usually_ inferior to the MLE.
    + If statistics converge on population values as $n$ grows, and the features uniquely identify parameters, then matching is consistent.
    + Matching estimates are more precise when the summary statistics have less variance (the statistics we're matching aren't very noisy)
    + Matching estimates are more precise when the features are very sensitive to the parameters (small changes to parameters imply big changes to the features).
    + If we know the variances of the summary statistics, and how sensitive the features are to the parameters, we can calculate (approximate, large-$n$) variances for our matching estimators
    + Matching estimators will have at least as much variance as the MLE
- There are many situations, especially in social statistics, where we _only_ have summary statistics, so we need to rely on matching.

\clearpage

# Checking Goodness-of-Fit

We fit our model by estimating the parameter(s) $\theta$.  The model however
may or may not be right --- it's responsible to check the over-all **goodness
of fit**.  (Or, if we're pessimists, the **lack of fit**.)  There are two
big ways of doing this:

- See if the model, with the estimated parameters, can predict _other_ features of the data
- See if the model, with the estimated parameters, gives a distribution that matches the over-all distribution of the data.

If you like, you can think of _both_ of these approaches as instances of a more
general rule, "see if the model, with the estimated parameters, can predict
something beyond those parameters".  But the practicalities are somewhat
different.

## Checking by Calculating New Features

This works rather like estimating by matching features.  On the one hand, we
have some summary statistic of the data, let's say $g_{test}(X_1, \ldots X_n)$.
This, as the notation indicates, is used as the "test statistic".  On the other
hand, we can calculate what, theoretically, that should be, through some
function $\gamma_{test}(\theta)$.  We can then just _compare_ the observed
value on our actual data, $g_{test}(x_1, \ldots x_n)$, to the theoretically
predicted value with our estimated parameters, $\gamma_{test}(\hat{\theta})$.
The bigger the gap, the worse the fit.  We can do this _whether or not_ our
estimate comes from matching features, from maximizing likelihood, from a
dream, etc.
Basically _any_ feature of the data can be used as test statistic in this sense,
though not all of them will work equally well.

There are some important considerations here.

- _Do not expect a perfect match_.  In the first place, the data are random, so $G_{test}=g_{test}(X_1, \ldots X_n)$ is going to be a random variable.  Even if we somehow knew $\TrueParam$, we shouldn't expect $G_{test} = \gamma_{test}(\TrueParam)$ all the time.  In the second place, we've only got an estimate $\hat{\theta}$, not $\TrueParam$.
   + As with estimating by matching, we prefer test statistics with small variance, and we prefer features which are very sensitive to the underlying parameters.
- _Deciding how big a mis-match is too big requires more work_.  We'll look, later in the course, at some methods for working how much the test statistic can differ from the theoretical prediction when the model is right.  (One way to do this is called "bootstrapping".)  For now, though, we often prefer a more
qualitative, touchy-feely evaluation.
- _Some features can be more important than others_ We might, in some applications, really _care_ about some features more than others.  (For instance, we might really want to get the income shares of high percentiles right.)  In those situations, it can be legitimate to focus about whether those features get predicted well, and not worry so much about other features which might be missed.

### Do Not Use Features Twice

At this point however I need to issue an important warning.  If we estimated
our model by matching on some features, we cannot check the fit of the model by
seeing if it predicts those same features.  We've just made sure that it does!
At best, this just makes sure we didn't botch any of our calculations.  So, for
example, if we estimate a log-normal distribution by matching the mean and the
median, we learn nothing by seeing that the fitted model successfully predicts
the mean and the median --- what else could it do?  The same goes for any
statistic which is a function of the features we matched on (say the ratio of
mean to median).  What we need is to look at _other_ features, which
are not implied by what we matched on[^sufficiencyandchecking].

[^sufficiencyandchecking]: You might think this wouldn't be a concern if we used maximum likelihood instead of feature-matching.  However, _if_ the model as sufficiency statistics (as defined above), especially if it's an "exponential family" model (as defined above), the MLE will automatically match those statistics.  So for instance using the MLE to estimate the log-normal will automatically match the mean of the log and the variance of the log, and so anything which is a function of those (like the "coefficient of variation" of the log, $\Var{\log{X}}/\Expect{\log{X}}$).

## More Global Checks of Fit for Distributions

We can also check the fit of a model by comparing the distribution with the
estimated parameters to some _other_, less-model-dependent way of estimating
the distribution.  For instance, we might compare the CDF implied by the model
with the "empirical" CDF of the raw data.  We can divide the range of $X$ up
into bins or brackets, and compare the actual number of samples in each bin to
the number we'd expect based on the model.  This in turn is very similar to
comparing the pdf implied by the model to a histogram of the data, which
suggests we can just compare the model's pdf to other estimates of the pdf
which don't rely on the model (such as the ones plotted in the last few
lectures).

Some of these comparisons can be turned into formal statistical tests, like
the Kolmogorov-Smirnov test (which compares the empirical CDF of the data to
the CDF implied by the model) or the $\chi^2$ test (which compares the observed
counts in bins to those expected under the model).  The distribution of the
test statistic is usually computed assuming that the theoretical model has no
adjustable parameters.  This is, of course, rarely the case.  If the model
_does_ have parameters estimated from the data, then the model's fit will be
apparently more impressive (the test statistic will be smaller) than the usual
calculation suggests, _just because_ the model was adjusted to fit _this_
particular data set.  The simplest way to avoid this problem is **data
splitting**: randomly split the data into two distinct parts, fit the model on
one part and check the fit on the other.  When data splitting is impractical,
though, there are other ways to correct for parameter estimation, such as the
bootstrap methods we'll look at later in the class.



## Good Fit $\not\Rightarrow$ Truth

Even if we've found a model which seems to fit well, we should be cautious
about concluding that it's _true_.  There are at least two reasons for this.

1. Other models might also fit the data in the ways we're checking.
2. The model might fail to fit the data in other ways, which we haven't looked at.

There are two ways to deal with the first problem.  We can _either_ try to
report the whole range of models which are (so far as we can tell) compatible
with the data, and admit to a certain amount of **model uncertainty**[^confidenceset].  _Or_ we
can try to find some _other_ aspect of the data, or new data, where the
different models make _different_ predictions, and test those predictions.
(That's science.)  Or, of course, we can pursue both approaches at once, if we
have the time and the energy.

[^confidenceset]: Within a model, the set of all parameter values which are compatible with the data, at a certain level of goodness-of-fit, forms a **confidence interval**, or more generally a **confidence set**.  In words: the confidence set contains all the parameter values we _can't_ confidently reject.  This idea is just extending that same approach to all the _models_ we can't reject.  (It was advocated long ago by @Haavelmo-probability-approach, and quite correctly in my opinion.)

That even the best-fitting model might fail to accommodate new data, or more
refined probes of old data, is just a fact of life: science is fallible, or
more precisely corrigible.  There are even many situations where we continue to
use models we know are not altogether accurate, if the inaccuracies are small
and/or un-important for our purposes.











\clearpage


# Complementary Problems

1. Prove Gibbs's inequality.
2. _Feature weighting_.  Consider doing feature matching by solving a _weighted_ least squares problem,
\begin{equation}
\hat{\theta}_w = \argmin_{\theta}{\sum_{k=1}^{m}{(G_k - \gamma_k(\theta))^2 w_k}}
\end{equation}
with weights $w_k > 0$.
    a. Find a formula for the asymptotic variance of this estimator, similar to the one derived for the unweighted estimator above.  Check that when all the weights are equal, your answer matches  the formula for the unweighted estimator.
    b. Using your formula from (a), find the weights $w_k$ which make $\Var{\hat{\theta}_{w}}$ as small as possible.  Is the answer unique?  How would you describe it in words?
    c. Convince yourself that using weights $w_k$ is the same as replacing the function $g_k(x_{1:n})$ by $\sqrt{w_k}g_k(x_{1:n})$, and similarly for $\gamma_k(\theta)$.  Why then does using weights make any difference?
3. _Propagation of error_ Read the handout on "propagation of error" linked to on the class website for HW 1.
    a. Use this technique to find formulas for the standard errors of estimating the log-normal parameters from the sample mean and the sample median.
    b. Do your results from (a) formulas match the ones obtained above, from the "sandwich variance" formula?  Should they?
    c. What's the general relationship between the propagation-of-error approach and the sandwich variance?  _Hint_: If $f(x)$ has an inverse, how is the derivative of $f^{-1}$ related to the derivative of $f$?
	

\clearpage

# Further reading

Estimation theory is a huge subject, which can be pursued at many levels of
mathematical and computational sophistication.  Some of the standard references,
like @Lehmann-estimation, are pretty stiff going.

There is a good discussion of likelihood-based inference with minimal
technicalities, in @Cox-principles-of-statistical-inference, and, with a _bit_
more mathematical demands, in @Pitman-basic-theory.  The latter gives the
simplest proof I've seen of the Cramer-Rao inequality, a key result which says
that the Fisher information gives a lower bound on the variance of any
estimator.  (This result goes back to @Cramer and @Rao-introduces-Cramer-Rao.)

More abstractly, maximum likelihood and feature-matching are both instances of
a general strategy for estimation, called "$M$-estimation", based on minimizing
(or maximizing) some function which says how well the model matches the data.
I have tried to present both feature-matching and the MLE in a way which
suggests this common structure.  In doing so I've been very influenced by the
presentation in @van-der-Vaart-asymptotic-stats and @van-de-Geer-empirical, but
both of these call for a lot more math than this course.  I'm not really aware
of another presentation of _general_ $M$-estimation at our level.

On the importance of test models, and using summary statistics to do so,
see @Gelman-Shalizi-philosophy.



# References
