---
title: Density estimation
author:
date: 21 March 2019
output: slidy_presentation
---

```{r, echo=FALSE}
library(knitr)
opts_chunk$set(size="small",background="white", highlight=FALSE, cache=TRUE,
               autodep=TRUE, tidy=TRUE,
               warning=FALSE, message=FALSE)
options(np.messages=FALSE)
```

\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Indicator}[1]{\mathbb{1}\left( #1 \right)}
\]

## Data: body weights of cats

```{r}
library(MASS)
data(cats)
head(cats$Bwt)
summary(cats$Bwt)
```

## The empirical cumulative distribution function

- ECDF:
\[
\hat{F}_n(a) = \frac{1}{n}\sum_{i=1}^{n}{\Indicator{x_i \leq a}}
\]
    + Makes a jump at each distinct value seen in the data
	+ Constant in between the jumps


## Computationally: `ecdf`

```{r}
plot(ecdf(cats$Bwt), xlab="Body weight (kg)", ylab="Cumulative probability")
```

## The ECDF estimates the CDF


- Suppose $x_1, x_2, \ldots x_n$ are IID with CDF $F$.  Then for _any_ $a$,
\[
\hat{F}_n(a) \rightarrow F(a)
\]
by law of large numbers
    + Also: see backup slides

## The CDF is not enough

```{r, echo=FALSE}
plot(ecdf(cats$Bwt), xlab="Body weight (kg)", ylab="Cumulative probability")
```

- What's the probability of a cat weighing between $3.21$ and $3.22$ kg?
    + ECDF says: exactly 0

- Q from discussion: What about running a smoother (e.g., a spline) through the ECDF?
- A: You could do that, but you'd need to make sure the curve was always monotone, between 0 and 1, etc.
    + Once we estimate the density, we can recover the CDF by integrating

## The CDF is not enough

- More generally, pdf is derivative of CDF, $f(x) = dF/dx$
- Derivative of an _empirical_ CDF is either 0 or $\infty$
    + Where is it $\infty$?
- How do we estimate a pdf? 



## What about a histogram?

```{r}
hist(cats$Bwt, xlab="Body weight (kg)", ylab="Number of cats", breaks=20)
```

- Gives us an idea of where we have more or fewer data points...


## Fitting a parametric model

- Pick a parametric pdf that looks sensible
- For example, the gamma distribution with shape parameter $a$ and scale $s$ is
\[
f(x) = \frac{1}{s^a \Gamma(a)} x^{a-1} e^{-x/s}
\]
    + Exponential is the special case $a=1$
    + $a > 1$ is a hump with an exponential right tail

```{r, echo=FALSE}
curve(dgamma(x, shape=1, scale=1), lty="solid", from=0, to=4,
      ylab="gamma pdf")
curve(dgamma(x, shape=2, scale=1), lty="dashed", add=TRUE)
curve(dgamma(x, shape=2, scale=0.5), lty="dotted", add=TRUE)
legend("topright", legend=c(expression(list(a==1, s==1)),
                            expression(list(a==2, s==1)),
                            expression(list(a==2, s==0.5))),
       lty=c("solid", "dashed", "dotted"))
```

## Fitting a parametric model (cont'd)

- Now maximize the likelihood 
- For some (but not all) standard distributions, the function `fitdistr` in the package MASS
  has pre-programmed maximum likelihood estimates
    + You saw this in homework 3 for the $t$ distribution
	+ As it happens, `fitdistr` knows about gamma distributions too
- It's also easy to maximize the log-likelihood yourself

```{r}
# Define the (negative, normalized) log-likelihood function
gamma.loglike <- function(params, data=cats$Bwt) {
    -mean(log(dgamma(x=data, shape=params[1], scale=params[2])))
}
# Ask R to optimize it, with a starting guess about the parameters
cats.gamma <- optim(fn=gamma.loglike, par=c(1,1))
# What is the best (negative, normalized) log-likelihood?
cats.gamma$value
# What parameters did optim find for us?
cats.gamma$par
```

## Compare the fitted density to the histogram

```{r}
plot(hist(cats$Bwt,plot=FALSE, breaks=20), freq=FALSE,
     xlab="Body weight (kg)")
curve(dgamma(x, shape=cats.gamma$par[1], scale=cats.gamma$par[2]), add=TRUE)
```

## The histogram is already a  density estimate

- Estimated density is _constant_ on bins
- $N_i =$ number of samples in bin $i$
- $i(x)=$ which bin $x$ falls into
- $h =$ width of the bins

\[
\hat{f}_{hist}(x) \equiv \frac{1}{h} \frac{N_{i(x)}}{n}
\]

## Why this works

Probability of falling into a region of length $h$ around $x$ will be
\[
\int_{x-h/2}^{x+h/2}{f(u) du} \approx f(x) h
\]
for very small $h$.  So
\[
N_i \sim \mathrm{Binom}(n, f(x) h)
\]
and
\[
\hat{f}_{hist}(x) = \frac{N_i}{nh} \sim \frac{1}{nh} \mathrm{Binom}(n, f(x) h)
\]
which $\rightarrow f(x)$ as $n \rightarrow \infty$

## How bad is a histogram?

Error from variance: $\propto 1/nh$

Error from bias: $\propto h^2 (f^{\prime}(x))^2$

Total error for estimating the density, $O(h^2) + O(1/nh)$

Best $h = O(n^{-1/3})$

Best total error, $O(n^{-2/3})$

vs. $O(n^{-1})$ if we've got a correct parametric model

Q from discussion: What about variable-width bins?

A: That can help, but it's even better to get away from fixed bins altogether

## Kernels

Idea: "Smooth out" each observation into a little distribution

Kernel functions _are_ pdfs

Put a copy of the kernel at each observation; add them up

\[
\hat{f}_{kern}(x) = \frac{1}{n} \sum_{i=1}^{n}{\frac{1}{h} K\left(\frac{x-x_i}{h}\right)}
\]

## How this works for cats

Say $h=0.01$

```{r, results="hide", echo=FALSE}
# Artificially add a very small amount of noise to each value, just for
# plotting purposes
jittered.cats <- jitter(cats$Bwt)
plot(0, xlim=range(cats$Bwt), ylim=c(0,4), type="n",
     xlab="Body weight (kg)", ylab="Density") # Establish empty plot
# Tick mark on the axis for each data point
rug(jittered.cats, side=1, col="grey")
# Put a little, scaled-down copy of the Gaussian pdf around each data point
sapply(jittered.cats, function(xi) {curve(dnorm(x, mean=xi, sd=0.01)/nrow(cats), n=1000, add=TRUE, col="grey")})
# Show what we get when we sum up all the little grey curves
lines(density(cats$Bwt, bw=0.01))
```


## How do we pick the bandwidth?

- Cross-validation
- Cross-validate _what_?
-  **Integrated squared error**: $\int{(f(x) - \hat{f}(x))^2 dx}$
- Or (negative, normalized) **log-likelihood** $\int{-\log{(\hat{f}(x))} f(x) dx}$

## Evaluating the log-likelihood

- The integral is an expected value:
\[
\int{-\log{(\hat{f}(x))} f(x) dx} = \Expect{-\log{\hat{f}(X)}}
\]
- With $n$ IID samples $x_i$,
\[
\Expect{-\log{\hat{f}(X)}} \approx \frac{1}{n}\sum_{i=1}^{n}{-\log{\hat{f}(x_i)}}
\]
by LLN
    + A similar but more complicated trick for ISE; see backup slides



## Computationally: `npudens`

```{r}
library(np)
bwt.np <- npudens(cats$Bwt)
plot(bwt.np, xlab="Body weight (kg)")
```

## Joint PDFs

- What about multiple variables?  We want $f(x_1, x_2, \ldots x_p) = f(\vec{x})$

```{r}
plot(Hwt ~ Bwt, data=cats, xlab="Body weight (kg)", ylab="Heart weight (g)")
```

- Parametric models work just the same way as before
- Kernel density estimates work almost the same way


## Joint KDEs

- Use a product of kernels, one for each variable
\[
\hat{f}_{kern}(x_1, x_2) = \frac{1}{n} \sum_{i=1}^{n}{\frac{1}{h_1 h_2} K\left(\frac{x_{1}-x_{i1}}{h_1}\right)K\left(\frac{x_{2}-x_{i2}}{h_2}\right)}
\]
- Need as many bandwidths as there are variables
- Pick bandwidths for best prediction of _both_ variables

## Computationally: `npudens` again

```{r}
cats.np <- npudens(~Bwt+Hwt, data=cats)
plot(cats.np, view="fixed", theta=10, phi=45)
```

**EXERCISE**: Where do the ridges come from?

- Q1: What happened to the bandwidth on bodyweight (`Bwt`)?
    + It got a lot smaller
- Q2: Why did that bandwidth change?
   + There's a strong relationship between `Bwt` and `Hwt`; if we want to get _both_ right, we can't average over as wide a range of `Bwt`


## Conditional PDFs

- Conditional pdf:
\[
f(x_2|x_1) = \frac{f(x_1, x_2)}{f(x_1)}
\]
so find the joint and marginal pdfs and divide
- Conditional KDE:
\[
\hat{f}_{kern}(x_2|x_1) = \frac{\frac{1}{n} \sum_{i=1}^{n}{\frac{1}{h_1 h_2} K\left(\frac{x_{1}-x_{i1}}{h_1}\right)K\left(\frac{x_{2}-x_{i2}}{h_2}\right)}}{\frac{1}{n} \sum_{i=1}^{n}{\frac{1}{h_1} K\left(\frac{x-x_i}{h_1}\right)}}
\]
- Pick bandwidths for best prediction of $x_2$
    + i.e., look at the conditional log-likelihood of $x_2$ given $x_1$



## Computationally: `npcdens`

```{r}
cats.hwt.on.bwt <- npcdens(Hwt ~ Bwt, data=cats)
plot(cats.hwt.on.bwt, view="fixed", theta=10, phi=65)
```


## How well does kernel density estimation work?

- The same sort of analysis we did for kernel regression will work here too
    + More book-keeping --- see the text
- Integrated squared error of KDE is $O(n^{-4/5})$ in 1D
    + Worse than a well-specified parametric model, $O(n^{-1})$
	+ Better than a mis-specified parametric model, $O(1)$
	+ Better than a histogram, $O(n^{-2/3})$
- In $p$ dimensions, ISE of joint KDE is $O(n^{-4/(4+p)})$
    + Nothing which is _universally_ consistent does any better
    + Curse of dimensionality again


## Summing up

- It's easy to learn the CDF
    + use the empirical CDF
- It's easy to estimate parametric distributions
    + maximize the log-likelihood
- Histograms are a first step to a _non-parametric_ density estimate
    + Kind of ugly and statistically inefficient, though
- Kernel density estimation lets us learn marginal pdfs, joint pdfs, conditional pdfs...
    + The curse of dimensionality is always with us though

READING: Chapter 14




## Backup: More about the ECDF

- **Glivenko-Cantelli** theorem: If $X$'s are IID with CDF $F$, then
\[
\max_{a}{|\hat{F}_n(a) - F(a)|} \rightarrow 0
\]
    + This is a stronger result than just saying "for _each_ $a$, $|\hat{F}_n(a)-F(a)| \rightarrow 0$"
- **Kolmogorov-Smirnov** (KS) test: To test whether $X$'s are IID with CDF $F_0$, use $d_{KS} = \max_{a}{|\hat{F}_n(a)-F_0(a)|}$
    + If $X \sim F_0$, then $d_{KS} \rightarrow 0$ as $n\rightarrow \infty$
	+ If $X \not\sim F_0$, then $d_{KS} \rightarrow$ some positive number
    + The sampling distribution of $d_{KS}$ under the null turns out to be the _same_ for all $F_0$ (at least for big $n$)
	+ Needs modification if you estimate the $F_0$ from data as well


## Backup: Parametric density estimates

- We have a parametric model for the density $f(x;\theta)$ with a _finite_-dimensional vector of parameters $\theta$
- We maximize the log-likelihood to get $\hat{\theta}$
- Equivalently we minimize the negative normalized log-likelihood
\[
L(\theta) = -\frac{1}{n}\sum_{i=1}^{n}{\log{f(X_i;\theta)}}
\]
- This is a sample mean so LLN says $L(\theta) \rightarrow \Expect{L(\theta)} \equiv \lambda(\theta)$
- You can show that _if_ the data came from the model with some true parameter $\theta_0$, then $\lambda(\theta_0)$ is the unique global minimum

## Backup: More about parametric density estimates

- Quite generally, the MLE $\hat{\theta} \rightarrow \theta_0$
- Almost as generally, for large $n$, $\hat{\theta} \sim \mathcal{N}(\theta_0, \mathbf{J}(\theta_0)/n)$
    + Hence $\hat{\theta} = \theta_0 + O(1/\sqrt{n})$
    + Here $\mathbf{J}(\theta) =$ matrix of 2nd partial derivatives of $\lambda$ = **Fisher information**
	+ Can plug in $\mathbf{J}(\hat{\theta})$ to get approximate confidence intervals, etc.
	+ Doesn't matter if the _data_ is Gaussian --- $L(\theta)$ is a sample average so by the CLT _it's_ (approximately)  Gaussian, and that propagates
- Things are more complicated if the model is mis-specified

## Backup: Evaluating ISE

\begin{eqnarray}
\int{(f(x) - \hat{f}(x))^2 dx} & = & \int{f^2(x) - 2f(x)\hat{f}(x) \hat{f}^2(x) dx}\\
& = & \int{f^2(x) dx} - 2\int{\hat{f}(x) f(x) dx} + \int{\hat{f}^2(x) dx}\\
& = & (A) + (B) + (C)\\
& = & (\text{same for all estimates}) - (\text{depends on truth and estimate}) + \text{(depends only on estimate)}
\end{eqnarray}

- If we're comparing different estimates, term (A) doesn't matter, so we can ignore it
- Term (C) depends only on our estimate, and we can get it by (say) numerical integration
- That leaves term (B), and with $n$ independent samples $x_i$,
\[
\int{\hat{f}(x) f(x) dx} \approx \frac{1}{n}\sum_{i=1}^{n}{\hat{f}(x_i)}
\]
by LLN

