Description Usage Arguments Details Value References Author(s) See Also Examples
otrimle
searches for G
approximately Gaussianshaped
clusters with/without noise/outliers. The method's tuning controlling the noise
level is adaptively chosen based on the data.
1 2 
data 
A numeric vector, matrix, or data frame of observations. Rows correspond
to observations and columns correspond to variables. Categorical
variables and 
G 
An integer specifying the number of clusters. 
initial 
An integer vector specifying the initial cluster
assignment with 
logicd 
A vector defining a grid of log(icd) values, where icd
denotes the improper constant density. If 
npr.max 
A number in 
erc 
A number 
beta 
A nonnegative constant. This is the beta penalty coefficient introduced in Coretto and Hennig (2016). 
iter.max 
An integer value specifying the maximum number of iterations allowed in the underlying ECMalgorithm. 
tol 
Stopping criterion for the underlying ECMalgorithm. An ECM iteration
stops if two successive improper loglikelihood values are within

ncores 
an integer value defining the number of cores used for parallel
computing. When 
monitor 
logical. If 
The otrimle
function computes the OTRIMLE solution based on the
ECMalgorithm (expectation conditional maximization algorithm)
proposed in Coretto and Hennig (2017).
The otrimle criterion is minimized over the logicd
grid of
log(icd)
values using parallel computing based on the
foreach
.
Note that, depending on the BLAS/LAPACK setting, increasing ncores
may not
produce the desired reduction in computing time. The latter
happens when optimized linear algebra routines are in use (e.g.
OpenBLAS, Intel Math Kernel Library (MKL), etc.). These optimized shared
libraries already implement multithreading. Therefore, in this case
increasing ncores
may only reduce the computing time marginally.
Occasionally, there may be datasets for which the function does not provide a
solution based on default arguments. This corresponds to
code=0
and flag=1
or flag=2
in the output (see
Valuesection below). This usually happens when some (or all) of the
following circumstances occur: (i) erc
is too
large; (ii) npr.max
is too large; (iii) choice of the initial
partition. Regarding (i) and (ii) it is not possible to give numeric
references because whether these numbers are too large/small
strongly depends on the sample size and the dimensionality of the
data. References given below explain the relationship between
these quantities.
It is suggested to try the following whenever a code=0
nonsolution occurs. Set the logicd
range wide enough
(e.g. logicd=seq(500,5, length=50)
), choose erc=1
,
and a low choice of npr.max
(e.g. npr.max=0.02
).
Monitor the solution with the criterion profiling plot
(plot.otrimle
). According to the criterion profiling
plot change logicd
, and increase erc
and npr.max
up to the point when a "clear" minimum in the criterion profiling plot
is obtained. If this strategy does not work it is suggested to
experiment with a different initial partitions (see initial
above).
TBA: Christian may add something about the beta
here.
The pi
object returned by the rimle
function (see
Value) corresponds to the vector of pi
parameters in
the underlying pseudomodel (1) defined in Coretto and Hennig (2017).
With logicd = Inf
the rimle
function approximates the
MLE for the plain Gaussian mixture model with eigenratio
covariance regularization, in this case the the first element of the
pi
vector is set to zero because the noise component is not
considered. In general, for iid sampling from finite mixture models
context, these pi parameters define expected clusters'
proportions. Because of the noise proportion constraint in the RIMLE,
there are situations where this connection may not happen in
practice. The latter is likely to happen when both logicd
and
npr.max
are large. Therefore, estimated expected clusters' proportions
are reported in the exproportion
object of the rimle
output, and these are computed based on the improper posterior
probabilities given in tau
. See Coretto and Hennig (2017) for
more discussion on this.
An earlier approximate version of the algorithm was originally proposed in Coretto and Hennig (2016). Software for the original version of the algorithm can be found in the supplementary materials of Coretto and Hennig (2016).
An S3 object of class 'otrimle'
providing the optimal (according to
the OTRIMLE criterion) clustering. Output components are as follows:
code 
An integer indicator for the convergence.

flag 
A character string containing one or more flags related to
the EM iteration at the optimal icd.

iter 
Number of iterations performed in the underlying EMalgorithm at the
optimal 
logicd 
Resulting value of the optimal 
iloglik 
Resulting value of the improper likelihood. 
criterion 
Resulting value of the OTRIMLE criterion. 
pi 
Estimated vector of the 
mean 
A matrix of dimension 
cov 
An array of size 
tau 
A matrix of dimension 
smd 
A matrix of dimension 
cluster 
A vector of integers denoting cluster assignments for each
observation. It's 
size 
A vector of integers with sizes (counts) of each cluster. 
exproportion 
A vector of estimated expected clusters' proportions (see Details). 
optimization 
A data.frame with the OTRIMLE optimization profiling. For each
value of 
Coretto, P. and C. Hennig (2016). Robust improper maximum likelihood: tuning, computation, and a comparison with other methods for robust Gaussian clustering. Journal of the American Statistical Association, Vol. 111(516), pp. 16481659. doi: 10.1080/01621459.2015.1100996
P. Coretto and C. Hennig (2017). Consistency, breakdown robustness, and algorithms for robust improper maximum likelihood clustering. Journal of Machine Learning Research, Vol. 18(142), pp. 139. https://jmlr.org/papers/v18/16382.html
Pietro Coretto pcoretto@unisa.it https://pietrocoretto.github.io
plot.otrimle
,
InitClust
,
rimle
,
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75  ## Load Swiss banknotes data
data(banknote)
x < banknote[,1]
## Perform otrimle clustering with default arguments
set.seed(1)
a < otrimle(data=x, G=2, logicd=c(Inf, 50, 10), ncores=1)
## Plot clustering
plot(a, data=x, what="clustering")
## Plot OTRIMLE criterion profiling
plot(a, what="criterion")
## Plot Improper loglikelihood profiling
plot(a, what="iloglik")
## PP plot of the clusterwise empirical weighted squared Mahalanobis
## distances against the target distribution pchisq(, df=ncol(data))
plot(a, what="fit")
plot(a, what="fit", cluster=1)
## Perform the same otrimle as before with nonzero penalty
set.seed(1)
b < otrimle(data=x, G=2, beta = 0.5, logicd=c(Inf, 50, 10), ncores=1)
## Plot clustering
plot(b, data=x, what="clustering")
## Plot OTRIMLE criterion profiling
plot(b, what="criterion")
## Plot Improper loglikelihood profiling
plot(b, what="iloglik")
## PP plot of the clusterwise empirical weighted squared Mahalanobis
## distances against the target distribution pchisq(, df=ncol(data))
plot(b, what="fit")
plot(b, what="fit", cluster=1)
## Not run:
## Perform the same example using the finer default grid of logicd
## values using multiple cores
##
a < otrimle(data = x, G = 2)
## Inspect the otrimle criterionvslogicd
plot(a, what = 'criterion')
## The minimum occurs at a$logicd=9, and exploring a$optimization it
## cane be seen that the interval [12.5, 4] brackets the optimal
## solution. We search with a finer grid located around the minimum
##
b < otrimle(data = x, G = 2, logicd = seq(12.5, 4, length.out = 25))
## Inspect the otrimle criterionvslogicd
plot(b, what = 'criterion')
## Check the difference between the two clusterings
table(A = a$cluster, B = b$cluster)
## Check differences in estimated parameters
##
colSums(abs(a$mean  b$mean)) ## L1 distance for mean vectors
apply({a$covb$cov}, 3, norm, type = "F") ## Frobenius distance for covariances
c(Noise=abs(a$nprb$npr), abs(a$cprb$cpr)) ## Absolute difference for proportions
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.