\documentclass[11pt]{article}
\usepackage{indentfirst}
\usepackage{amsmath}
\usepackage{natbib}
\usepackage{url}
\begin{document}
\title{Cauchy Example for Le Cam Made Simple}
\author{Charles J. Geyer}
\maketitle
\section{Introduction}
We do a very simple one-parameter model for which everything is trivial,
the Cauchy location model. This model is mildly notorious among theoreticians
because the likelihood is multimodal and the multimodality does not go
away as $n$ goes to infinity \citep{reeds}. It's not that there is any theory
that leads one to expect that multimodality should go away as $n$ goes
to infinity, but just that the Cauchy model is simple enough to analyze
and in it one can see that indeed what theory doesn't lead one to expect
actually doesn't happen.
So let's make up some data. Since the model is translation invariant,
it doesn't matter what we choose at the ``simulation truth'' parameter value
<>=
set.seed(42)
n <- 30
theta0 <- 0
x <- theta0 + rcauchy(n)
@
To do maximum likelihood using local optimization we need a good starting
point for Newton's method. The obvious starting point for the Cauchy
location model is the sample median, easily calculated, highly robust,
guaranteed $\sqrt{n}$-consistent.
Minus the log likelihood (minus because \verb@nlm@ does minimization rather
than maximization)
<>=
mlogl <- function(theta, x) sum(- dcauchy(x, theta, log = TRUE))
@
and the MLE
<>=
out <- nlm(mlogl, median(x), typsize = n, hessian = TRUE, x = x)
print(out)
@
So our MLE, observed Fisher information (we will pretend we
don't know how to calculate expected Fisher information), and asymptotic
confidence interval calculated from them are
<>=
theta.hat <- out$estimate
info <- out$hessian
conf.level <- 0.95
crit <- qnorm((1 + conf.level) / 2)
theta.hat + crit * c(-1, 1) / sqrt(info)
@
Couldn't be simpler.
\section{Single Bootstrap}
By ``bootstrap'' we mean \emph{parametric bootstrap}, since we are working
with maximum likelihood.
<>=
nboot <- 200
save.seed <- .Random.seed
theta.star <- double(nboot)
theta.star.code <- double(nboot)
info.star <- double(nboot)
for (iboot in 1:nboot) {
x.star <- theta.hat + rcauchy(n)
out <- nlm(mlogl, median(x.star), typsize = n, hessian = TRUE, x = x.star)
theta.star[iboot] <- out$estimate
theta.star.code[iboot] <- out$code
info.star[iboot] <- out$hessian
}
all(theta.star.code <= 3)
@
The pivotal quantity used to construct the confidence interval
(or, more precisely, its bootstrap approximation) is
<>=
z.star <- (theta.star - theta.hat) * sqrt(info.star)
@
Figure~\ref{fig:one} shows its Q-Q plot
<