


\documentclass{chapman}

%%% copy Sweave.sty definitions

%%% keeps `sweave' from adding `\usepackage{Sweave}': DO NOT REMOVE
%\usepackage{Sweave}


\RequirePackage[T1]{fontenc}
\RequirePackage{graphicx,ae,fancyvrb}
\IfFileExists{upquote.sty}{\RequirePackage{upquote}}{}
\usepackage{relsize}

\DefineVerbatimEnvironment{Sinput}{Verbatim}{}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontfamily=courier,
                                              fontshape=it,
                                              fontsize=\relsize{-1}}
\DefineVerbatimEnvironment{Scode}{Verbatim}{}
\newenvironment{Schunk}{}{}

%%% environment for raw output
\newcommand{\SchunkRaw}{\renewenvironment{Schunk}{}{}
    \DefineVerbatimEnvironment{Soutput}{Verbatim}{fontfamily=courier,
                                                  fontshape=it,
                                                  fontsize=\small}
    \rawSinput
}

%%% environment for labeled output
\newcommand{\nextcaption}{}
\newcommand{\SchunkLabel}{
  \renewenvironment{Schunk}{\begin{figure}[ht] }{\caption{\nextcaption}
  \end{figure} }
  \DefineVerbatimEnvironment{Sinput}{Verbatim}{frame = topline}
  \DefineVerbatimEnvironment{Soutput}{Verbatim}{frame = bottomline,
                                                samepage = true,
                                                fontfamily=courier,
                                                fontshape=it,
                                                fontsize=\relsize{-1}}
}


%%% S code with line numbers
\DefineVerbatimEnvironment{Sinput}
{Verbatim}
{
%% numbers=left
}

\newcommand{\numberSinput}{
    \DefineVerbatimEnvironment{Sinput}{Verbatim}{numbers=left}
}
\newcommand{\rawSinput}{
    \DefineVerbatimEnvironment{Sinput}{Verbatim}{}
}


%%% R / System symbols
\newcommand{\R}{\textsf{R}}
\newcommand{\rR}{{R}}
\renewcommand{\S}{\textsf{S}}
\newcommand{\SPLUS}{\textsf{S-PLUS}}
\newcommand{\rSPLUS}{{S-PLUS}}
\newcommand{\SPSS}{\textsf{SPSS}}
\newcommand{\EXCEL}{\textsf{Excel}}
\newcommand{\ACCESS}{\textsf{Access}}
\newcommand{\SQL}{\textsf{SQL}}
%%\newcommand{\Rpackage}[1]{\hbox{\rm\textit{#1}}}
%%\newcommand{\Robject}[1]{\hbox{\rm\texttt{#1}}}
%%\newcommand{\Rclass}[1]{\hbox{\rm\textit{#1}}}
%%\newcommand{\Rcmd}[1]{\hbox{\rm\texttt{#1}}}
\newcommand{\Rpackage}[1]{\index{#1 package@\textit{#1} package}\textit{#1}}
\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\Rclass}[1]{\index{#1 class@\textit{#1} class}\textit{#1}}
\newcommand{\Rcmd}[1]{\index{#1 function@\texttt{#1} function}\texttt{#1}}
\newcommand{\Roperator}[1]{\texttt{#1}}
\newcommand{\Rarg}[1]{\texttt{#1}}
\newcommand{\Rlevel}[1]{\texttt{#1}}


%%% other symbols
\newcommand{\file}[1]{\hbox{\rm\texttt{#1}}}
%%\newcommand{\stress}[1]{\index{#1}\textit{#1}}
\newcommand{\stress}[1]{\textit{#1}}
\newcommand{\booktitle}[1]{`#1'} %%'

%%% Math symbols
\newcommand{\E}{\mathsf{E}}
\newcommand{\Var}{\mathsf{Var}}
\newcommand{\Cov}{\mathsf{Cov}}
\newcommand{\Cor}{\mathsf{Cor}}
\newcommand{\x}{\mathbf{x}}
\newcommand{\y}{\mathbf{y}}
\renewcommand{\a}{\mathbf{a}}
\newcommand{\W}{\mathbf{W}}
\newcommand{\C}{\mathbf{C}}
\renewcommand{\H}{\mathbf{H}}
\newcommand{\X}{\mathbf{X}}
\newcommand{\B}{\mathbf{B}}
\newcommand{\V}{\mathbf{V}}
\newcommand{\I}{\mathbf{I}}
\newcommand{\D}{\mathbf{D}}
\newcommand{\bS}{\mathbf{S}}
\newcommand{\N}{\mathcal{N}}
\renewcommand{\P}{\mathsf{P}}
\usepackage{amstext}

%%% links
\usepackage{hyperref}

%%% captions & tables
%% <FIXME>: conflics with figure definition in chapman.cls
%%\usepackage[format=hang,margin=10pt,labelfont=bf]{caption}
%% </FIMXE>
\usepackage{longtable}
\usepackage{rotating}

%%% Bibliography
\usepackage[round,comma]{natbib}
\renewcommand{\refname}{References \addcontentsline{toc}{chapter}{References}}
\citeindexfalse

%%% texi2dvi complains that \newblock is undefined, hm...
\def\newblock{\hskip .11em plus .33em minus .07em}

%%% Example sections
\newcounter{exercise}[chapter]
\setcounter{exercise}{0}
\newcommand{\exercise}{\item{\stepcounter{exercise} Ex.
                       \arabic{chapter}.\arabic{exercise} }}


%% URLs
\newcommand{\curl}[1]{\begin{center} \url{#1} \end{center}}

%%% for manual corrections
%\renewcommand{\baselinestretch}{2}

%%% plot sizes
\setkeys{Gin}{width=0.95\textwidth}

%%% color
\usepackage{color}

%%% hyphenations
\hyphenation{drop-out}

%%% new bidirectional quotes need
\usepackage[utf8]{inputenc}
\begin{document}
%% Title page
\title{A Handbook of Statistical Analyses Using \R}
\author{Brian S. Everitt and Torsten Hothorn}
\maketitle
%%\VignetteIndexEntry{Chapter Density Estimation}
%%\VignetteDepends{mclust,flexmix}
\setcounter{chapter}{6}
\SweaveOpts{prefix.string=figures/HSAUR,eps=FALSE}
<<setup, echo = FALSE, results = hide>>=
rm(list = ls())
if (!file.exists("tables")) dir.create("tables")
set.seed(290875)
options(prompt = "R> ", width = 63, # digits = 4,
    SweaveHooks = list(leftpar = function()
        par(mai = par("mai") * c(1, 1.05, 1, 1))))
HSAURpkg <- require("HSAUR")
if (!HSAURpkg) stop("cannot load package ", sQuote("HSAUR"))
rm(HSAURpkg)
@
\pagestyle{headings}
<<DE-setup, echo = FALSE, results = hide>>=
x <- library("KernSmooth")
### don't attach mclust since it masks stats:::density
### x <- library("mclust")
x <- library("flexmix")
x <- library("boot")
@
\chapter[Density Estimation]{Density Estimation: Erupting Geysers and Star
Clusters \label{DE}}
\section{Introduction}
\section{Density Estimation}
\section{Analysis Using \R{}}
\numberSinput
\begin{figure}
\begin{center}
<<DE-faithful-density, echo = TRUE, fig = TRUE, height = 4>>=
data("faithful", package = "datasets")
x <- faithful$waiting
layout(matrix(1:3, ncol = 3))
hist(x, xlab = "Waiting times (in min.)", ylab = "Frequency",
     probability = TRUE, main = "Gaussian kernel", border = "gray")
lines(density(x, width = 12), lwd = 2)
rug(x)
hist(x, xlab = "Waiting times (in min.)", ylab = "Frequency",
     probability = TRUE, main = "Rectangular kernel", border = "gray")
lines(density(x, width = 12, window = "rectangular"), lwd = 2)
rug(x)
hist(x, xlab = "Waiting times (in min.)", ylab = "Frequency",
     probability = TRUE, main = "Triangular kernel", border = "gray")
lines(density(x, width = 12, window = "triangular"), lwd = 2)
rug(x)
@
\caption{Density estimates of the geyser eruption data imposed on a histogram
of the data. \label{DE:faithfuldens}}
\end{center}
\end{figure}
\rawSinput
\begin{figure}
\begin{center}
<<DE-CYGOB1-contour, echo = TRUE, fig = TRUE>>=
library("KernSmooth")
data("CYGOB1", package = "HSAUR")
CYGOB1d <- bkde2D(CYGOB1, bandwidth = sapply(CYGOB1, dpik))
contour(x = CYGOB1d$x1, y = CYGOB1d$x2, z = CYGOB1d$fhat,
        xlab = "log surface temperature", ylab = "log light intensity")
@
\caption{A contour plot of the bivariate density estimate of the \Robject{CYGOB1} data,
         i.e., a two-dimensional graphical display for a three-dimensional
         \newline problem.
         \label{DE:CYGOB12Dcontour}}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
<<DE-CYGOB1-persp, echo = TRUE, fig = TRUE>>=
persp(x = CYGOB1d$x1, y = CYGOB1d$x2, z = CYGOB1d$fhat,
      xlab = "log surface temperature", ylab = "log light intensity",
      zlab = "estimated density", theta = -35, axes = TRUE, box = TRUE)
@
\caption{The bivariate density estimate of the \Robject{CYGOB1} data, here shown in a
         three-dimensional fashion using the \Rcmd{persp} function.
         \label{DE:CYGOB12Dpersp}}
\end{center}
\end{figure}
\subsection{A Parametric Density Estimate for the Old Faithful Data
\label{DE-waiting}}
<<DE-faithful-optim, echo = TRUE, results = hide>>=
logL <- function(param, x) {
    d1 <- dnorm(x, mean = param[2], sd = param[3])
    d2 <- dnorm(x, mean = param[4], sd = param[5])
    -sum(log(param[1] * d1 + (1 - param[1]) * d2))
}
startparam <- c(p = 0.5, mu1 = 50, sd1 = 3, mu2 = 80, sd2 = 3)
opp <- optim(startparam, logL, x = faithful$waiting, method = "L-BFGS-B",
             lower = c(0.01, rep(1, 4)),
             upper = c(0.99, rep(200, 4)))
opp
@
<<DE-faithful-optim-print, echo = FALSE>>=
print(opp[names(opp) != "message"])
@
Of course, optimising the appropriate likelihood `by hand' %'
is not very convenient. In fact, (at least) two packages offer high-level
functionality for estimating mixture models. The first one is package
\Rpackage{mclust} \citep{PKG:mclust} implementing the methodology described
in \cite{HSAUR:FraleyRaftery2002}. Here, a Bayesian information criterion
(BIC) is applied to choose the form of the mixture model:
\index{Bayesian Information Criterion (BIC)}
<<DE-attach-mclust, echo = FALSE, results = hide>>=
library("mclust")
@
<<DE-faithful-mclust, echo = TRUE>>=
library("mclust")
mc <- Mclust(faithful$waiting)
mc
@
and the estimated means are
<<DE-faithful-mclust-mu, echo = TRUE>>=
mc$mu
@
with estimated standard deviation (found to be equal within both groups)
<<DE-faithful-mclust-para, echo = TRUE>>=
sqrt(mc$sigmasq)
@
<<DE-detach-mclust, echo = FALSE, results = hide>>=
detach(package:mclust)
@
The proportion is $\hat{p} = \Sexpr{round(mc$pro[1], 2)}$. The second package is called
\Rpackage{flexmix} whose functionality is described by
\cite{HSAUR:Leisch2004}.
A mixture of two normals can be fitted using
<<DE-faithful-flexmix, echo = TRUE>>=
library("flexmix")
fl <- flexmix(waiting ~ 1, data = faithful, k = 2)
@
with $\hat{p} = \Sexpr{round(fl@prior, 2)}$ and estimated parameters
<<DE-faithful-flexmix-parameters, echo = TRUE>>=
parameters(fl, component = 1)
parameters(fl, component = 2)
@
\begin{figure}
\begin{center}
<<DE-faithful-2Dplot, echo = TRUE, fig = TRUE>>=
opar <- as.list(opp$par)
rx <- seq(from = 40, to = 110, by = 0.1)
d1 <- dnorm(rx, mean = opar$mu1, sd = opar$sd1)
d2 <- dnorm(rx, mean = opar$mu2, sd = opar$sd2)
f <- opar$p * d1 + (1 - opar$p) * d2
hist(x, probability = TRUE, xlab = "Waiting times (in min.)",
     border = "gray", xlim = range(rx), ylim = c(0, 0.06), main = "")
lines(rx, f, lwd = 2)
lines(rx, dnorm(rx, mean = mean(x), sd = sd(x)), lty = 2, lwd = 2)
legend(50, 0.06, legend = c("Fitted two-component mixture density",
                            "Fitted single normal density"),
       lty = 1:2, bty = "n")
@
\caption{Fitted normal density and two-component normal mixture for geyser
eruption data. \label{DE:2Dplot}}
\end{center}
\end{figure}
\index{Bootstrap approach|(}
We can get standard errors for the five parameter estimates
by using a bootstrap approach \citep[see][]{HSAUR:EfronTibshirani1993}.
The original data are slightly perturbed by drawing $n$ out of $n$
observations \stress{with replacement} and those artificial replications of
the original data are called \stress{bootstrap samples}. Now, we can fit the
mixture for each bootstrap sample and assess the variability of the
estimates, for example using confidence intervals.
\index{Confidence interval!derived from bootstrap samples}
Some suitable \R{} code based on the \Rcmd{Mclust}
function follows. First, we define a function that,
for a bootstrap sample \Robject{indx}, fits a two-component mixture model
and returns $\hat{p}$ and the estimated means (note that we need to make
sure that we always get an estimate of $p$, not $1 - p$):
<<DE-faithful-boot, echo = TRUE>>=
library("boot")
fit <- function(x, indx) {
    a <- Mclust(x[indx], minG = 2, maxG = 2)
    if (a$pro[1] < 0.5)
        return(c(p = a$pro[1], mu1 = a$mu[1], mu2 = a$mu[2]))
    return(c(p = 1 - a$pro[1], mu1 = a$mu[2], mu2 = a$mu[1]))
}
@
The function \Rcmd{fit} can now be fed into the \Rcmd{boot} function \citep{PKG:boot}
for bootstrapping (here $1000$ bootstrap samples are drawn)
\begin{Schunk}
\begin{Sinput}
R> bootpara <- boot(faithful$waiting, fit, R = 1000)
\end{Sinput}
\end{Schunk}
<<DE-faithful-bootrun, echo = FALSE>>=
bootparafile <- file.path(.find.package("HSAUR"), "cache", "DE-bootpara.rda")
if (file.exists(bootparafile)) {
    load(bootparafile)
} else {
    bootpara <- boot(faithful$waiting, fit, R = 1000)
}
@
We assess the variability of our estimates $\hat{p}$ by
means of adjusted bootstrap percentile (BCa) confidence intervals, which for
$\hat{p}$ can be obtained from
<<DE-faithful-p-ci, echo = TRUE>>=
boot.ci(bootpara, type = "bca", index = 1)
@
We see that there is a reasonable variability in the mixture model, however,
the means in the two components are rather stable, as can be seen from
<<DE-faithful-mu1-ci, echo = TRUE>>=
boot.ci(bootpara, type = "bca", index = 2)
@
for $\hat{\mu}_1$ and for $\hat{\mu}_2$ from
<<DE-faithful-mu2-ci, echo = TRUE>>=
boot.ci(bootpara, type = "bca", index = 3)
@
Finally, we show a graphical representation of both the bootstrap
distribution of the mean estimates \stress{and} the corresponding confidence
intervals. For convenience, we define a function for plotting, namely
<<DE-bootplot, echo = TRUE>>=
bootplot <- function(b, index, main = "") {
    dens <- density(b$t[,index])
    ci <- boot.ci(b, type = "bca", index = index)$bca[4:5]
    est <- b$t0[index]
    plot(dens, main = main)
    y <- max(dens$y) / 10
    segments(ci[1], y, ci[2], y, lty = 2)
    points(ci[1], y, pch = "(")
    points(ci[2], y, pch = ")")
    points(est, y, pch = 19)
}
@
The element \Robject{t} of an object created by \Rcmd{boot} contains the
bootstrap replications of our estimates, i.e., the values computed by
\Rcmd{fit} for each of the $1000$ bootstrap samples of the geyser data.
First, we plot a simple density estimate and then construct a line
representing the confidence interval.
We apply this function to the bootstrap distributions of our estimates
$\hat{\mu}_1$ and $\hat{\mu}_2$ in Figure~\ref{DE-bootplot}.
\begin{figure}
\begin{center}
<<DE-faithful-boot-plot, echo = TRUE, fig = TRUE, height = 4>>=
layout(matrix(1:2, ncol = 2))
bootplot(bootpara, 2, main = expression(mu[1]))
bootplot(bootpara, 3, main = expression(mu[2]))
@
\caption{Bootstrap distribution and confidence intervals for the mean estimates
         of a two-component mixture for the geyser data. \label{DE-bootplot}}
\end{center}
\end{figure}
\index{Bootstrap approach|)}
\bibliographystyle{../LaTeX/refstyle}
\bibliography{../LaTeX/HSAUR}
\end{document}
