\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}
\date{\today}
\maketitle
\begin{abstract}
The \pkg{rebmix} package provides R functions for random univariate and multivariate finite mixture model generation, estimation, clustering and classification. Variables can be continuous, discrete, independent or dependent and may follow normal, lognormal, Weibull, gamma, binomial, Poisson or Dirac parametric families.
\end{abstract}
\section{Introduction}
To cite the REBMIX algorithm please refer to \citep{Nagode_Fajdiga_2011a, Nagode_Fajdiga_2011b, Nagode_2015}. For theoretical backgrounds please upload also \url{http://doi.org/10.5963/JAO0302001}.
\section{What's new in version \texttt{2.9.1}}
This version is further debugged version 2.8.4. The R code is extended and rewritten in S4 class system. The background C code is extended and rewritten as object-oriented C\texttt{++} code, too. The package can easier be extended to other parametric families. Multivariate normal mixtures with unrestricted variance-covariance matrices are added. Clustering is added and classification is 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 <- list(pdf1 = "gamma",
  theta1.1 = c(1/100, 1/100, 1/100, 1/100),
  theta2.1 = c(200, 400, 600, 800))

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

n <- c(40, 360)

Theta <- list(pdf1 = "gamma",
  theta1.1 = c(1/27, 1/270),
  theta2.1 = c(9, 90))

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

n <- c(80, 240, 80)

Theta <- list(pdf1 = "gamma",
  theta1.1 = c(1/20, 1, 1/20),
  theta2.1 = c(40, 6, 200))

gamma3 <- RNGMIX(Dataset.name = "gamma3", rseed = -4, n = n, 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 = gamma1@Dataset,
  Preprocessing = "Parzen window",
  cmax = 8,
  Criterion = c("AIC", "BIC"),
  pdf = "gamma")

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

gamma3est <- REBMIX(Dataset = gamma3@Dataset,
  Preprocessing = "histogram",
  cmax = 8,
  Criterion = "BIC",
  pdf = "gamma",
  K = 23:27)
@
\subsubsection{Plot method}
\begin{figure}[h]\centering
<<gamma3-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 2.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 and coef methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(gamma2est)

coef(gamma1est, pos = 2)
@
\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 <- list(pdf1 = rep("Poisson", 2),
  theta1.1 = c(3, 2),
  theta2.1 = c(NA, NA),
  pdf2 = rep("Poisson", 2),
  theta1.2 = c(9, 10),
  theta2.2 = c(NA, NA),
  pdf3 = rep("Poisson", 2),
  theta1.3 = c(15, 16),
  theta2.3 = c(NA, NA))

poisson <- RNGMIX(Dataset.name = paste("Poisson_", 1:10, sep = ""), n = n, 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 = poisson@Dataset,
  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 = 4.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}
\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 = poisson@Zt)

plot(poissonclu)
@
\caption{Poisson dataset. Predictive cluster membership (coloured circles), error (black circles).}
\end{figure}
\subsubsection{Summary and coef methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
## Visualize results.

summary(poissonest)

coef(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 and coef methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(wreathest)

coef(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}
<<rebmix-code, split = FALSE, results = hide, keep.source = FALSE>>=
## Estimate number of components, component weights and component parameters.

ex4.1est <- REBMIX(model = "REBMVNORM",
  Dataset = list(as.data.frame(ex4.1)),
  Preprocessing = "Parzen window",
  cmax = 10,
  Criterion = "AIC")
@
\subsubsection{Plot method}
\begin{figure}[h]\centering
<<ex4_1-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 3.0, 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}
\subsubsection{Clustering}
\begin{figure}[h]\centering
<<ex4_1-clu-fig, fig = TRUE, pdf = TRUE, png = FALSE, eps = FALSE, height = 3.0, 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}
\subsubsection{Summary method}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
summary(ex4.1est)
@
\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 = Iris@train,
  Preprocessing = "Parzen window",
  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 = Iris@test,
  Zt = Iris@Zt)
@
\subsubsection{Show and summary methods}
<<rebmix-code, split = FALSE, keep.source = FALSE>>=
iriscla

summary(iriscla)
@
\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 = chunk(Adult, i)@train,
    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 = Adult@test,
  Zt = Adult@Zt)
@
\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}
