\documentclass[11pt, oneside, a4paper]{article}

\usepackage[utf8]{inputenc}
\usepackage[sectionbib, round]{natbib}
\usepackage[textwidth = 17cm, top = 2cm, bottom = 2cm]{geometry}
\usepackage{nameref}
\usepackage{bm}
\usepackage{amsmath, amssymb, amsfonts}
\usepackage{graphicx}
\usepackage{hyperref}

\newcommand{\pkg}[1]{\texorpdfstring%
{{\normalfont\fontseries{b}\selectfont #1}}%
{#1}}

\begin{document}
\SweaveOpts{engine = R}
%\VignetteIndexEntry{rebmix: The Rebmix Package}
%\VignetteKeyword{continuous variable}
%\VignetteKeyword{discrete variable}
%\VignetteKeyword{finite mixture}
%\VignetteKeyword{parameter estimation}
%\VignetteKeyword{REBMIX algorithm}
%\VignettePackage{rebmix}
\bibliographystyle{abbrvnat}
\title{\pkg{rebmix}: Finite Mixture Modeling, Clustering \& Classification}
\author{Marko Nagode \& Branislav Pani\'{c}}
\date{\today}
\maketitle

\begin{abstract}
The \pkg{rebmix} package provides R functions for random univariate and multivariate finite mixture model generation, estimation, clustering, latent class analysis and classification. Variables can be continuous, discrete, independent or dependent and may follow normal, lognormal, Weibull, gamma, binomial, Poisson, Dirac or von Mises parametric families.
\end{abstract}

\section{Introduction}
To cite the REBMIX algorithm please refer to \citep{Nagode_Fajdiga_2011a, Nagode_Fajdiga_2011b, Nagode_2015, Nagode_2018}. For theoretical backgrounds please upload \url{http://doi.org/10.5963/JAO0302001}.

\section{What's new in version \texttt{2.11.0}}

Version \texttt{2.11.0} introduces the Expectation-Maximization (EM) algorithm for the improved estimation of Gaussian mixture model parameters (with diagonal and unrestriced covariance matrices). Here the REBMIX algorithm is used to asses the initial parameters of the EM algorithm. Two different variants of the EM algorithm are implemented, namely the original EM algorithm from \citep{Dempster_1977} and a $k$-means like variant of the EM algorithm (Classification EM) as described in \citep{Celeux_1992}. As the REBMIX algorithm estimates a wide range of parameters for the Gaussian mixture model for different numbers of components, three different strategies, named \textbf{exhaustive}, \textbf{best} and \textbf{single} have been implemented. The \textbf{exhaustive} strategy is used to run the EM algorithm (or variant) on each solution of Gaussian mixture model parameters provided by the REBMIX algorithm. The \textbf{best} strategy utilizes a voting scheme for the estimated parameters from the REBMIX algorithm and runs the EM algorithm only on selected optimal parameters. Best candidates are chosen based on the value of the likelihood function (the highest one) for each number of components $c$ from a minimum specified \texttt{cmin} to a maximum specified \texttt{cmax}. The \textbf{single} strategy is useful when the single value of, for example, the number of bins in histogram preprocessing is supplied as input for the REBMIX algorithm. Otherwise, when multiple numbers of bins $k$ are supplied, this strategy is the same as the \textbf{exhaustive} strategy. To tackle the slow linear convergence of the EM algorithm, simple acceleration methods are implemented which can be controled with parameter \texttt{acceleration} and \texttt{acceleration.multiplier}. The increment of the EM algorithm in each iteration can be written as

\begin{equation}
\varDelta\boldsymbol{\Theta} = \boldsymbol{\Theta}^{(i+1)} - \boldsymbol{\Theta}^{(i)}
\end{equation}

Insted of using a standard EM increment $\varDelta\boldsymbol{\Theta}$ to reduce the number of iterations needed for the EM algorithm, this increment can be multiplied with some multiplier $a_{\textmd{EM}}$ which is referred to as \texttt{acceleration.multiplier}. Therefore the update in each EM iteration now becomes

\begin{equation}
\boldsymbol{\Theta}^{i+1} = \boldsymbol{\Theta}^{(i)} + a_{\textmd{EM}}\varDelta\boldsymbol{\Theta}
\end{equation}

The safe range for the $a_{\textmd{EM}}$ multiplier lies between 1.0 and 2.0, where 1.0 gives a standard EM increment and 2.0 doubles the EM increment. However, this does not necessarily mean that multiplication by a value of 2.0 will double the speed of the EM algorithm (i.e. by reducing the required number of iterations by 2). Here, 1.5 is a safe value which mostly speeds up the EM algorithm whilst retaining good results for the estimated parameters. A value of 1.9 can significantly speed up the estimation process, yet it can also deteriorate the quality of the resulting estimated parameters. Therefore, the value of the multiplicator needs to be set carefully. This value is set with \texttt{acceleration.multiplier} parameter. The other parameter \texttt{acceleration} controls how the $a_{\textmd{EM}}$ multiplier is handled and can be one of \textbf{fixed}, \textbf{line} and \textbf{golden}. Selecting the \textbf{fixed} option means that the $a_{\textmd{EM}}$ multiplier is specified via the \texttt{acceleration.multiplier} parameter and for each iteration of the EM algorithm the increment is increased by a specified value of $a_{\textmd{EM}}$. The \textbf{line} and \textbf{golden} options perform a, line and golden search (respectively) for the optimal value of $a_{\textmd{EM}}$ for which the highest increase in the likelihood function of each EM iteration is achieved.

EM handling is  carried out using the newly introduced class \texttt{"EM.Control"}. Classes \texttt{"REBMIX"} and \texttt{"REBMVNORM"} and its signature method \texttt{REBMIX} now accept the \texttt{"EM.Control"} object via the argument called \texttt{"EMcontrol"}. The class \texttt{EM.Control} has the same name convection for slots as the input argument \texttt{EMcontrol} (\texttt{strategy}, \texttt{variant}, \texttt{acceleration}, \texttt{tolerance}, \texttt{acceleration.multiplier} and \texttt{maximum.iterations}) as well as all accessor functions with the same name convention as \texttt{a.\textit{slot name}} and setter function \texttt{a.\textit{slot name}<-}.

\section{Previous versions}
Methods \texttt{Zp} and \texttt{coef} have been replaced by \texttt{a.Zp}, \texttt{a.theta1.all} and \texttt{a.theta2.all} getters. All slots can be accessed via accessors. Their names are  generally composed of \texttt{a.} followed by the slot name and are used to read the slots. Class \texttt{"RNGMIX.Theta"} has been added to simplify random finite mixture model generation. Method \texttt{show} has been added for \texttt{"RCLS.chunk"} class. The minimum number of components \texttt{cmin} was added to \texttt{REBMIX} arguments and to the \texttt{"REBMIX"} class.  The \texttt{"Parzen window"} preprocessing has been renamed to more commonly known \texttt{"kernel density estimation"}. Rough parameter estimation for binomial and Poisson parametric families has also been improved and the package is now broadened to latent class analysis in version 2.10.3. Method \texttt{split} has been improved and examples for its proper use are added.

GCC 8.1 notes and warnings in C++ functions have been eliminated in version 2.10.2. Cholesky decomposition is now used to calculate the logarithm of the determinant and inverse of variance-covariance matrices instead of LU decomposition. Special attention has been paid to resolving numerical problems related to high dimensional datasets.

Version 2.10.1 is the further debugged version of 2.10.0. Large \texttt{K} in combination with large dimension $d$ can lead to histograms with numerous nonempty bins $v$. In order to restrain $v$, the well known RootN rule \citep{Velleman_1976} may intuitively be extended to multidimensions
\begin{equation}
v_{\mathrm{max}} = \frac{1 + d}{d} n^\frac{d}{1 + d}.
\end{equation}
If $d = \infty$, then $v_{\mathrm{max}} = n$. If $d = 1$, then $v_{\mathrm{max}} = 2 \sqrt{n}$. Minor debugging and function improvements have also been carried out in version 2.10.0. The acceleration rate is now progressively increasing. Each time the inner loop starts, the counter $I_{2}$ \citep[see][for details]{Nagode_2015} is initiated and constant
\begin{equation}
A = \left. \frac{1 - a_{\mathrm{r}}}{a_{\mathrm{r}} (D_{l} w_{l} - D_{\mathrm{min}})} \right|_{I_{2} = 1}
\end{equation}
is calculated. The acceleration rate $a_{\mathrm{r}}$ at $I_{2} = 1$ always equals the value stored in the input argument \texttt{ar}. Otherwise
\begin{equation}
a_{\mathrm{r}} = \left. \frac{1}{A (D_{l} w_{l} - D_{\mathrm{min}}) + 1} \right|_{I_{2} > 1}.
\end{equation}
The Newton-Raphson root finding in C\texttt{++} functions was improved in version 2.9.3. This affects only Weibull, gamma and von Mises parametric families. A circular von Mises parametric family has been added and further debugging carried out in version 2.9.2. Version 2.9.1 is a further debugged version 2.8.4. The R code has been extended and rewritten in S4 class system. The background C code has also been extended and rewritten as object-oriented C\texttt{++} code. The package can now more easily be extended to other parametric families. Multivariate normal mixtures with unrestricted variance-covariance matrices have been added. Clustering has also been added and classification improved.

\section{Examples}
To illustrate the use of the REBMIX algorithm, univariate and multivariate datasets are considered. The \pkg{rebmix} is loaded and the prompt before starting new page is set to \texttt{TRUE}.
<<rebmix-code, split = FALSE, echo = FALSE, keep.source = FALSE>>=
##############################################
## R sources for reproducing the results in ##
##   Marko Nagode:                          ##
##   rebmix: The Rebmix Package             ##
##############################################

options(prompt = "R> ", continue = "+  ", width = 80,
  useFancyQuotes = FALSE, digits = 3)
@
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
###################
## Preliminaries ##
###################

## load package and set prompt before starting new page to TRUE.

library("rebmix")
devAskNewPage(ask = TRUE)
@

\subsection{Gamma datasets}
Three gamma mixtures are considered \citep{Wiper_2001}. The first has four well-separated components with means $2$, $4$, $6$ and $8$, respectively
\begin{center}
\(\begin{array}{lll}
\theta_{1} = 1/100 & \beta_{1} = 200 & n_{1} = 100 \\
\theta_{2} = 1/100 & \beta_{2} = 400 & n_{2} = 100 \\
\theta_{3} = 1/100 & \beta_{3} = 600 & n_{3} = 100 \\
\theta_{4} = 1/100 & \beta_{4} = 800 & n_{4} = 100.
\end{array}\)
\end{center}
The second has equal means but different variances and weights
\begin{center}
\(\begin{array}{lll}
\theta_{1} = 1/27 & \beta_{1} = 9 & n_{1} = 40 \\
\theta_{2} = 1/270 & \beta_{2} = 90 & n_{2} = 360.
\end{array}\)
\end{center}
The third is a mixture of a rather diffuse component with mean $6$ and two lower weighted components with smaller variances and means of $2$ and $10$, respectively
\begin{center}
\(\begin{array}{lll}
\theta_{1} = 1/20 & \beta_{1} = 40 & n_{1} = 80 \\
\theta_{2} = 1 & \beta_{2} = 6 & n_{2} = 240 \\
\theta_{3} = 1/20 & \beta_{3} = 200 & n_{3} = 80.
\end{array}\)
\end{center}

\subsubsection{Finite mixture generation}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
######################
##  Gamma datasets  ##
######################

## Generate gamma datasets.

n <- c(100, 100, 100, 100)

Theta <- new("RNGMIX.Theta", c = 4, pdf = "gamma")

a.theta1(Theta) <- rep(1/100, 4)
a.theta2(Theta) <- c(200, 400, 600, 800)

gamma1 <- RNGMIX(Dataset.name = "gamma1", n = n, Theta = a.Theta(Theta))

n <- c(40, 360)

Theta <- new("RNGMIX.Theta", c = 2, pdf = "gamma")

a.theta1(Theta) <- c(1/27, 1 / 270)
a.theta2(Theta) <- c(9, 90)

gamma2 <- RNGMIX(Dataset.name = "gamma2", n = n, Theta = a.Theta(Theta))

n <- c(80, 240, 80)

Theta <- new("RNGMIX.Theta", c = 3, pdf = "gamma")

a.theta1(Theta) <- c(1/20, 1, 1/20)
a.theta2(Theta) <- c(40, 6, 200)

gamma3 <- RNGMIX(Dataset.name = "gamma3", rseed = -4, n = n, Theta = a.Theta(Theta))
@

\subsubsection{Finite mixture estimation}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
## Estimate number of components, component weights and component parameters.

gamma1est <- REBMIX(Dataset = a.Dataset(gamma1),
  Preprocessing = "kernel density estimation",
  cmax = 8,
  Criterion = "BIC",
  pdf = "gamma")

gamma2est <- REBMIX(Dataset = a.Dataset(gamma2),
  Preprocessing = "histogram",
  cmax = 8,
  Criterion = "BIC",
  pdf = "gamma")

gamma3est <- REBMIX(Dataset = a.Dataset(gamma3),
  Preprocessing = "histogram",
  cmax = 8,
  Criterion = "BIC",
  pdf = "gamma",
  K = 23:27)
@

\subsubsection{Plot method}
\begin{figure}[htb]\centering
<<gamma3-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 3.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
plot(gamma3est, pos = 1, what = c("den", "dis"), ncol = 2, npts = 1000)
@
\caption{Gamma 3 dataset. Empirical density (circles) and predictive gamma mixture density in black solid line.}
\end{figure}

\subsubsection{Summary, a.theta1.all and a.theta2.all methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(gamma2est)

a.theta1.all(gamma1est, pos = 1)

a.theta2.all(gamma1est, pos = 1)
@

\subsubsection{Bootstrap methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
## Bootstrap finite mixture.

gamma3boot <- boot(x = gamma3est, pos = 1, Bootstrap = "p", B = 10)

gamma3boot

summary(gamma3boot)
@

\subsection{Poisson dataset}
Dataset consists of $n = 600$ two~dimensional observations obtained by generating data points separately from each of three Poisson distributions. The component~dataset sizes and parameters, which are those studied in \citet{Jinwen_2009}, are displayed below
\begin{center}
\(\begin{array}{ll}
\bm{\theta}_{1} = (3, 2)^{\top} & n_{1} = 200 \\
\bm{\theta}_{2} = (9, 10)^{\top} & n_{2} = 200 \\
\bm{\theta}_{3} = (15, 16)^{\top} & n_{3} = 200
\end{array}\)
\end{center}
For the dataset \citet{Jinwen_2009} conduct $100$ experiments by selecting different initial values of the mixing proportions. In all the cases, the adaptive gradient BYY learning algorithm leads to the correct model selection, i.e., finally allocating the correct number of Poissons for the dataset. In the meantime, it also results in an estimate for each parameter in the original or true Poisson mixture which
generated the dataset. As the dataset of \citet{Jinwen_2009} can not exactly be reproduced, $10$ datasets are generated with random seeds $r_{\mathrm{seed}}$ ranging from $-1$ to $-10$.

\subsubsection{Finite mixture generation}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
#########################
##   Poisson dataset   ##
#########################

## Generate the Poisson dataset.

n <- c(200, 200, 200)

Theta <- new("RNGMIX.Theta", c = 3, pdf = rep("Poisson", 2))

a.theta1(Theta, 1) <- c(3, 2)
a.theta1(Theta, 2) <- c(9, 10)
a.theta1(Theta, 3) <- c(15, 16)

poisson <- RNGMIX(Dataset.name = paste("Poisson_", 1:10, sep = ""), n = n, Theta = a.Theta(Theta))
@

\subsubsection{Finite mixture estimation}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
## Estimate number of components, component weights and component parameters.

poissonest <- REBMIX(Dataset = a.Dataset(poisson),
  Preprocessing = "histogram",
  cmax = 10,
  Criterion = "MDL5",
  pdf = rep("Poisson", 2),
  K = 1)
@

\subsubsection{Plot method}
\begin{figure}[h]\centering
<<poisson-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 6.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
plot(poissonest, pos = 9, what = c("dens", "marg", "IC", "D", "logL"), nrow = 2, ncol = 3, npts = 1000)
@
\caption{Poisson dataset. Empirical densities (coloured large circles), predictive multivariate Poisson-Poisson mixture density (coloured small circles), empirical densities (circles), predictive univariate marginal Poisson mixture densities and progress charts (solid line).}
\end{figure}

\pagebreak
\subsubsection{Clustering}
\begin{figure}[h]\centering
<<poisson-clu-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 4.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
poissonclu <- RCLRMIX(x = poissonest, pos = 9, Zt = a.Zt(poisson))

plot(poissonclu)
@
\caption{Poisson dataset. Predictive cluster membership (coloured circles), error (black circles).}
\end{figure}

\subsubsection{Summary, a.theta1.all and a.theta2.all methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
## Visualize results.

summary(poissonest)

a.theta1.all(poissonest, pos = 9)

a.theta2.all(poissonest, pos = 9)
@

\subsection{Multivariate normal \texttt{wreath} dataset}
A \texttt{wreath} dataset \citep{Fraley_2005} consist of $1000$ observations drawn from a $14$-component normal mixture in which the covariances of the components have the same size and shape but differ in orientation.
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
data("wreath", package = "mclust")
@

\subsubsection{Finite mixture estimation}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
## Estimate number of components, component weights and component parameters.

wreathest <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(wreath)),
  Preprocessing = "histogram",
  cmax = 20,
  Criterion = "BIC")
@

\subsubsection{Plot method}
\begin{figure}[h]\centering
<<wreath-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 5.0, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
plot(wreathest)
@
\caption{Dataset \texttt{wreath}. Empirical densities (coloured circles), predictive multivariate normal mixture density (coloured lines).}
\end{figure}

\subsubsection{Clustering}
\begin{figure}[h]\centering
<<wreath-clu-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 4.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
wreathclu <- RCLRMIX(model = "RCLRMVNORM", x = wreathest)

plot(wreathclu, s = 14)
@
\caption{Dataset \texttt{wreath}. Predictive cluster membership (coloured circles).}
\end{figure}

\subsubsection{Summary method}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(wreathest)
@

\subsubsection{Summary method}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(wreathclu)
@

\subsection{Multivariate normal \texttt{ex4.1} dataset}
A \texttt{ex4.1} dataset \citep{Baudry_2010, mclust_20__} consist of $600$ two~dimensional observations.
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
data("Baudry_etal_2010_JCGS_examples", package = "mclust")
@

\subsubsection{Finite mixture estimation using exhaustive REBMIX\&EM strategy}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
## Estimate number of components, component weights and component parameters.

## Create the EM.Control object to utilize one of the REBMIX&EM strategies

EM <- new("EM.Control",
  strategy = "exhaustive",
  variant = "EM",
  acceleration = "fixed",
  acceleration.multiplier = 1.0,
  tolerance = 1e-4,
  maximum.iterations = 1000)

ex4.1est.dens <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(ex4.1)),
  Preprocessing = "histogram",
  cmax = 10,
  Criterion = "AIC",
  EMcontrol = EM)
@

\subsubsection{Plot method}
\begin{figure}[h]\centering
<<ex4_1_dens-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 4.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
plot(ex4.1est.dens, pos = 1, what = c("dens"), nrow = 1, ncol = 1)
@
\caption{Dataset \texttt{ex4.1}. Empirical densities (coloured circles), predictive multivariate normal mixture density (coloured lines).}
\end{figure}

\subsubsection{Clustering}
\begin{figure}[h]\centering
<<ex4_1_dens-clu-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 4.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
ex4.1clu.dens <- RCLRMIX(model = "RCLRMVNORM", x = ex4.1est.dens)

plot(ex4.1clu.dens)
@
\caption{Dataset \texttt{ex4.1}. Predictive cluster membership (coloured circles).}
\end{figure}

\subsubsection{Summary method}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(ex4.1est.dens)
@

\subsubsection{Summary method}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(ex4.1clu.dens)
@

\subsubsection{Clustering with exhaustive REBMIX\&ECM strategy and ICL criterion}

<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
## Estimate number of components, component weights and component parameters.

## Create the EM.Control object to utilize one of the REBMIX&EM strategies

CEM <- new("EM.Control",
  strategy = "exhaustive",
  variant = "ECM",
  acceleration = "fixed",
  acceleration.multiplier = 1.0,
  tolerance = 1e-4,
  maximum.iterations = 1000)

ex4.1est <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(ex4.1)),
  Preprocessing = "histogram",
  cmax = 10,
  Criterion = "ICL",
  EMcontrol = CEM)

@
\begin{figure}[htb]\centering
<<ex4_1-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 5.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
plot(ex4.1est, pos = 1, what = c("dens"), nrow = 1, ncol = 1)
@
\caption{Dataset \texttt{ex4.1}. Empirical densities (coloured circles), predictive multivariate normal mixture density (coloured lines).}
\end{figure}

<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(ex4.1est)
@


\begin{figure}[htb]\centering
<<ex4_1-clu-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 5.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
ex4.1clu <- RCLRMIX(model = "RCLRMVNORM", x = ex4.1est)

plot(ex4.1clu)
@
\caption{Dataset \texttt{ex4.1}. Predictive cluster membership (coloured circles).}
\end{figure}

<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(ex4.1clu)
@


\subsubsection{Acceleration of EM algorithm}

Standard EM algorithm with fixed acceleration.multiplier of $a_{\texttt{EM}} = 1.0$:

<<rebmix-code, split = FALSE, keep.source = FALSE>>=
## Estimate number of components, component weights and component parameters.

## Create the EM.Control object to utilize one of the REBMIX&EM strategies

EM.normal <- new("EM.Control",
  strategy = "exhaustive",
  variant = "EM",
  acceleration = "fixed",
  acceleration.multiplier = 1.0,
  tolerance = 1e-4,
  maximum.iterations = 1000)

ex4.1est.em.normal <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(ex4.1)),
  Preprocessing = "histogram",
  cmax = 15,
  Criterion = "AIC",
  EMcontrol = EM.normal)

cat("Total number of EM algorithm iterations: ",
  a.summary.EM(ex4.1est.em.normal, pos=1, col.name = "total.iterations.nbr"),
  ". Value of AIC: ", a.summary(ex4.1est.em.normal, pos = 1, col.name = "IC"))
@

Standard EM algorithm with fixed acceleration.multiplier of $a_{\texttt{EM}} = 1.5$::

<<rebmix-code, split = FALSE, keep.source = FALSE>>=
EM.fixed1.5 <- new("EM.Control",
  strategy = "exhaustive",
  variant = "EM",
  acceleration = "fixed",
  acceleration.multiplier = 1.5,
  tolerance = 1e-4,
  maximum.iterations = 1000)

ex4.1est.em.fixed1.5 <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(ex4.1)),
  Preprocessing = "histogram",
  cmax = 15,
  Criterion = "AIC",
  EMcontrol = EM.fixed1.5)

cat("Total number of EM algorithm iterations: ",
  a.summary.EM(ex4.1est.em.fixed1.5, pos=1, col.name = "total.iterations.nbr"),
  ". Value of AIC: ", a.summary(ex4.1est.em.fixed1.5, pos = 1, col.name = "IC"))
@

Standard EM algorithm with line search for optimal increment $a_{\texttt{EM}}$ in each iteration:

<<rebmix-code, split = FALSE, keep.source = FALSE>>=
EM.line <- new("EM.Control",
  strategy = "exhaustive",
  variant = "EM",
  acceleration = "line",
  acceleration.multiplier = 1.0,
  tolerance = 1e-4,
  maximum.iterations = 1000)

ex4.1est.em.line <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(ex4.1)),
  Preprocessing = "histogram",
  cmax = 15,
  Criterion = "AIC",
  EMcontrol = EM.line)

cat("Total number of EM algorithm iterations: ",
  a.summary.EM(ex4.1est.em.line, pos=1, col.name = "total.iterations.nbr"),
  ". Value of AIC: ", a.summary(ex4.1est.em.line, pos = 1, col.name = "IC"))
@

Standard EM algorithm with golden search for optimal increment $a_{\texttt{EM}}$ in each iteration:

<<rebmix-code, split = FALSE, keep.source = FALSE>>=
EM.golden <- new("EM.Control",
  strategy = "exhaustive",
  variant = "EM",
  acceleration = "golden",
  acceleration.multiplier = 1.0,
  tolerance = 1e-4,
  maximum.iterations = 1000)

ex4.1est.em.golden <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(ex4.1)),
  Preprocessing = "histogram",
  cmax = 15,
  Criterion = "AIC",
  EMcontrol = EM.golden)

cat("Total number of EM algorithm iterations: ",
  a.summary.EM(ex4.1est.em.golden, pos=1, col.name = "total.iterations.nbr"),
  ". Value of AIC: ", a.summary(ex4.1est.em.golden, pos = 1, col.name = "IC"))
@

\subsection{Multivariate \texttt{iris} dataset}
The well known set of iris data as collected originally by \citet{Anderson_1936} and first analysed by \citet{Fisher_1936} is considered here. It is available at \citet{Asuncion_Newman_2007} consisting of the measurements of the length and width of both sepals and petals of $50$ plants for each of the three types of iris species setosa, versicolor and virginica.  The iris dataset is loaded, split into three subsets for the three classes and the \texttt{Class} column is removed.
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
data("iris")

# Show level attributes discrete variables.

levels(iris[["Class"]])

# Split dataset into train (75%) and test (25%) subsets.

set.seed(5)

Iris <- split(p = 0.75, Dataset = iris, class = 5)
@

\subsubsection{Finite mixture estimation}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
# Estimate number of components, component weights and component
# parameters for train subsets.

irisest <- REBMIX(model = "REBMVNORM",
  Dataset = a.train(Iris),
  Preprocessing = "kernel density estimation",
  cmax = 10,
  Criterion = "ICL-BIC")
@

\subsubsection{Classification}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
# Selected features.

iriscla <- RCLSMIX(model = "RCLSMVNORM",
  x = list(irisest),
  Dataset = a.test(Iris),
  Zt = a.Zt(Iris))
@

\subsubsection{Show and summary methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
iriscla

summary(iriscla)
@
\pagebreak
\subsubsection{Plot method}
\begin{figure}[h]\centering
<<iris-cla-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 5.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
# Plot selected features.

plot(iriscla, nrow = 3, ncol = 2)
@
\caption{Dataset \texttt{iris}. Predictive class membership (coloured circles), error (black circles).}
\end{figure}

\subsection{Multivariate \texttt{adult} dataset}
The \texttt{adult} dataset containing $48842$ instances with $16$ continuous, binary and discrete variables was extracted from the census bureau database \citet{Asuncion_Newman_2007}. Extraction was done by Barry Becker from the 1994 census bureau database. The \texttt{adult} dataset is loaded, complete cases are extracted and levels are replaced with numbers.
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
data("adult")

# Find complete cases.

adult <- adult[complete.cases(adult),]

# Replace levels with numbers.

adult <- as.data.frame(data.matrix(adult))
@
Numbers of unique values for variables are determined and displayed.
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
# Find numbers of levels.

cmax <- unlist(lapply(apply(adult[, c(-1, -16)], 2, unique), length))

cmax
@
The  dataset is split into train and test subsets for the two incomes and the \texttt{Type} and \texttt{Income} columns are removed.
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
# Split adult dataset into train and test subsets for two Incomes
# and remove Type and Income columns.

Adult <- split(p = list(type = 1, train = 2, test = 1),
  Dataset = adult, class = 16)
@

\subsubsection{Finite mixture estimation}
Number of components, component weights and component parameters are estimated assuming that the variables are independent for the set of chunks $y_{1j}, y_{2j}, \ldots, y_{14j}$.
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
# Estimate number of components, component weights and component parameters
# for the set of chunks 1:14.

adultest <- list()

for (i in 1:14) {
  adultest[[i]] <- REBMIX(Dataset = a.train(chunk(Adult, i)),
    Preprocessing = "histogram",
    cmax = min(120, cmax[i]),
    Criterion = "BIC",
    pdf = "Dirac",
    K = 1)
}
@

\subsubsection{Classification}
The class membership prediction is based upon the best first search algorithm.
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
# Class membership prediction based upon the best first search algorithm.

adultcla <- BFSMIX(x = adultest,
  Dataset = a.test(Adult),
  Zt = a.Zt(Adult))
@

\subsubsection{Show and summary methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
adultcla

summary(adultcla)
@

\subsubsection{Plot method}
\begin{figure}[h]\centering
<<adult-cla-fig, fig = TRUE, pdf = FALSE, png = TRUE, eps = FALSE, height = 5.5, width = 5.5, echo = TRUE, results = hide, keep.source = FALSE>>=
# Plot selected chunks.

plot(adultcla, nrow = 5, ncol = 2)
@
\caption{Dataset \texttt{adult}. Predictive class membership (coloured circles), error (black circles).}
\end{figure}
<<rebmix-code, split = FALSE, echo = FALSE, results = hide, keep.source = FALSE>>=
rm(list = ls())
@

\section{Summary}\label{sec:summary}
The users of the \texttt{rebmix} package are kindly encouraged to inform the author about bugs and wishes.

\bibliography{rebmix}
\vspace{\baselineskip}\noindent\emph{Marko Nagode\\
University of Ljubljana\\
Faculty of Mechanical Engineering\\
A\v{s}ker\v{c}eva 6\\
1000 Ljubljana\\
Slovenia}\\
\href{mailto:Marko.Nagode@fs.uni-lj.si}{Marko.Nagode@fs.uni-lj.si}.
\end{document}
