%\VignetteIndexEntry{Mixture Analysis of the Old Faithful Geyser Data}

\documentclass[a4paper, 11pt]{article}

\newcommand{\FIGKEEPDIR}{./figuresKeep}
\newcommand{\FIGDIR}{./figures}

\usepackage{amssymb, amsmath}
\usepackage{color}
\usepackage{graphicx} 
\usepackage{psfrag}
\usepackage{natbib}
\usepackage{SweaveAK}
%\usepackage{Sweave}
\usepackage[breaklinks=true]{hyperref}

\setlength{\textwidth}{6.5in} 
\setlength{\textheight}{9in}
\setlength{\hoffset}{-0.6in} 
\setlength{\parskip}{1ex}
\setlength{\parindent}{0pt}

\newcommand{\AKpack}{\textsf{mixAK}}
\newcommand{\Ritem}{\textcolor{blue}{\textsf{R}$\boldsymbol{\Rightarrow}$}\ }

\title{Mixture Analysis of the \textsf{Old Faithful Geyser} Data Using the Package {\AKpack}}
\author{Arno\v{s}t Kom\'arek}
%\address{Dept. of Probability and Mathematical Statistics, 
%         Faculty of Mathematics and Physics,
%         Charles University in Prague, 
%         Sokolovsk\'a 83, CZ--186$\;$75, Praha 8, Czech Republic}

\begin{document}

\maketitle
\par\rule{\textwidth}{1pt}

\bigskip
This document supplements a~paper \citet{KomarekNMix} and shows an~analysis of the Old Faithful 
Geyser data introduced in \citet{Hardle91} using the {\Rko} package \Rfun{\AKpack}. The data have
been analysed using mixtures by several researchers, e.g., \citet{Stephens00}, 
\citet{DellaportasPapageorgiou06}. 

\section{Introduction}
\textbullet\ 
Due to the fact that some code (especially MCMC) is time consuming, some code chunks found in this
vignette are not run when compiling the package. You should set the variable \Rvar{RUN.TIMECONSUMING.CODE}
to \Rvar{TRUE} to run full MCMC and related code.

\textbullet\ Having run full MCMC and related code, setting the variable \Rvar{RUN.ALLOUT}
to \Rvar{TRUE} will cause that all output shown in this vignette is re-created and not taken
from previously computed results.

\textbullet\ Results based on full MCMC runs are stored (chains excluded) in \Rfile{Faithful-Result.RData}
in \Rfile{./RESULT\_{}OBJ} directory and are used to create majority of this vignette
at the time of compilation of the package.

\bigskip
\Ritem Setting variables \Rvar{RUN.ALLOUT} and \Rvar{RUN.TIMECONSUMING.CODE}.
<<What should be created in this Sweave run>>=
RUN.TIMECONSUMING.CODE <- FALSE
RUN.ALLOUT <- FALSE
@ 

\Ritem Directory to store postscript files with figures.
Figures which require chains are stored in \Rvar{FIGKEEPDIR} directory all other figures
are stored in \Rvar{FIGDIR} directory.
<<Directory to store figures>>=
FIGDIR <- "./figures/"
FIGKEEPDIR <- "./figuresKeep/"
@ 

\Ritem Directories with results computed in past.
Objects with chains will be stored in directory specified by variable \Rvar{RESULTDIR}.
All other objects will be stored in directory \Rvar{RESULT2DIR}.
%%% Scripts to compute main objects:  /home/komarek/Rlib/mixAK/Vignette/RMCMC/Faithful-MCMC.R
<<Directory with results computed in past, keep.source=TRUE>>=
RESULTDIR <- "/home/komarek/RESULT_OBJ/mixAK-Faithful-S081115/"
RESULT2DIR <- "./RESULT_OBJ/"      ### available in package as /inst/doc/RESULT_OBJ
@

\Ritem Display options.
<<Options>>=
options(width=80)
@ 

\Ritem Load results computed in past. Variable \Rvar{Kshow} determines which value for fixed
number of components $K$ is used in section~\ref{sec:ModelFixedK}.
<<Load results computed in past, keep.source=TRUE>>=
Kshow <- 3
if (RUN.ALLOUT){
  load(paste(RESULT2DIR, "Faithful-Result.RData", sep=""))            
     ## contains ModelK (without chains), MPDensModelK, JPDensModelK
  load(paste(RESULTDIR, "Faithful-Model0", Kshow, ".RData", sep=""))  
     ## contains Model0=ModelK[[Kshow]] (chains included)
  MPDensModel0 <- MPDensModelK[[Kshow]]
  JPDensModel0 <- JPDensModelK[[Kshow]]  
}else{
  load(paste(RESULT2DIR, "Faithful-Result.RData", sep=""))  
     ## contains ModelK (without chains), MPDensModelK, JPDensModelK
  Model0 <- ModelK[[Kshow]]
  MPDensModel0 <- MPDensModelK[[Kshow]]
  JPDensModel0 <- JPDensModelK[[Kshow]]    
}  
@ 

\Ritem Load the package {\AKpack} and packages \Rfun{coda} and \Rfun{colorspace}. 
Package \Rfun{coda} is used to perform some basic convergence diagnostics, package \Rfun{colorspace}
is used to draw nicer image plots with estimated bivariate densities.
<<Load needed packages>>=
library("mixAK")
library("coda")
library("colorspace")
@ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Exploration of the data}

\Ritem The data are read and summarized as follows. 
<<Read the data and compute a brief summary>>=
data("Faithful", package="mixAK")
summary(Faithful)
@

\Ritem Additionally, Figure~\ref{fig01} shows the scatterplot and histograms of the data.
<<Draw scatterplot and histogram, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful01.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(bty="n")
layout(matrix(c(0,1,1,1,1,0, 2,2,2,3,3,3), nrow=2, byrow=TRUE))
plot(Faithful, col="red", pch=16, 
     xlab="Eruptions (min)", ylab="Waiting (min)")
hist(Faithful$eruptions, prob=TRUE, col="sandybrown", 
     xlab="Eruptions (min)", ylab="Density", main="", 
     breaks=seq(1.4, 5.6, by=0.3))
hist(Faithful$waiting, prob=TRUE, col="sandybrown", 
     xlab="Waiting (min)", ylab="Density", main="")
dev.off()
@ 
\begin{figure}[p!]\centering
\includegraphics[width=\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful01.ps}
\caption{Scatterplot and histograms of Faithful data.}\label{fig01}
\end{figure}


%%%%% %%%%%%%%%%%%%%%
\clearpage
\section{Preparation of the MCMC}
\Ritem Length of the MCMC simulation for all models in this document 
(burn-in of 100\,000 iterations, additional 500\,000 iterations
are kept for the inference, thinning of 1:10):
<<Length of the MCMC, keep.source=TRUE>>=
nMCMC <- c(burn=100000, keep=500000, thin=10, info=10000)
@ 
\Ritem Grid of values where we evaluate and subsequently plot the predictive density 
for all models in this document:
<<Grid of values for the predictive density, keep.source=TRUE>>=
ygrid <- list(eruptions=seq(1, 6, length=100), 
              waiting=seq(40, 100, length=100))
@ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Model with three mixture components}\label{sec:ModelFixedK}
In this section, we will fit a~mixture model with $K=3$ components.

<<TESTING: Prior distribution and model with three components and posterior predictive density, echo=FALSE, results=hide>>=
set.seed(770328)
Prior0 <- list(priorK="fixed", Kmax=3)
TEST.Model0 <- NMixMCMC(y0=Faithful, prior=Prior0, nMCMC=c(burn=100, keep=100, thin=2, info=50))
TEST.MPDensModel0 <- NMixPredDensMarg(TEST.Model0[[1]], grid=ygrid)
TEST.JPDensModel0 <- NMixPredDensJoint2(TEST.Model0[[1]], grid=ygrid)
@ 

%%%%% %%%%%%%%%%%%%%%
\subsection{Specification of the prior distributions and MCMC simulation}
\Ritem The minimal specification of the prior distribution and running MCMC with default values 
for all prior parameters and $K=3$:
<<Prior distribution>>=
Prior0 <- list(priorK="fixed", Kmax=3)
@ 

\Ritem Two chains will be generated since the argument \Rarg{PED} is set to \Rvar{TRUE}
(output is shown from MCMC simulation performed by author):
<<Model with three components, keep.source=TRUE, results=hide>>=
if (RUN.TIMECONSUMING.CODE){
  set.seed(777988621)  
  Model0 <- NMixMCMC(y0=Faithful, prior=Prior0, nMCMC=nMCMC, PED=TRUE)
}  
@ 
%%% CP time: Snehurka
\begin{Schunk}
\begin{Soutput}
Chain number 1
==============
MCMC sampling started on Sat Nov 15 21:18:35 2008.
Burn-in iteration 100000
Iteration 600000
MCMC sampling finished on Sat Nov 15 21:38:44 2008.

Chain number 2
==============
MCMC sampling started on Sat Nov 15 21:38:52 2008.
Burn-in iteration 100000
Iteration 600000
MCMC sampling finished on Sat Nov 15 21:58:44 2008.

Computation of penalized expected deviance started on Sat Nov 15 21:58:53 2008.
Computation of penalized expected deviance finished on Sat Nov 15 22:09:28 2008.
\end{Soutput}
\end{Schunk}
\Ritem The prior distribution for the function \Rfun{NMixMCMC} was the same as with
<<The prior distribution was the same as with, keep.source=TRUE, eval=FALSE>>=
Prior0 <- list(priorK="fixed", Kmax=3, 
               delta=1, 
               priormuQ="independentC",
               xi=c(-0.1207, -0.1028), D=diag(c(9.4033, 15.1983)),
               zeta=3, g=0.2, h=c(1.0635, 0.6580))
@ 
Note that due to the fact that the argument \Rarg{scale} has not been specified in function 
\Rfun{NMixMCMC}, the MCMC has been run on standardized data. Adding 
\Rarg{scale=list(shift=0, scale=1)} would lead to the MCMC run on the original data.


%%%%% %%%%%%%%%%%%%%%
\clearpage
\subsection{Posterior inference}

% --------------------------------------
\Ritem Basic posterior summary of the fitted model is obtained using the command \Rfun{print(Model0)}.

\Ritem Quantities shown in section labeled 
``\textsf{\color{blue} Penalized expected deviance}'' are computed from two sampled chains and
have the following meaning: \Rvar{D.expect} is $\hat{D}_e$ from \citet{KomarekNMix},
\Rvar{p(opt)} is estimated optimism $\hat{p}_{opt}$ computed with unit weights,
\Rvar{PED} equals \Rvar{D.expect} $+$ \Rvar{p(opt)} gives the estimate of penalized
expected deviance with optimism computed without the use of importance sampling.
Further, \Rvar{wp(opt)} is estimated optimism $\hat{p}_{opt}$ computed as described
in \citet{KomarekNMix}, i.e., using importance sampling. Finally, \Rvar{wPED} equals
\Rvar{D.expect} $+$ \Rvar{wp(opt)} gives the estimate of penalized expected deviance
as described in \citet{KomarekNMix}.

\Ritem Quantities shown in section labeled
``\textsf{\color{blue} Deviance information criteria}'' are computed separately from the first
and the second sampled chain. They have the following meaning:
\Rvar{DIC} is deviance information criterion denoted as DIC in \citet{KomarekNMix},
\Rvar{pD} is the effective dimension $p_D$, \Rvar{D.bar} is approximated posterior mean
$\overline{D}$ of the deviance and \Rvar{D.in.bar} is $\tilde{D}$ -- deviance evaluated
in the ``estimate''.

\Ritem Section ``\textsf{\color{blue} Posterior summary statistics for moments of mixture
for original data}'' gives posterior summary statistics for 
$\mbox{E}(\boldsymbol{Y}) = \boldsymbol{m} + \boldsymbol{S}\mbox{E}(\boldsymbol{Y}^*)$
and quantities derived from $\mbox{var}(\boldsymbol{Y}) = \boldsymbol{S}\,\mbox{var}(\boldsymbol{Y}^*)\,\boldsymbol{S}'$
in the notation of \citet{KomarekNMix}, separately for each generated chain.

<<Model0: Basic posterior summary>>=
print(Model0)
@

% --------------------------------------
\clearpage
\Ritem Computation of the marginal (univariate) predictive densities
(separately for chain~1 and chain~2):
<<Model0: Marginal predictive densities, keep.source=TRUE>>=
if (RUN.TIMECONSUMING.CODE){
  MPDensModel0 <- list()
  MPDensModel0[[1]] <- NMixPredDensMarg(Model0[[1]], grid=ygrid)  
  MPDensModel0[[2]] <- NMixPredDensMarg(Model0[[2]], grid=ygrid)  
}  
@

% --------------------------------------
\Ritem Default \Rfun{plot} method for the computed object (see Figure~\ref{fig02}):
<<Model0: Plot of the marginal predictive density, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful02.ps", sep=""), width=7, height=5, 
           horizontal=FALSE)
plot(MPDensModel0[[1]])
dev.off()
@ 

% --------------------------------------
\Ritem Computation of the joint (bivariate) predictive densities
(separately for chain~1 and chain~2):
<<Model0: Joint predictive density, keep.source=TRUE>>=
if (RUN.TIMECONSUMING.CODE){
  JPDensModel0 <- list()
  JPDensModel0[[1]] <- NMixPredDensJoint2(Model0[[1]], grid=ygrid)  
  JPDensModel0[[2]] <- NMixPredDensJoint2(Model0[[2]], grid=ygrid)  
}  
@


% --------------------------------------
\Ritem Default \Rfun{plot} method for the computed object (see Figure~\ref{fig03}):
<<Model0: Plot of the joint predictive density, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful03.ps", sep=""), width=7, height=5, 
           horizontal=FALSE)
plot(JPDensModel0[[1]])
dev.off()
@ 

\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.45\textheight]{\FIGDIR/figFaithful02.ps}
\caption{Default \Rfun{plot} method for the marginal predictive densities based on the model 
with three mixture components (results from chain~1).}\label{fig02}
\includegraphics[width=1.0\textwidth, height=0.45\textheight]{\FIGDIR/figFaithful03.ps}
\caption{Default \Rfun{plot} method for the joint predictive density based on the model 
with three mixture components (results from chain~1).}\label{fig03}
\end{figure}


%%%%% %%%%%%%%%%%%%%%
\clearpage
\subsection{Nicer figures of posterior predictive densities}

% --------------------------------------
\Ritem Plot of the marginal predictive densities together with the histograms of the data,
separate figures for two generated chains (see Figure~\ref{fig04}):
<<Model0: Plot of the marginal predictive density with the histogram, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful04.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(mfrow=c(2, 2), bty="n")
ylimE <- c(0, max(c(MPDensModel0[[1]]$dens[[1]], MPDensModel0[[2]]$dens[[1]])))
ylimW <- c(0, max(c(MPDensModel0[[1]]$dens[[2]], MPDensModel0[[2]]$dens[[2]])))
for (CH in 1:2){
  hist(Faithful$eruptions, prob=TRUE, col="sandybrown", 
       xlab="Eruptions (min)", ylab="Density", main=paste("Chain", CH), 
       breaks=seq(1.4, 5.6, by=0.3), ylim=ylimE)
  lines(MPDensModel0[[CH]]$x$x1, MPDensModel0[[CH]]$dens[[1]], col="darkblue", lwd=2)
  hist(Faithful$waiting, prob=TRUE, col="sandybrown", 
       xlab="Waiting (min)", ylab="Density", main="", ylim=ylimW)
  lines(MPDensModel0[[CH]]$x$x2, MPDensModel0[[CH]]$dens[[2]], col="darkblue", lwd=2)
}  
dev.off()
@ 

% --------------------------------------
\Ritem Contour plot of the joint predictive density together with the scatterplot of the data,
separate figures for two generated chains (see Figure~\ref{fig05}):
<<Model0: Contour plot of the joint predictive density with the scatterplot, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful05.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(mfrow=c(2, 1), bty="n")
for (CH in 1:2){
  plot(Faithful, col="red", xlab="Eruptions (min)", ylab="Waiting (min)", 
       main=paste("Chain", CH))
  contour(JPDensModel0[[CH]]$x$x1, JPDensModel0[[CH]]$x$x2, 
          JPDensModel0[[CH]]$dens[["1-2"]], 
          col="darkblue", add=TRUE)
}  
dev.off()
@
%$

\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful04.ps}
\caption{Marginal predictive densities (blue line) based on the model with three mixture 
components, separately for chain~1 and chain~2.}\label{fig04}
\end{figure}

\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful05.ps}
\caption{Contour plot of the joint predictive density based on the model with three mixture 
components, separately for chain~1 and chain~2.}\label{fig05}
\end{figure}

% --------------------------------------
\clearpage
\Ritem Image plot of the joint predictive density together with the scatterplot of the data,
separate figures for two generated chains (see Figure~\ref{fig06}):
<<Model0: Image plot of the joint predictive density with the scatterplot, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful06.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(mfrow=c(2, 1), bty="n")
for (CH in 1:2){
  plot(Faithful, col="darkblue", xlab="Eruptions (min)", ylab="Waiting (min)", 
       main=paste("Chain", CH))
  image(JPDensModel0[[CH]]$x$x1, JPDensModel0[[CH]]$x$x2, 
        JPDensModel0[[CH]]$dens[["1-2"]], add=TRUE, 
        col=rev(heat_hcl(33, c=c(80, 30), l=c(30, 90), power=c(1/5, 1.3))))  
  points(Faithful, col="darkblue")
}  
dev.off()
@ 
%$
Note that package \Rfun{colorspace} is needed to specify the colors in the plot.

\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful06.ps}
\caption{Image plot of the joint predictive density based on the model with three mixture 
components, separately for chain~1 and chain~2.}\label{fig06}
\end{figure}


%%%%% %%%%%%%%%%%%%%%
\clearpage
\subsection{Convergence diagnostics}
\Ritem Single chain convergence diagnostics using chain~1 will be shown here.
<<RJMCMC: Choose chain>>=
CH <- 1
@ 

\Ritem Converting the chains into \Rfun{mcmc} objects to be used in the package \Rfun{coda}:
<<Model0: Create mcmc objects, keep.source=TRUE>>=
if (RUN.ALLOUT){
  start <- Model0[[CH]]$nMCMC["burn"] + 1
  end <- Model0[[CH]]$nMCMC["burn"] + Model0[[CH]]$nMCMC["keep"]
  chgammaInv <- mcmc(Model0[[CH]]$gammaInv, start=start, end=end)
  chmixture <- mcmc(Model0[[CH]]$mixture, start=start, end=end)
  chdeviance <- mcmc(Model0[[CH]]$deviance, start=start, end=end)
}  
@ 

% --------------------------------------
\Ritem Traceplots for selected parameters (only the first set of traceplots is shown, see Figure~\ref{fig07}).

\Ritem Traceplots will be drawn for last $5\,000$ iterations only:
<<Choose iters to draw traceplots, keep.source=TRUE>>=
if (RUN.ALLOUT){
  tstart <- 495001
  tend <- 500000
  titers <- tstart:tend
  chgammaInv2 <- mcmc(Model0[[CH]]$gammaInv[titers,], start=tstart, end=tend)
  chmixture2 <- mcmc(Model0[[CH]]$mixture[titers,], start=tstart, end=tend)
  chdeviance2 <- mcmc(Model0[[CH]]$deviance[titers,], start=tstart, end=tend)
}
@ 
%$

<<Model0: Traceplots, keep.source=TRUE, results=hide>>=
if (RUN.ALLOUT){
  lwd <- 0.5
  postscript(paste(FIGKEEPDIR, "figFaithful07.ps", sep=""), width=7, height=10, 
             horizontal=FALSE)
  par(mfrow=c(2, 3), bty="n")
  traceplot(chmixture2[, "y.Mean.1"], smooth=FALSE, col="darkblue", lwd=lwd, 
            main="E(eruptions)")
  traceplot(chmixture2[, "y.Mean.2"], smooth=FALSE, col="darkblue", lwd=lwd, 
            main="E(waitings)")
  traceplot(chgammaInv2[, "gammaInv1"], smooth=FALSE, col="brown", lwd=lwd, 
            main="gamma^{-1}")
  traceplot(chmixture2[, "y.SD.1"], smooth=FALSE, col="darkgreen", lwd=lwd, 
            main="sd(eruptions)")
  traceplot(chmixture2[, "y.SD.2"], smooth=FALSE, col="darkgreen", lwd=lwd, 
            main="sd(waitings)")
  traceplot(chmixture2[, "y.Corr.2.1"], smooth=FALSE, col="red", lwd=lwd, 
            main="corr(eruptions, waitings)")
  dev.off()
}
@ 

\newpage
<<Model0: Traceplots, keep.source=TRUE, results=hide>>=
if (RUN.ALLOUT){
  postscript(paste(FIGKEEPDIR, "figFaithful08.ps", sep=""), width=7, height=10, 
             horizontal=FALSE)
  par(mfrow=c(2, 2), bty="n")
  traceplot(chdeviance2[, "LogL0"], smooth=FALSE, col="red", lwd=lwd, 
            main="Log(L0)")
  traceplot(chdeviance2[, "LogL1"], smooth=FALSE, col="red", lwd=lwd, 
            main="Log(L1)")
  traceplot(chdeviance2[, "dev.complete"], smooth=FALSE, col="red", lwd=lwd, 
            main="D(complete)")
  traceplot(chdeviance2[, "dev.observed"], smooth=FALSE, col="red", lwd=lwd, 
            main="D(observed)")
  dev.off()
}  
@ 

% --------------------------------------
\Ritem Posterior density estimates for selected parameters (see Figure~\ref{fig09}):
<<Model0: Density plots, keep.source=TRUE, results=hide>>=
if (RUN.ALLOUT){
  postscript(paste(FIGKEEPDIR, "figFaithful09.ps", sep=""), width=7, height=10, 
             horizontal=FALSE)
  par(mfrow=c(2, 3), bty="n")
  densplot(chmixture[, "y.Mean.1"], show.obs=FALSE, col="darkblue", 
           main="E(eruptions)")
  densplot(chmixture[, "y.Mean.2"], show.obs=FALSE, col="darkblue", 
           main="E(waitings)")
  densplot(chgammaInv[, "gammaInv1"], show.obs=FALSE, col="brown", 
           main="gamma^{-1}")
  densplot(chmixture[, "y.SD.1"], show.obs=FALSE, col="darkgreen", 
           main="sd(eruptions)")
  densplot(chmixture[, "y.SD.2"], show.obs=FALSE, col="darkgreen", 
           main="sd(waitings)")
  densplot(chmixture[, "y.Corr.2.1"], show.obs=FALSE, col="red", 
           main="corr(eruptions, waitings)")
  dev.off()
} 
@

\begin{figure}[p!]\centering
\rotatebox{270}{
  \includegraphics[width=0.9\textheight, height=\textwidth]{\FIGKEEPDIR/figFaithful07.ps}
}
\caption{Model with three mixture components. Traceplots for selected parameters.}\label{fig07}
\end{figure}

\begin{figure}[p!]\centering
\rotatebox{270}{
  \includegraphics[width=0.9\textheight, height=\textwidth]{\FIGKEEPDIR/figFaithful09.ps}
}
\caption{Model with three mixture components. Posterior density estimates for selected parameters.}\label{fig09}
\end{figure}


% --------------------------------------
\clearpage
\Ritem Autocorrelation plots for selected parameters (see Figure~\ref{fig10}):
<<Model0: Autocorrelation plots, keep.source=TRUE, results=hide>>=
if (RUN.ALLOUT){
  postscript(paste(FIGKEEPDIR, "figFaithful10.ps", sep=""), width=7, height=10, 
             horizontal=FALSE)
  par(mfrow=c(2, 3), bty="n")
  autocorr.plot(chmixture[, "y.Mean.1"], auto.layout=FALSE, 
                ask=FALSE, col="darkblue", main="E(eruptions)")
  autocorr.plot(chmixture[, "y.Mean.2"], auto.layout=FALSE, 
                ask=FALSE, col="darkblue", main="E(waitings)")
  autocorr.plot(chgammaInv[, "gammaInv1"], auto.layout=FALSE, 
                ask=FALSE, col="brown", main="gamma^{-1}")
  autocorr.plot(chmixture[, "y.SD.1"], auto.layout=FALSE, 
                ask=FALSE, col="darkgreen", main="sd(eruptions)")
  autocorr.plot(chmixture[, "y.SD.2"], auto.layout=FALSE, 
                ask=FALSE, col="darkgreen", main="sd(waitings)")
  autocorr.plot(chmixture[, "y.Corr.2.1"], auto.layout=FALSE, 
                ask=FALSE, col="red", main="corr(eruptions, waitings)")
  dev.off()
}
@

\begin{figure}[p!]\centering
\rotatebox{270}{
  \includegraphics[width=0.9\textheight, height=\textwidth]{\FIGKEEPDIR/figFaithful10.ps}
}
\caption{Model with three mixture components. Autocorrelation plots for selected parameters.}\label{fig10}
\end{figure}


% --------------------------------------
\clearpage
\Ritem In advance, we know that the posterior distribution has $3!$ symmetric modes which should all be visited
by the Markov chain. Indication of failing to visit all the modes is provided by exploration of posterior distribution
of mixture weights, means and variance-covariance matrices. If no identifiability constraints are imposed (as is done
in our implementation of MCMC), posterior distributions of mixture weights, means and variance-covariance matrices
should be identical for all mixture components if the chain suceeded to visit all the modes of the posterior.
The following figures show possibility for exploration of these posterior distributions
(see Figures \ref{fig16}, \ref{fig17}, \ref{fig18} for results). To make the size of the resulting pdf file
reasonable, all scatterplots are drawn for $5\,000$ randomly selected iterations only.
<<Model0: Scatterplots and histograms of mixture elements, keep.source=TRUE, results=hide>>=
if (RUN.ALLOUT){
  PCH <- 1;  CEX <- 0.5
  set.seed(770328)
  SELECT <- sample(end-start+1, size=5000, replace=FALSE)

  MuDens1 <- MuDens2 <- list()
  for (j in 1:3){
    MuDens1[[j]] <- density(Model0[[CH]]$mu[,(j-1)*2+1])
    MuDens2[[j]] <- density(Model0[[CH]]$mu[,j*2])
  }  
  LTY <- c(1, 2, 4)  
  COL <- c("blue", "red", "darkgreen")
  
  XLIM <- c(-2, 2);  YLIM <- c(-2, 2);  
  XLAB <- "Mean, margin 1";  YLAB <- "Mean, margin 2"  
  postscript(paste(FIGKEEPDIR, "figFaithful16.ps", sep=""), width=7, height=10, 
             horizontal=FALSE)
  par(bty="n")
  layout(matrix(c(1,0, 1,4, 2,4, 2,5, 3,5, 3,6), ncol=2, byrow=TRUE))  
  plot(Model0[[CH]]$mu[SELECT, "mu.1.1"], Model0[[CH]]$mu[SELECT, "mu.1.2"], 
       col="blue", pch=PCH, cex=CEX, xlim=XLIM, ylim=YLIM, 
       xlab=XLAB, ylab=YLAB, main=expression(mu[1]))
  plot(Model0[[CH]]$mu[SELECT, "mu.2.1"], Model0[[CH]]$mu[SELECT, "mu.2.2"], 
       col="blue", pch=PCH, cex=CEX, xlim=XLIM, ylim=YLIM, 
       xlab=XLAB, ylab=YLAB, main=expression(mu[2]))
  plot(Model0[[CH]]$mu[SELECT, "mu.3.1"], Model0[[CH]]$mu[SELECT, "mu.3.2"], 
       col="blue", pch=PCH, cex=CEX, xlim=XLIM, ylim=YLIM, 
       xlab=XLAB, ylab=YLAB, main=expression(mu[3]))
  #
  plot(MuDens1[[1]]$x, MuDens1[[1]]$y, type="l", lty=LTY[1], col=COL[1], 
       xlab="Mean, margin 1", ylab="Posterior density", xlim=XLIM, 
       main="Margin 1")
  for (j in 2:3) lines(MuDens1[[j]]$x, MuDens1[[j]]$y, lty=LTY[j], col=COL[j])  
  #
  plot(MuDens2[[1]]$x, MuDens2[[1]]$y, type="l", lty=LTY[1], col=COL[1], 
       xlab="Mean, margin 2", ylab="Posterior density", xlim=YLIM, 
       main="Margin 2")
  for (j in 2:3) lines(MuDens2[[j]]$x, MuDens2[[j]]$y, lty=LTY[j], col=COL[j])
  #
  plot(0:100, 0:100, type="n", xlab="", ylab="", xaxt="n", yaxt="n")
  legend(10, 99, paste("Comp.", 1:3), lty=LTY, col=COL, bty="n", y.intersp=1)    
  dev.off()
  
  XLIM <- c(0, 1)
  YLIM <- c(0, 1.5)
  XLAB <- "Variance, margin 1"
  YLAB <- "Variance, margin 2"   
  postscript(paste(FIGKEEPDIR, "figFaithful17.ps", sep=""), width=7, height=10, 
             horizontal=FALSE)
  par(bty="n")
  layout(matrix(c(0,1,1,0, 2,2,3,3), nrow=2, byrow=TRUE))
  plot(Model0[[CH]]$Sigma[SELECT, "Sigma1.1.1"], 
       Model0[[CH]]$Sigma[SELECT, "Sigma1.2.2"], 
       col="red", pch=PCH, cex=CEX, xlim=XLIM, ylim=YLIM, 
       xlab=XLAB, ylab=YLAB, main=expression(Sigma[1]))
  plot(Model0[[CH]]$Sigma[SELECT, "Sigma2.1.1"], 
       Model0[[CH]]$Sigma[SELECT, "Sigma2.2.2"], 
       col="red", pch=PCH, cex=CEX, xlim=XLIM, ylim=YLIM, 
       xlab=XLAB, ylab=YLAB, main=expression(Sigma[2]))
  plot(Model0[[CH]]$Sigma[SELECT, "Sigma3.1.1"], 
       Model0[[CH]]$Sigma[SELECT, "Sigma3.2.2"], 
       col="red", pch=PCH, cex=CEX, xlim=XLIM, ylim=YLIM, 
       xlab=XLAB, ylab=YLAB, main=expression(Sigma[3]))
  dev.off()
  
  XLIM <- c(0, 1)
  YLIM <- c(0, 4.5)
  XLAB <- "Weight"
  YLAB <- "Density"
  postscript(paste(FIGKEEPDIR, "figFaithful18.ps", sep=""), width=7, height=10, 
             horizontal=FALSE)
  par(bty="n")
  layout(matrix(c(0,1,1,0, 2,2,3,3), nrow=2, byrow=TRUE))
  hist(Model0[[CH]]$w[, "w1"], prob=TRUE, col="sandybrown", 
       xlim=XLIM, ylim=YLIM, xlab=XLAB, ylab=YLAB, main=expression(w[1]))  
  hist(Model0[[CH]]$w[, "w2"], prob=TRUE, col="sandybrown", 
       xlim=XLIM, ylim=YLIM, xlab=XLAB, ylab=YLAB, main=expression(w[2]))  
  hist(Model0[[CH]]$w[, "w3"], prob=TRUE, col="sandybrown", 
       xlim=XLIM, ylim=YLIM, xlab=XLAB, ylab=YLAB, main=expression(w[3]))  
  dev.off()
}  
@ 
%$

\begin{figure}[p!]\centering
\includegraphics[width=\textwidth, height=0.9\textheight]{\FIGKEEPDIR/figFaithful16.ps}
\caption{Model with three mixture components. Scatterplots of sampled mixture means ($5\,000$ randomly selected
  iterations) without imposing any identifiability constraints and posterior densities
  for marginal mixture means in the three components.}\label{fig16}
\end{figure}

\begin{figure}[p!]\centering
\includegraphics[width=\textwidth, height=0.9\textheight]{\FIGKEEPDIR/figFaithful17.ps}
\caption{Model with three mixture components. Scatterplots of sampled mixture variances ($5\,000$ randomly selected
  iterations) without imposing any identifiability constraints.}\label{fig17}
\end{figure}

\begin{figure}[p!]\centering
\includegraphics[width=\textwidth, height=0.9\textheight]{\FIGKEEPDIR/figFaithful18.ps}
\caption{Model with three mixture components. Histograms of sampled mixture weights 
  without imposing any identifiability constraints.}\label{fig18}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Models with different fixed numbers of components}
In this section, we will fit a~mixture model for $K=1,\dots, 10$, compare the deviance based 
quantities and predictive densities.

% --------------------------------------
\Ritem Running the MCMC simulation for $K = 1,\dots, 10$ 
(output printed during the MCMC run on the screen not shown),
computation of predictive densities. \par
\Ritem After predictive densities are computed, we remove all chains from resulting objects 
(to save some memory).
<<Fixed K MCMC, keep.source=TRUE, results=hide>>=
if (RUN.TIMECONSUMING.CODE){
  Seed <- c(777988621, 777988621, 777988621, 777988621, 780830,
            780830, 777988621, 777988621, 777988621, 777988621)
#
  Keep <- c("iter", "nMCMC", "dim", "prior", "init", "RJMCMC", 
            "scale", "state", "freqK", "propK", "DIC", "moves", 
            "pm.y", "pm.z", "pm.indDev", "pred.dens", "summ.y.Mean", 
            "summ.y.SDCorr", "summ.z.Mean", "summ.z.SDCorr")
#
  ModelK <- list()
  MPDensModelK <- list()
  JPDensModelK <- list()
  for (k in 1:10){
    set.seed(Seed[k])    
    cat(paste("K = ", k, "\n-------------------------------\n", sep=""))  
    PriorNow <- Prior0
    PriorNow$Kmax <- k
    ModelK[[k]] <- NMixMCMC(y0=Faithful, prior=PriorNow, nMCMC=nMCMC, PED=TRUE)
#
    cat(paste("\nComputation of marginal pred. densities started on ", date(), 
              "\n", sep=""))  
    MPDensModelK[[k]] <- list()
    MPDensModelK[[k]][[1]] <- NMixPredDensMarg(ModelK[[k]][[1]], grid=ygrid)
    MPDensModelK[[k]][[2]] <- NMixPredDensMarg(ModelK[[k]][[2]], grid=ygrid)
    cat(paste("Computation of marginal pred. densities finished on ", date(), 
              "\n\n", sep=""))  
#
    cat(paste("Computation of joint pred. densities started on ", date(), 
              "\n", sep=""))    
    JPDensModelK[[k]] <- list()
    JPDensModelK[[k]][[1]] <- NMixPredDensJoint2(ModelK[[k]][[1]], grid=ygrid)
    JPDensModelK[[k]][[2]] <- NMixPredDensJoint2(ModelK[[k]][[2]], grid=ygrid)
    cat(paste("Computation of joint pred. densities finished on ", date(), 
              "\n\n\n", sep=""))      
#
    ModelK[[k]][[1]] <- ModelK[[k]][[1]][Keep]
    ModelK[[k]][[2]] <- ModelK[[k]][[2]][Keep]  
    class(ModelK[[k]][[1]]) <- class(ModelK[[k]][[2]]) <- "NMixMCMC"
  }
}  
@ 
%$


% --------------------------------------
%\clearpage
\Ritem Summary of PED and DIC's for the fitted models:
<<Fixed K: PED and DIC>>=
PED <- ModelK[[1]]$PED
DIC <- list(Chain1=ModelK[[1]][[1]]$DIC, Chain2=ModelK[[1]][[2]]$DIC)
for (k in 2:length(ModelK)){
  PED <- rbind(PED, ModelK[[k]]$PED)
  DIC[[1]] <- rbind(DIC[[1]], ModelK[[k]][[1]]$DIC)
  DIC[[2]] <- rbind(DIC[[2]], ModelK[[k]][[2]]$DIC)  
}
rownames(PED) <- rownames(DIC[[1]]) <- rownames(DIC[[2]]) <- paste("K=", 1:length(ModelK), sep="")
@
<<Fixed K: print PED>>=
print(PED)
@ 
<<Fixed K: print DIC>>=
print(DIC)
@ 


% --------------------------------------
\Ritem Plot of the marginal predictive densities (\Rvar{eruptions}) for different values of $K$,
densities computed from chain~1 (see Figures~\ref{fig11} and \ref{fig12}):
<<Fixed K MCMC: Eruptions - plot of the marginal predictive density, keep.source=TRUE, results=hide>>=
CH <- 1
postscript(paste(FIGDIR, "figFaithful11.ps", sep=""), width=6, height=6, 
           horizontal=FALSE)
par(mfrow=c(1, 1), bty="n")
hist(Faithful$eruptions, prob=TRUE, col="grey90", 
     xlab="Eruptions (min)", ylab="Density", main="", 
     breaks=seq(1.4, 5.6, by=0.3))
for (k in 1:10){
  lines(MPDensModelK[[k]][[CH]]$x$x1, MPDensModelK[[k]][[CH]]$dens[[1]], 
        col="red")
}
dev.off()
@ 

% --------------------------------------
<<Fixed K MCMC: Eruptions - plot of the marginal predictive density, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful12.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(mar=c(3, 2, 2, 1)+0.1)
par(mfrow=c(5, 2), bty="n")
for (k in 1:10){
  hist(Faithful$eruptions, prob=TRUE, col="lightblue", 
       xlab="", ylab="", main=paste("K = ", k, sep=""), 
       breaks=seq(1.4, 5.6, by=0.3))
  lines(MPDensModelK[[k]][[CH]]$x$x1, MPDensModelK[[k]][[CH]]$dens[[1]], 
        col="red", lwd=2)
}
dev.off()
@ 

% --------------------------------------
\newpage
\Ritem Plot of the marginal predictive densities (\Rvar{waiting}) for different values of $K$,
densities computed from chain~1 (see Figures~\ref{fig13} and \ref{fig14}):
<<Fixed K MCMC: Waiting - plot of the marginal predictive density, keep.source=TRUE, results=hide>>=
CH <- 1
postscript(paste(FIGDIR, "figFaithful13.ps", sep=""), width=6, height=6, 
           horizontal=FALSE)
par(mfrow=c(1, 1), bty="n")
hist(Faithful$waiting, prob=TRUE, col="grey90", 
     xlab="Waiting (min)", ylab="Density", main="")
for (k in 1:10){
  lines(MPDensModelK[[k]][[CH]]$x$x2, MPDensModelK[[k]][[CH]]$dens[[2]], 
        col="red")
}
dev.off()
@ 

% --------------------------------------
<<Fixed K MCMC: Waiting - plot of the marginal predictive density, keep.source=TRUE, results=hide>>=
postscript(paste(FIGDIR, "figFaithful14.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(mar=c(3, 2, 2, 1)+0.1)
par(mfrow=c(5, 2), bty="n")
for (k in 1:10){
  hist(Faithful$waiting, prob=TRUE, col="lightblue", 
       xlab="", ylab="", main=paste("K = ", k, sep=""))
  lines(MPDensModelK[[k]][[CH]]$x$x2, MPDensModelK[[k]][[CH]]$dens[[2]], 
        col="red", lwd=2)
}  
dev.off()
@

% --------------------------------------
\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.45\textheight]{\FIGDIR/figFaithful11.ps}
\caption{Marginal predictive density for \Rvar{eruptions} based on the models with a~fixed number
of mixture components, results from chain~1.}\label{fig11}
\includegraphics[width=1.0\textwidth, height=0.45\textheight]{\FIGDIR/figFaithful13.ps}
\caption{Marginal predictive density for \Rvar{waiting} based on the models with a~fixed number 
of mixture components, results from chain~1.}\label{fig13}
\end{figure}

% --------------------------------------
\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful12.ps}
\caption{Marginal predictive density for \Rvar{eruptions} based on the models with a~fixed number
of mixture components, results from chain~1.}\label{fig12}
\end{figure}

% --------------------------------------
\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful14.ps}
\caption{Marginal predictive density for \Rvar{waiting} based on the models with a~fixed number
of mixture components, results from chain~1..}\label{fig14}
\end{figure}

% --------------------------------------
\clearpage
\Ritem Plots of the joint predictive density (contour plots) for different values of $K$,
densities computed from chain~1 (see Figure~\ref{fig15}):
<<Fixed K MCMC: Waiting - plot of the joint predictive density, keep.source=TRUE, results=hide>>=
CH <- 1
postscript(paste(FIGDIR, "figFaithful15.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(mar=c(3, 2, 2, 1)+0.1)
par(mfrow=c(5, 2), bty="n")
for (k in 1:10){
  plot(Faithful, col="red", xlab="", ylab="", 
       main=paste("K = ", k, sep=""))
  contour(JPDensModelK[[k]][[CH]]$x$x1, JPDensModelK[[k]][[CH]]$x$x2, 
          JPDensModelK[[k]][[CH]]$dens[["1-2"]], 
          col="darkblue", add=TRUE)
}  
dev.off()
@
%$

\Ritem Plots of the joint predictive density (image plots) for different values of $K$,
densities computed from chain~1 (see Figure~\ref{fig19}):
<<Fixed K MCMC: Waiting - plot of the joint predictive density, keep.source=TRUE, results=hide>>=
CH <- 1
postscript(paste(FIGDIR, "figFaithful19.ps", sep=""), width=7, height=10, 
           horizontal=FALSE)
par(mar=c(3, 2, 2, 1)+0.1)
par(mfrow=c(5, 2), bty="n")
for (k in 1:10){
  plot(Faithful, col="darkblue", xlab="", ylab="", 
       main=paste("K = ", k, sep=""))
  image(JPDensModelK[[k]][[CH]]$x$x1, JPDensModelK[[k]][[CH]]$x$x2, 
        JPDensModelK[[k]][[CH]]$dens[["1-2"]], add=TRUE,
        col=rev(heat_hcl(33, c=c(80, 30), l=c(30, 90), power=c(1/5, 1.3))))  
#  points(Faithful, col="darkblue")
}  
dev.off()
@
%$

% --------------------------------------
\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful15.ps}
\caption{Joint predictive density based on the models with a~fixed number of mixture components,
results from chain~1.}\label{fig15}
\end{figure}

% --------------------------------------
\begin{figure}[p!]\centering
\includegraphics[width=1.0\textwidth, height=0.9\textheight]{\FIGDIR/figFaithful19.ps}
\caption{Joint predictive density based on the models with a~fixed number of mixture components,
results from chain~1.}\label{fig19}
\end{figure}

% --------------------------------------
\clearpage
\Ritem Save results for future use:
<<Save results, keep.source=TRUE>>=
if (RUN.TIMECONSUMING.CODE){
  save(list="Model0", 
       file=paste(RESULTDIR, "/Faithful-Model0", Kshow, ".RData", sep=""))        
  save(list=c("ModelK", "MPDensModelK", "JPDensModelK"), 
       file=paste(RESULT2DIR, "/Faithful-Result.RData", sep=""))    
}  
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\bibliographystyle{akplainnat}\bibliography{mixAK}


\end{document}
