Multivariate normal mixture modeling, clustering and classification with the rebmix package
Marko Nagode
stat.ML cs.LG 62H30, 62F10














exampleVerbatim
example*




=[1]#1








=


abbrvnat




Multivariate Normal Mixture Modeling, Clustering and Classification with the rebmix Package
    Marko Nagode
    December 30, 2023
===========================================================================================




The rebmix package provides R functions for random univariate and multivariate finite mixture model generation, estimation, clustering and classification. The paper is focused on multivariate normal mixture models with unrestricted variance-covariance matrices. The objective is to show how to generate datasets for a known number of components, numbers of observations and component parameters, how to estimate the number of components, component weights and component parameters and how to predict cluster and class membership based upon a model trained by the REBMIX algorithm. The accompanying plotting, bootstrapping and other features of the package are dealt with, too. For demonstration purpose a multivariate normal dataset with unrestricted variance-covariance matrices is studied.




§ INTRODUCTION

Mixture models are a commonly employed tool in statistical modeling. There are quite some packages in the R statistical environment enabling number of components, component weights and component parameter estimation, clustering, classification or all <cit.>. The normal mixture models are implemented, e.g., in the following packages listed in the chronological order: mclust <cit.>, flexmix <cit.>, mixtools <cit.>, HDclassif <cit.>, bgmm <cit.>, EMCluster <cit.>, Rmixmod <cit.>, mixture <cit.> and GMCM <cit.>. Some of them can also handle non-normal components and they all rely on the expectation-maximization (EM) algorithm or on one of its variants.

The paper aims to present the rebmix package. REBMIX stands for the Rough and Enhanced component parameter estimation that is followed by the Bayesian classification of the remaining observations for the finite MIXture estimation. REBMIX originating in <cit.> is an alternative algorithm to parameter estimation for finite mixture models. It is based on the hypothesis that the problem of finite mixture estimation can be broken into multiple problems of component parameter estimations for basic parametric family types. Here the governing equations preexist and originate in the maximum likelihood. REBMIX thus estimates component weight and component parameters for each component individually as distinguished from EM where component weights and component parameters for all components are estimated simultaneously.

The rest of the paper is organized as follows. First, the REBMIX algorithm with its input and output arguments is presented. Second, the set of equations for multivariate normal mixture models with unrestricted variance-covariance matrices is derived. Next, a multivariate normal dataset is studied to demonstrate how the rebmix package can be applied in the area of finite mixture estimation, clustering and classification. The last section is the conclusion.



§ ALGORITHM

Let y_1, …, y_n be an observed d dimensional dataset of size n of continuous
vector observations y_j = (y_1j, …, y_ij, …, y_dj)^⊤. Each observation is assumed to follow the predictive mixture density

    f(y | c, w, ) = ∑_l = 1^c w_l f(y | _l)

with multivariate normal component densities

    f(y | _l) = 1/√((2 π)^dl)exp{-1/2 (y - μ_l)^⊤l^-1 (y - μ_l) }.

The objective is to obtain the number of components c, component weights w_l summing to 1 and component parameters _l = (μ_l, l)^⊤. The REBMIX algorithm for conditionally independent normal, lognormal, Weibull, gamma, binomial, Poisson and Dirac component densities is explained in detail in <cit.>. The paper is therefore focused on the mixtures of multivariate normal component densities with unrestricted variance-covariance matrices.



 §.§ Algorithm flow


[1]Mandatory argument.
The REBMIX algorithm implemented in R package rebmix is depicted in Fig. <ref>. It may be affected by maximum fourteen input arguments. Depending on the parametric families, three or four of them are mandatory while the rest are optional. For mixtures of conditionally independent component densities, argument model can be omitted. If model = "REBMVNORM", then arguments pdf and theta1 are superfluous and the output for mixtures of multivariate normal component densities with unrestricted variance-covariance matrices is returned. Dataset is a list of data frames of size n × d containing d dimensional datasets. Each of the d columns represents one random variable. Number of observations n equals the number of rows in the datasets. Preprocessing is a character vector giving one of the preprocessing types "histogram", "Parzen window" or "k-nearest neighbour". Maximum number of components c_max∈ℕ. The default value is cmax = 15. Criterion is a character vector giving one of the information criterion types default Akaike "AIC", "AIC3", "AIC4" or "AICc", Bayesian "BIC", consistent Akaike "CAIC", Hannan-Quinn "HQC", minimum description length "MDL2" or "MDL5", approximate weight of evidence "AWE", classification likelihood "CLC", integrated classification likelihood "ICL" or "ICL-BIC", partition coefficient "PC", total of positive relative deviations "D" or sum of squares error "SSE". Argument pdf is a character vector of length d composed of one of continuous or discrete parametric family types "normal", "lognormal", "Weibull", "gamma", "vonMises", "binomial", "Poisson" or "Dirac". Initial component parameters theta1 is a vector of length d and equals n_il = Number of categories - 1 for "binomial" distribution or "NA" otherwise. Argument theta2 is a vector of length d and is currently not used.

K is a vector or a list of vectors containing numbers of bins v for the histogram and the Parzen window or numbers of nearest
neighbours k for the k nearest neighbour. There is no genuine rule to identify v or k. Consequently, the algorithm identifies them from the set K of input values by minimizing the information criterion. The <cit.> rule v = 1 + log_2(n), Log_10 rule v = 10 log_10(n) or RootN
rule v = 2 √(n) can be applied to estimate the limiting numbers of bins or the rule of thumb k = √(n) to guess the intermediate number of nearest neighbours. If, e.g., K = c(10, 20, 40, 60) and minimum IC coincides with 40, brackets are set to 20 and 60 and the golden section is applied to refine the minimum search. The default value is "auto". Arguments y0, ymin and ymax are vectors of length d that hold origins, minimum and maximum observations respectively. The default values are numeric(). Argument ar stands for acceleration rate 0 < a_r≤ 1. The default value is 0.1 and does not have to be altered. Restraints is a character giving one of the restraints type "rigid" or default "loose". The rigid restraints are obsolete and applicable for well separated components only.

The algorithm is to a large extent independent of parametric families. This means that when it is extended to other parametric families, only the set of equations for rough and enhanced component parameter estimation, first and second moment calculation and Bayes classification of the unassigned observations have to be derived. The rest remains untouched.

An object of class "REBMIX" or "REBMVNORM" is returned by the REBMIX method. The object includes input and output arguments. Output Variables is a character vector of length d containing one of the types of variables "continuous" or "discrete". Argument w is a list of vectors of length c that holds component weights w_l summing to 1. Theta is a list of lists each containing c parametric family types pdfl. Each character vector pdfl is composed of "normal", "lognormal", "Weibull", "gamma", "vonMises", "binomial", "Poisson" or "Dirac" parametric family types. Component parameters theta1.l follow the parametric family types. Each vector theta1.l is composed of μ_il for normal, lognormal and von Mises distributions or θ_il for Weibull, gamma, binomial, Poisson and Dirac distributions. Component parameters theta2.l follow theta1.l. Each vector theta2.l is composed of σ_il for normal and lognormal distributions, β_il for Weibull and gamma distributions, p_il for binomial distribution and κ_il for von Mises distribution. Argument summary is a data frame with additional information about dataset, preprocessing, c_max, information criterion type, a_r, restraints type, optimal c, optimal v or k, K, y_i0, y_imin, y_imax, optimal bin widths h_i, information criterion IC, log likelihood log L and degrees of freedom M.

Output arguments opt.c, opt.IC, opt.logL and opt.D are lists of vectors that contain numbers of components, information criteria, log likelihoods and totals of positive relative deviations for optimal v for the histogram and the Parzen window or for optimal number of nearest neighbours k for the k nearest neighbour. Output all.K and all.IC are lists of vectors with all processed numbers of bins v for the histogram and the Parzen window or for all processed numbers of nearest neighbours k for the k nearest neighbour and the corresponding information criteria.



 §.§ Preprocessing of observations

The observations have to be preprocessed initially. In other words, empirical densities have to be assigned to the observations. For this purpose the histogram, Parzen window or k nearest neighbour may be used. Here number of bins v or nearest neighbours k plays an important role. To find the number of bins or nearest neighbours resulting in the lowest value of the information criterion, the loop 4 in Fig. <ref> runs for all K. Package KernSmooth may, e.g., be used to attain the optimal number of bins for d = 1.



 §.§ Global mode detection

Fundamental assumption of the REBMIX algorithm is that at least one component of the mixture should appear at the vicinity of the global mode, where for the Parzen window and k nearest neighbour the global mode stands for the d dimensional observation y_m with the highest empirical density. For the histogram the global mode corresponds to the mean of a bin y̅_m with the highest empirical density. Index m is used here to denote this particular value in the d dimensional space. The loops 3 and 2 in Fig. <ref> are executed iteratively. The algorithm presumes only one component initially and calculates the information criterion for a mixture with that component. Onwards the number of components increases gradually and the information criterion is calculated for the corresponding mixtures. The loops 3 and 2 stop when number of components c ≥ c_max or c ≥ v or c ≥ k.



 §.§ Rough component parameter estimation

Rough multivariate normal component parameter estimation is based on the component conditional densities

    f(y_i | y_î, _il) = 1/√(2 π)σ_ilexp{ -1/2(y_i - μ_il)^2)/σ_il^2},  i = 1, …, d.

If the equality of the predictive and empirical component conditional densities

    f(ŷ_im | ŷ_îm, _il) = f_i | î.lm

at the global mode ŷ_m is met, then

    μ_il = ŷ_im, σ_il = 1/√(2 π) f_i | î.lm.

Index î = 1, …, i - 1, i + 1, …, d. Observation ŷ_m equals the mean of a bin y̅_m for the histogram and ŷ_m = y_m for the Parzen window and k nearest neighbour <cit.>. For the histogram, the variance-covariance matrix is given by

    _l = 1/n_l∑_j = 1^v k_lj (ŷ_j - μ_l) (ŷ_j - μ_l)^⊤,

where k_lj denote frequencies and μ_l = (μ_1l, …, μ_dl)^⊤. Number of bins v is replaced by n for the Parzen window or k nearest neighbour. Once the variance-covariance matrix is known, the correlation matrix

    _l = l^-1/2_ll^-1/2 = {ρ_iĩl}

and its inverse

    _l^-1 = {g_iĩl}

are determined. It can be proved that

    σ_iil = g_iilσ_il^2,  i = 1, …, d,

which yields

    l = l^1/2_ll^1/2 = {σ_iĩl}.

The idea is to prevent the component from flowing away from the global mode as at least one component is supposed to be in its vicinity. This yields

    f(y = ŷ_m | _l) = 1/√((2 πε)^dl) = f_lm,

wherefrom

    ε = max{1,  (f_lm√((2 π)^dl) )^-2/d}.

The lower limit of ε is set to 1 in order to prevent f(y = ŷ_m | _l) to be less than the empirical density at the global mode f_lm. This finally yields

    l = εl = {εσ_iĩl}.

The loose restraints are treated exactly the same way as in <cit.>. The loop 1 in Fig. <ref> splits the dataset into two clusters, the one corresponding to the currently observed component and the residue. The former is used to estimate the component weight and rough component parameters. The latter is split into two clusters repeatedly until c ≥ c_max or c ≥ v or c ≥ k. When the total of positive relative deviations attains its minimum, the enhanced component parameters are estimated and the loop 1 stops.



 §.§ Enhanced component parameter estimation

Maximum likelihood is employed to obtain enhanced component parameters. For the histogram, enhanced multivariate normal component
parameters are given by

    μ_l = 1/n_l∑_j = 1^v k_ljŷ_j and l = 1/n_l∑_j = 1^v k_lj (ŷ_j - μ_l) (ŷ_j - μ_l)^⊤.

Index v is replaced by n for the Parzen window or k nearest neighbour.


 §.§ First and second moment calculation

The first and second moment

    m_l = μ_l and V_l = l + μ_lμ_l^⊤ = {V_iĩl}

of the multivariate normal distribution are required for classification of the unassigned observations.


 §.§ Bayes classification of the unassigned observations

When the weight of the unassigned observations w_l≤ D_min (l - 1), then the assignment of new components stops. The constant 0 < D_min≤ 1 is optimized by one of the information criteria. Unassigned observations k_lj are then assumed to belong to the existing classes. The classification of unassigned observations is accomplished by the Bayes decision rule <cit.>

    l =  largmax w_l f(y_j | _l) 
     w_l = w_l + k_lj/n,  m_il = m_il + k_lj (y_ij - m_il)/n w_l and  V_iĩl = V_iĩl + k_lj (y_ij y_ĩj - V_iĩl)/n w_l,

where k_lj is added to the lth class and the component weight and both moments are recalculated <cit.>. Once all v bin means or all n observations are processed, the predictive mixture parameters are gained by inverting Equation (<ref>).

The overall optimal number of components, component weights and component parameters are those resulting in the lowest value of the information criterion.



§ EXAMPLE

For demonstration purpose a multivariate normal dataset with unrestricted variance-covariance matrices is studied. By varying dataset dimension d, dataset size n, number of components c, the seed, range of means μ_l, range of eigenvalues Λ_l of variance-covariance matrix l = P_lΛ_lP_l^⊤ and eigenvectors stored in P_l, a great variety of very different datasets may be generated and the package can be tested.



 §.§ Random dataset generation

This section demonstrates how multivariate normal datasets with unrestricted variance-covariance matrices are generated. Usually the d dimensional datasets of size n are generated for some chosen component weights w_l and component parameters μ_l and l, where l ranges from 1 to c. Here the diversity of datasets is achieved by varying

d <- 2; n <- 50000; c <- 20; set.seed(123)

In order to increase the diversity of datasets, the chosen component weights are replaced by randomly generated ones

w <- runif(c, 0.1, 0.9); w <- w / sum(w)

Other parameters affecting dataset generation are ranges of means mu and eigenvalues lambda

mu <- c(-100, 100); lambda <- c(1, 100)

Component means Mu and diagonal matrices of eigenvalues Lambda are uniformly distributed. The limits of the distributions are in mu and lambda. Component variance-covariance matrices Sigma are generated by the singular-value decomposition of uniformly distributed rectangular matrices. Here the limits of the distribution are -1 and 1.

Mu <- list(); Sigma <- list()
for (l in 1:c) 
  Mu[[l]] <- runif(d, mu[1], mu[2])
  Lambda <- diag(runif(d, lambda[1], lambda[2]), nrow = d, ncol = d)
  P <- svd(matrix(runif(d * d, -1, 1), nc = d))u
  Sigma[[l]] <- P 

Numbers of observations in classes n, Mu and Sigma are strored in a format required by the rebmix package

n <- round(w * n); Theta <- list()
for (l in 1:c) 
  Theta[[paste0("pdf", l)]] <- rep("normal", d)
  Theta[[paste0("theta1.", l)]] <- Mu[[l]]
  Theta[[paste0("theta2.", l)]] <- as.vector(Sigma[[l]])

Now, all arguments entering the RNGMIX method are ready for finite mixture generation

mvnorm <- RNGMIX(model = "RNGMVNORM", Dataset.name = "mvnorm_1",
  rseed = -1, n = n, Theta = Theta)

The method returns an object of class "RNGMVNORM" containing Dataset.name, rseed, n, Theta, Dataset, Zt, w, Variables, ymin and ymax slots. Refer to help("RNGMIX-class") for details. Only three slots are shown if mvnorm is called

> mvnorm
An object of class "RNGMVNORM"
Slot "w":
 [1] 0.0305 0.0676 0.0395 0.0746 0.0788 0.0126 0.0483 0.0753 0.0500 0.0430
[11] 0.0800 0.0428 0.0594 0.0516 0.0169 0.0758 0.0275 0.0124 0.0335 0.0799
Slot "ymin":
[1] -116 -100
Slot "ymax":
[1] 113 119

To plot the dataset, call

plot(mvnorm)

In Fig. <ref> the clustered multivariate normal dataset with considerably overlapped components is given. True cluster membership Zt is known. Each cluster is depicted in a different colour.

If the code of this section is called anew by only changing

set.seed(124); lambda <- c(1, 10)

a very different and less overlapped dataset in Fig. <ref> is obtained.


 §.§ Finite mixture estimation

This section shows how normal finite mixture with unrestricted variance-covariance matrices is estimated given that the dataset from the previous section is known. By calling the REBMIX method, finite mixture for the data frame mvnorm_1 stored in the list mvnorm@Dataset is estimated. From experience it is reasonable to set Preprocessing to "histogram" forn > 1000. Otherwise "Parzen window" or "k-nearest neighbour" may be preferable. The algorithm always searches the optimal number of components between1andc_max. Maximum number of componentsc_maxis, as a rule, chosen to be larger than the maximum expected number of components. For the complete list of arguments refer to help("REBMIX-methods").

mvnormest <- REBMIX(model = "REBMVNORM", Dataset = mvnorm@Dataset,
  Criterion = "BIC", Preprocessing = "histogram", cmax = ceiling(1.2 * c))

Object mvnormest of class "REBMVNORM" contains Dataset, Preprocessing, cmax, Criterion, Variables, pdf, theta1, theta2, K, y0, ymin, ymax, ar, Restraints, w, Theta, summary, pos, opt.c, opt.IC, opt.logL, opt.D, all.K and all.IC slots. Refer to help("REBMIX-class") for details. To produce object summaries, the summary method is called.

> summary(mvnormest)
   Dataset Preprocessing Criterion  c v/k     IC    logL   M
1 mvnorm_1     histogram       BIC 24  46 931191 -464822 143
Maximum logL = -464822 at pos = 1.

In the particular case only one dataset is stored in mvnorm@Dataset. However, it is possible to store more than one dataset to mvnorm@Dataset. All datasets are then processed with the same arguments when the REBMIX method is called. Refer to the example section of help("RNGMIX-methods") for details. To print the coefficients, the coef method is called, where pos stands for the desired row number in mvnormest@summary to be printed or plotted.

> coef(mvnormest, pos = 1)
   comp1  comp2 comp3  comp4  comp5  comp6
w 0.0328 0.0335 0.071 0.0339 0.0712 0.0124
 ⋮
              1      2
theta1.1  -38.6 -17.28
theta1.2  -38.4 -17.43
theta1.3  -24.5  21.91
theta1.4  -24.9  97.55
theta1.5  -68.3 -72.67
theta1.6  -32.6 -20.44
 ⋮
            1-1      1-2      2-1    2-2
theta2.1   27.1  -3.0596  -3.0596   5.38
theta2.2   11.1  -0.9241  -0.9241   3.79
theta2.3   24.0 -15.9239 -15.9239  82.41
theta2.4   11.7  -2.8936  -2.8936  15.97
theta2.5   42.6  13.9595  13.9595  48.01
theta2.6   10.6   3.8294   3.8294  14.21
 ⋮
The plot method is called to plot the mvnormest object in Fig <ref>. Refer to help("plot-methods") for details.

plot(mvnormest, pos = 1, nrow = 3, ncol = 1, what = c("density", "marginal"))


 §.§ Bootstrapping

Bootstrapping is dealt with in this section. It is the practice of estimating standard errors and coefficients of variation by measuring those properties when sampling from approximating distribution stored in the mvnormest object. If Bootstrap = "parametric", the bootstrap datasets are generated by executing the RNGMIX method within the boot method based on the component weights and component parameters of the mvnormest object. The number of bootstrap datasets is controlled by the B argument. If Bootstrap = "nonparametric", the bootstrap datasets are generated directly from the mvnorm_1 data frame by applying the sample.int built-in R method. The boot method

mvnormboot <- boot(mvnormest, pos = 1, Bootstrap = "parametric", B = 10)

returns object mvnormboot of class "REBMVNORM.boot". It contains x, rseed, pos, Bootstrap, B, n, replace, prob, c, c.se, c.cv, c.mode, c.prob, w, w.se, w.cv, Theta, Theta.se and Theta.cv slots. Refer to help("REBMIX.boot-class") for details. Only the most important five slots are shown if mvnormboot is called.

> mvnormboot
An object of class "REBMVNORM.boot"
Slot "c":
 [1] 24 23 21 23 23 23 24 22 24 23
Slot "c.se":
[1] 0.943
Slot "c.cv":
[1] 0.041
Slot "c.mode":
[1] 23
Slot "c.prob":
[1] 0.5

The optimal number of components in slot c is returned for each bootstrap dataset. When standard errors and coefficients of variation are calculated, only those datasets are taken into account for which c equals c.mode. In the particular case only five datasets out of ten are considered. The summary method returns coefficients of variation of component weights and component parameters. The probability of identifying exactly c.mode components is returned, too.

> summary(mvnormboot)
     comp1 comp2 comp3 comp4 comp5 comp6
w.cv  0.46 0.692 0.673 0.679 0.235 0.439
 ⋮
                    1       2
theta1.1.cv   -0.0306  -0.127
theta1.2.cv   -0.1466  -2.040
theta1.3.cv   -0.2214   2.056
theta1.4.cv   -0.2033   2.174
theta1.5.cv   -0.5526   3.060
theta1.6.cv    4.5885  -0.419
 ⋮
               1-1     1-2     2-1   2-2
theta2.1.cv  0.593   0.949   0.949 0.895
theta2.2.cv  0.534  -1.822  -1.822 1.715
theta2.3.cv  0.538  -1.194  -1.194 1.100
theta2.4.cv  0.461  -1.223  -1.223 1.026
theta2.5.cv  0.385   2.452   2.452 0.605
theta2.6.cv  0.495  -1.551  -1.551 0.567
 ⋮
Mode probability = 0.5 at c = 23 components.


 §.§ Clustering

Clustering is used here to group observations in such a way that those with similar features come together and the ones with dissimilar features go apart. It belongs to unsupervised learning. Clustering requires the mvnormest object and pos. It is supposed that components exactly match the clusters and that true cluster membership Zt is unknown. In the particular multivariate normal dataset Zt is known and therefore may enter the RCLRMIX method.

mvnormclu <- RCLRMIX(model = "RCLRMVNORM", mvnormest, pos = 1, Zt = mvnorm@Zt)

Object mvnormclu of class "RCLRMIX" is returned. It contains x, pos, Zt, Zp, c, prob, from, to, EN and ED slots. Refer to help("RCLRMIX-class") for details. Combining mixture components for clustering follows <cit.>. If Zt is known and enters the RCLRMIX method, then the probabilities of correct clustering mvnormclu@prob are calculated, which is convenient to test the package. The optimal number of clusters is given by

copt <- which.max(mvnormclu@prob))
> copt
[1] 16

and does not necessarily coincide with the true number of clusters. To plot the mvnormclu object in Fig <ref> for copt, the plot method is called.

plot(mvnormclu, s = copt)

The probability of correct clustering for copt is then

> mvnormclu@prob[copt]
[1] 0.779

The summary method returns entropy EN and entropy decrease ED<cit.> for all possible numbers of clusters. The information on how the clusters are merged is also printed, e.g., if number of clusters equals23, then cluster12is merged with cluster10.

> summary(mvnormclu)
Number of clusters  1           2          3          4          5
From cluster        4           13         5          10         7
To cluster          1           1          1          5          5
Entropy             -1.86e-13   2.39e+00   7.00e+00   1.28e+01   2.40e+01
Entropy decrease    2.39        4.61       5.79       11.21      12.15
 ⋮
Number of clusters  21         22
From cluster        19         8
To cluster          10         3
Entropy             1.21e+04   1.46e+04
Entropy decrease    2533.21    3150.78


 §.§ Classification

Classification is a process of categorization where observations are recognized, differentiated and understood on the basis of a train subset. It belongs to supervised learning, where for the train subset the true class membership Zt is available. The simulated multivariate normal dataset mvnorm_1 is first combined by columns with Zt from section Random dataset generation.

Dataset <- cbind(mvnorm@Zt, mvnorm@Datasetmvnorm_1)

Next, Dataset is split into train60 %and test40 %datasets.

Mvnorm <- split(p = 0.6, Dataset = Dataset, class = 1)

The split method returns the object Mvnorm of class "RCLS.chunk". Argument class stands for the column number in Dataset containing the class membership information. Refer to help("split-methods") for details. By calling the REBMIX method, finite mixtures for the train subsets in the list Mvnorm@train are estimated. Maximum number of components cmax is set here to5.

mvnormest <- REBMIX(model = "REBMVNORM", Dataset = Mvnorm@train,
  Preprocessing = "histogram", cmax = 5, Criterion = "BIC")

Classification requires the list of objects mvnormest of lengtho. In the particular example number of chunksoequals1. Dataset is a data frame containing test dataset Mvnorm@test. Factor of true class membership Zt for the test dataset enters the RCLSMIX method, too.

mvnormcla <- RCLSMIX(model = "RCLSMVNORM", x = list(mvnormest),
  Dataset = Mvnorm@test, Zt = Mvnorm@Zt)

Object mvnormcla of class "RCLSMIX" is returned. It contains x, o, Dataset, s, ntrain, P, ntest, Zt, Zp, CM, Accuracy, Error, Precision, Sensitivity, Specificity and Chunks slots. The summary method returns table containing confusion matrix CM for multiclass classifier and proportion of all test observations that are classified wrongly Error.

> summary(mvnormcla)
    Test Predictive Frequency
1      1          1       322
2      2          1         0
3      3          1         0
4      4          1         0
5      5          1         0
6      6          1         0
7      7          1         0
8      8          1         0
9      9          1         0
10    10          1         0
 ⋮
391   11         20         0
392   12         20         0
393   13         20         0
394   14         20         0
395   15         20         6
396   16         20         0
397   17         20         0
398   18         20         0
399   19         20         0
400   20         20      1484
Error = 0.0662.

Refer to help("RCLSMIX-class") for details. To plot the mvnormcla object in Fig <ref>, the plot method is called.

plot(mvnormcla)


§ SUMMARY

The rebmix package is aimed to generate, estimate, cluster and classify finite mixture models. Variables can be continuous, discrete or mixed, independent or dependent and may follow normal, lognormal, Weibull, gamma, von Mises, binomial, Poisson or Dirac parametric families. The package has been compared recently with the well established flexmix package by <cit.>. It has turned out that the rebmix package is especially favourable for large datasets in combination with the "histogram" preprocessing due to its numerical stability and computational speed. The preprocessing of observations is beneficial for large datasets. However, for small datasets it results in slightly poorer estimates as compared to the EM based packages. REBMIX can also be used to assess an initial set of unknown parameters and number of components, e.g., for the EM based packages.

In this paper the features for multivariate normal mixture models with unrestricted variance-covariance matrices are presented. For non-normal mixtures readers may address <cit.> where they can also find thorough theoretical backgrounds of the algorithm. New features for finite mixture modeling, bootstrapping, clustering and classification are presented. The example enables the study of a great variety of very different datasets by modifying d, n, c, the seed, mu and lambda. The number specifying the fraction of observations for training p may be changed as well. All this enables further evaluations of the package.



§ ACKNOWLEDGMENT

This paper is a part of research work within program Nr. P2-0182 entitled Development evaluation financed by the Slovenian Ministry of Education, Science and Sport.


Marko Nagode

University of Ljubljana

Faculty of Mechanical Engineering

Aškerčeva 6

1000 Ljubljana

Slovenia
mailto:Marko.Nagode@fs.uni-lj.siMarko.Nagode@fs.uni-lj.si.