MIXTURES GAUSSIANES I ALGORISME EM

Andrea Baena Espejo, 2022-01-24

library(rebmix)
library(stats)
library(ggplot2)
library(tidyr)
library(gridExtra)
library(cluster)
library(jpeg)
library(cluster.datasets)

1. Implementació de l'algorisme K-mitjanes

data(all.mammals.milk.1956)

Dins de la llibreria "cluster.datasets" de R, podem trobar diverses bases de dades que es poden emprar per fer anàlisi de clústers. En el nostre cas, hem triat les dades "all.mammals.milk.1956". Aquesta base de dades conté 25 observacions sobre les 6 variables següents:

• name: un vector de caràcters per al nom dels animals

• water: un vector numèric per al contingut d'aigua a la mostra de llet

• protein: un vector numèric per a la quantitat de proteïna a la mostra de llet

• fat: un vector numèric per al contingut de greix a la mostra de llet

• lactose: un vector numèric per a la quantitat de lactosa a la mostra de llet

• ash: un vector numèric per a la quantitat de minerals a la mostra de llet

#Veiem com són les dades
str(all.mammals.milk.1956)

'data.frame':	25 obs. of  6 variables:
 $ name   : chr  "Horse" "Orangutan" "Monkey" "Donkey" ...
 $ water  : num  90.1 88.5 88.4 90.3 90.4 87.7 86.9 82.1 81.9 81.6 ...
 $ protein: num  2.6 1.4 2.2 1.7 0.6 3.5 4.8 5.9 7.4 10.1 ...
 $ fat    : num  1 3.5 2.7 1.4 4.5 3.4 1.7 7.9 7.2 6.3 ...
 $ lactose: num  6.9 6 6.4 6.2 4.4 4.8 5.7 4.7 2.7 4.4 ...
 $ ash    : num  0.35 0.24 0.18 0.4 0.1 0.71 0.9 0.78 0.85 0.75 ...

#Mostrem les 6 primeres observacions
head(all.mammals.milk.1956)

A data.frame: 6 × 6

	name	water	protein	fat	lactose	ash
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	Horse	90.1	2.6	1.0	6.9	0.35
2	Orangutan	88.5	1.4	3.5	6.0	0.24
3	Monkey	88.4	2.2	2.7	6.4	0.18
4	Donkey	90.3	1.7	1.4	6.2	0.40
5	Hippo	90.4	0.6	4.5	4.4	0.10
6	Camel	87.7	3.5	3.4	4.8	0.71

Ara que hem vist com són les nostres dades, volem classificar aquestes espècies en grups, segons com és de semblant la composició de la seva llet respecte a les variables observades.

#Començem fem una representació dels valors que pren cada 
plot1 <- all.mammals.milk.1956 %>% 
    ggplot(aes(x = "all mammals", y = water)) + 
    geom_jitter(width = .025, height = 0, size = 2, alpha = .5, color = "blue") +
    labs(x = "", y="aigua")
  
plot2 <-  all.mammals.milk.1956 %>%
  ggplot(aes(x = "all mammals", y = protein)) + 
    geom_jitter(width = .02, height = 0, size = 2, alpha = .6,  color = "orange") +
    labs(x = "", y="proteïna")
  
plot3 <-  all.mammals.milk.1956 %>%
    ggplot(aes(x = "all mammals", y = fat)) + 
    geom_jitter(width = .02, height = 0, size = 2, alpha = .6,  color = "green") +
    labs(x = "", y="greix")
  
plot4 <-  all.mammals.milk.1956 %>%
    ggplot(aes(x = "all mammals", y = lactose)) + 
    geom_jitter(width = .02, height = 0, size = 2, alpha = .6,  color = "red") +
    labs(x = "", y="lactosa")
  
plot5 <-  all.mammals.milk.1956 %>%
    ggplot(aes(x = "all mammals", y = ash)) + 
    geom_jitter(width = .02, height = 0, size = 2, alpha = .6,  color = "violet") +
    labs(x = "", y="minerals")
  
grid.arrange(plot1, plot2, plot3, plot4, plot5)

Cada variable té un comportament diferent i podríem identificar grups de mamífers en cadascuna individualment, però el propòsit és fer una agrupació que tingui en compte totes les variables alhora. El primer que ens plantegem és: "quin és el valor K de clústers hem d'usar per aplicar l'algorisme K-mitjanes?"

Per estimar quin és el millor nombre de grups, apliquem l'algorisme per a K=$\left\lbrace{1,\ldots, 20}\right\rbrace$ i comparem el valor de la suma wss (within sum of squares (wss)) en cada cas. Aquest valor dona una mesura de l'error comès en l'agrupació. Per aplicar l'algorisme en aquest rang de K, usem la funció kmeans implementada en el paquet 'stats' de R.

obs <- all.mammals.milk.1956[,2:6]
wss <- (nrow(obs)-1)*sum(apply(obs,2,var))
for(i in 1:20){
    wss[i] <- sum(kmeans(obs, centers=i)$withinss)
    }
    plot(1:20, wss, type="b", xlab="No. of Clusters", ylab="wss")
round(wss,3)

1. 7004.36
2. 2186.735
3. 1388.795
4. 342.099
5. 212.646
6. 181.055
7. 109.657
8. 99.472
9. 92.687
10. 88.201
11. 86.034
12. 41.961
13. 62.676
14. 53.091
15. 59.361
16. 17.127
17. 16.037
18. 42.265
19. 12.614
20. 7.621

Si analitzem la gràfica, veiem un punt d'inflexió (colze) en K = 4. Per tant, agrupem les nostres dades en 4 grups, aplicant kmeans

set.seed(20)
allmamals.clu <- kmeans(obs, 4, nstart = 20)
allmamals.clu

K-means clustering with 4 clusters of sizes 7, 6, 10, 2

Cluster means:
             water           protein              fat          lactose
1 81.1857142857143  7.42857142857143  6.9000000000000 4.01428571428571
2 68.3333333333333  9.55000000000000 17.4166666666667 2.91666666666667
3 88.5000000000000  2.57000000000000  2.8000000000000 5.68000000000000
4 45.6500000000000 10.15000000000000 38.4500000000000 0.45000000000000
                ash
1 0.931428571428571
2 1.330000000000000
3 0.485000000000000
4 0.690000000000000

Clustering vector:
 [1] 3 3 3 3 3 3 3 1 1 1 1 3 3 1 3 1 1 2 2 2 2 2 2 4 4

Within cluster sum of squares by cluster:
[1]  63.5349142857143 191.9610000000000  59.4122500000000  27.1912000000000
 (between_SS / total_SS =  95.1 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"      

Veiem les mitjanes dels clústers, és a dir, els prototips (cluster means) i el vector d'assignacions (clustering vector).

#Nombre d'observacions assignades a cada clúster
table(allmamals.clu$cluster)

 1  2  3  4 
 7  6 10  2 

Per últim, validem com ha estat de bona l'agrupació. Per fer-ho, usem el mètode de la silueta, a través de la funció silhouette del paquet 'cluster' de R.

sil <- silhouette(allmamals.clu$cluster, dist(obs))
plot(sil)

El gràfic de la silueta anterior ens mostra que la nostra agrupació amb 4 grups és bona perquè no hi ha una amplada negativa de la silueta i la majoria dels valors són més grans que 0,5. La mitjana val 0,6.

A continuació, implementem l'algorisme K-mitjanes usant una funció pròpia

kmeanspropi <-function(data, K=4, tol=1e-5){
    # Inicialitzem els prototips (aleatòriament d'entre el conjunt dobservacions)
    prototips <-data[sample.int(nrow(data),K), 1:2]
    
    # La variable dif emmagatzema la diferència entre els prototips; la inicialitzem  a un valor arbitràriament gran (per entrar en el 'while')
    dif <-1000
    
    
    # Inicialitzem un vector per guardar les assignacions dels punts als clústers a cada iteració
    cluster <-rep(0,nrow(data))
    
   
    it<-1
    convergencia<-FALSE
    while(dif>=tol && convergencia==FALSE) {
        it=it+1
        if(dif<=tol){
            convergencia<-TRUE
        }
        prototips_old <-prototips
        
        # Assignem cada punt del conjunt de dades al clúster més proper (distància euclidiana)
        for (i in 1:nrow(data)){
            dist_min<-10e10
            for (prototip in 1:nrow(prototips)){
                dist_prototip <-sum((prototips[prototip,1:2]-data[i,1:2])^2)
                if (dist_prototip<=dist_min){
                    cluster[i] <-prototip
                    dist_min <-dist_prototip
                }
            }
        } 

        
        # Reestimem els nous prototips (mitjana dels punts que han estat assignats al clúster corresponent en el pas previ)
        for (alpha in 1:nrow(prototips)){
            prototips[alpha,1:2]=apply(as.matrix(data[cluster==alpha,1:2]),2,mean)
        }
        
        dif <-mean(abs(as.matrix(prototips_old-prototips)))
    }
    
    return(list(result=data.frame(data,cluster)))
}

Emprem la funció que acabem d'implementar per agrupar les dades anteriors (all.mammals.milk.1956).

allmamals.clu2 <-kmeanspropi(obs,4,1e-5)

Per últim, validem els resultats amb el mètode de la silueta.

sil <- silhouette(allmamals.clu2$result$cluster, dist(obs))
plot(sil)

allmamals.clu2$resul

A data.frame: 25 × 6

water	protein	fat	lactose	ash	cluster
<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
90.1	2.6	1.0	6.9	0.35	4
88.5	1.4	3.5	6.0	0.24	4
88.4	2.2	2.7	6.4	0.18	4
90.3	1.7	1.4	6.2	0.40	4
90.4	0.6	4.5	4.4	0.10	4
87.7	3.5	3.4	4.8	0.71	4
86.9	4.8	1.7	5.7	0.90	4
82.1	5.9	7.9	4.7	0.78	1
81.9	7.4	7.2	2.7	0.85	1
81.6	10.1	6.3	4.4	0.75	1
81.6	6.6	5.9	4.9	0.93	1
86.5	3.9	3.2	5.6	0.80	4
90.0	2.0	1.8	5.5	0.47	4
82.8	7.1	5.1	3.7	1.10	1
86.2	3.0	4.8	5.3	0.70	4
82.0	5.6	6.4	4.7	0.91	1
76.3	9.3	9.5	3.0	1.20	1
70.7	3.6	17.6	5.6	0.63	3
71.3	12.3	13.1	1.9	2.30	3
72.5	9.2	12.6	3.3	1.40	3
65.9	10.4	19.7	2.6	1.40	3
64.8	10.7	20.3	2.5	1.40	3
64.8	11.1	21.2	1.6	0.85	3
46.4	9.7	42.0	0.0	0.85	2
44.9	10.6	34.9	0.9	0.53	2

Observem que amb les dues funcions, les dades s'han agrupat en els mateixos grups

Clustering vector (kmeans) --> 3 3 3 3 3 3 3 1 1 1 1 3 3 1 3 1 1 2 2 2 2 2 2 4 4

Cluster (funció propia) --> 4 4 4 4 4 4 4 1 1 1 1 4 4 1 4 1 1 3 3 3 3 3 3 2 2

Aplicació de l'algorisme K-mitjanes en compressió d'imatges

Vegem com podem usar l'algorisme d'agrupació K-mitjanes per comprimir imatges.

#Llegim la imatge
img <- readJPEG("imagecompression.jpg")

#Mostrem la imatge
options(repr.plot.width=5,repr.plot.height=5)
plot(c(0,1),c(0,1),"n",xlab="",ylab="")
rasterImage(img,0,0,1,1)

#Extraiem la dimensió de la imatge
imgDm <- dim(img)
imgDm

1. 410
2. 710
3. 3

#Separem la imatge en les seves components RGB (colors primaris)
imgRGB <- data.frame(
  x = rep(1:imgDm[2], each = imgDm[1]),
  y = rep(imgDm[1]:1, imgDm[2]),
  R = as.vector(img[,,1]),
  G = as.vector(img[,,2]),
  B = as.vector(img[,,3])
)

#Creem una funció per personalitzar la representació de les imatges
plotTheme <- function() {
  theme(
    panel.background = element_rect(
      size = 3,
      colour = "black",
      fill = "white"),
    axis.ticks = element_line(
      size = 2),
    panel.grid.major = element_line(
      colour = "gray80",
      linetype = "dotted"),
    panel.grid.minor = element_line(
      colour = "gray90",
      linetype = "dashed"),
    axis.title.x = element_text(
      size = rel(1.2),
      face = "bold"),
    axis.title.y = element_text(
      size = rel(1.2),
      face = "bold"),
    plot.title = element_text(
      size = 20,
      face = "bold",
      vjust = 1.5)
  )
}

#Usem la funicó que hem definit per mostrar primer la nostra imatge original
ggplot(data = imgRGB, aes(x = x, y = y)) + 
    geom_point(colour = rgb(imgRGB[c("R", "G", "B")])) +
    labs(title = paste("Imatge original")) +
    xlab("x") +
    ylab("y") +
plotTheme()

Apliquem l'algorisme k-mitjanes amb K = 10 clústers sobre aquesta imatge. D'aquesta manera, agrupem tots els colors existents a la imatge en tan sols 10 colors, i visualitzem el resultat obtingut.

img10 <- kmeans(imgRGB[, c("R", "G", "B")], centers = 10)
colors <- rgb(img10$centers[img10$cluster,])

ggplot(data = imgRGB, aes(x = x, y = y)) + 
  geom_point(colour = colors) +
  labs(title = paste("k-mitjanes amb k = 10 colors")) +
  xlab("x") +
  ylab("y") + 
plotTheme()

2. Algorisme EM

Paquet REBMIX

El paquet rebmix proporciona funcions R per a generar i estimar models de mixtures aleatòries univariants i multivariants finites. Les variables poden ser contínues, discretes, independents o dependents i poden seguir una distribució normal, lognormal, Weibull, gamma, binomial, Poisson o von Mises. La versió 2.11.0 introdueix l'algorisme EM per a l'estimació millorada dels paràmetres del model de mixtures gaussianes

Generació d'un model de mixtura de gaussianes

# Definim un model de mixtura de gaussianes

n <- c(50, 150, 75, 125, 100) 

# La mixtura té cinc components. El total d'observacions és 500. La primera component en té 50, per tant π_1 = 0.1
# Anàlogament π_2 = 0.3 , π_3 = 0.15 , π_4 = 0.25 , π_5 = 0.2

Theta <- new("RNGMVNORM.Theta", c = 5, d = 2) 

#Theta contindrà els paràmetres de la mixtura de K = 5 (components) i  d = 2 (dimensió/variables observades)

# mu (mijtanes): vector d-dimensional
a.theta1(Theta, 1) <- c(2.7, 3.7)
a.theta1(Theta, 2) <- c(5.7, 9.1)
a.theta1(Theta, 3) <- c(2, 9)
a.theta1(Theta, 4) <- c(9.5, 6.6)
a.theta1(Theta, 5) <- c(6.3, 0.6)

# Sigma (variàncies): matriu real, simètrica i definida positiva, de dimensió d x d
a.theta2(Theta, 1) <- c(0.9, -0.1, -0.1, 0.4)
a.theta2(Theta, 2) <- c(2.8, -1.3, -1.3, 1.5)
a.theta2(Theta, 3) <- c(0.1, 0, 0, 0.3)
a.theta2(Theta, 4) <- c(1.3, -0.4, -0.4, 0.4)
a.theta2(Theta, 5) <- c(0.5, 0.3, 0.3, 2.5)

# Generem aleatòriament la mostra
mvnorm.simulated <- RNGMIX(model = "RNGMVNORM", Dataset.name = "mvnormdataset", rseed = -1, n = n, Theta = a.Theta(Theta))

RNGMIX Version 2.13.1

Dataset = mvnormdataset

#Vegem les primeres 6 observacions que ens ha generat aquest procés
head(mvnorm.simulated@Dataset$mvnormdataset)

A data.frame: 6 × 2

	1	2
	<dbl>	<dbl>
1	2.53656091344232	3.19628986315287
2	4.23254080375765	3.64640591888040
3	3.31962744405870	3.52091452476089
4	2.88418285207739	3.33881754963756
5	3.84251901887412	4.15033449722182
6	1.41400300222414	4.79735719197654

#Fem un gràfic dels punts (observacions) de la mostra aleatòria que hem generat
plot(mvnorm.simulated)

Estimació usant l'algorisme REBMIX

Usem la funció REMBIX del paquet 'rebmix' per donar una primera estimació màxim versemblant del model. La funció REMBIX es basa en el processament d'histogrames per determinar la funció de densitat.

#Usem la funicó REBMIX del paquet 'rebmix' per estimar el model de mixtura
mvnormest <- REBMIX(model = "REBMVNORM", Dataset = a.Dataset(mvnorm.simulated), Preprocessing = "histogram", cmax = 20, Criterion = "BIC")

REBMIX Version 2.13.1

Dataset = mvnormdataset

#Grafiquem el model de mixtura obtingut
plot(mvnormest)

mvnormest

An object of class "REBMVNORM"
Slot "w":
[1] 0.146895813464989 0.260051424291348 0.194002356457407 0.190199088554626
[5] 0.106088223921686 0.102763093309944
Slot "Theta":
$pdf1
[1] normal normal

$theta1.1
[1] 1.97479881153781 9.11311226419078

$theta2.1
[1] 0.1815090816359257 0.0171283190432554 0.0171283190432554 0.6059462229556800

$pdf2
[1] normal normal

$theta1.2
[1] 9.48473870286257 6.54916960016614

$theta2.2
[1]  2.107541502395406 -0.560617301737565 -0.560617301737565  0.470306432272665

$pdf3
[1] normal normal

$theta1.3
[1] 6.330867466978292 0.254603835905561

$theta2.3
[1] 0.555546717935158 0.461260108702528 0.461260108702528 2.196408687256392

$pdf4
[1] normal normal

$theta1.4
[1] 4.77224898612766 9.62350400277331

$theta2.4
[1]  1.919141324311902 -0.456668653190171 -0.456668653190171  0.852438489761127

$pdf5
[1] normal normal

$theta1.5
[1] 2.82166870872033 3.80276539727387

$theta2.5
[1] 1.375159426399008 0.102756529295187 0.102756529295187 0.586736883527212

$pdf6
[1] normal normal

$theta1.6
[1] 7.56192190706474 8.08442812866403

$theta2.6
[1]  1.798500367510549 -0.189029098392155 -0.189029098392155  0.386302896755296

Slot "summary":
        Dataset Preprocessing Criterion c v/k               IC
1 mvnormdataset     histogram       BIC 6  14 4251.10470170412
               logL  M
1 -2016.79670912967 35

Observem que en aquesta primera estimació, l'algorisme REBMIX ha donat com a resultat una mixtura de 6 components (quan, la mostra està extreta d'una mixtura de 5 components). En aquest sentit, ja podem intuir que aquest mètode té algunes carències. A continuació veurem com millorar aquesta estimació, usant conjuntament l'algorisme REBMIX seguit de l'algorisme EM.

Estratègies d'estimació REBMIX&EM

S'utilitza l'algorisme REBMIX per avaluar els paràmetres inicials de l'algorisme EM. Com l'algorisme REBMIX estima una àmplia gamma de paràmetres pel model de mixtura gaussiana per a diferents nombres de components, s'implementen tres estratègies diferents per aplicar REBMIX&EM:

• Exhaustive: s'usa per executar l'algorisme EM en cada solució dels paràmetres del model de mixtura gaussiana proporcionada per l'algorisme REBMIX.

• Best: s'utilitza un esquema de votació per als paràmetres estimats de l'algorisme REBMIX i s'executa l'algorisme EM només en paràmetres òptims seleccionats. Els millors candidats són triats en funció del valor de la funció de versemblança (la més alta) per a cada nombre de components c d'un mínim especificat cmin a un màxim especificat cmax.

• Single: és útil quan el nombre de bins (barres) en el preprocessament d'histogrames es proporciona com a entrada per a l'algorisme REBMIX. En cas contrari, aquesta estratègia és el mateix que l'estratègia exhaustiva.

SINGLE

Usem primer l'estratègia 'single' de REBMIX&EM. Amb la funció optbins estimem el valor òptim de bins (barres) per a usar en la funció REBMIX.

EM <- new("EM.Control", strategy = "single", variant = "EM", acceleration = "fixed", acceleration.multiplier = 1, tolerance = 1e-04, maximum.iterations = 1000, K = 0)
K <- optbins(Dataset = a.Dataset(mvnorm.simulated), Rule = "Knuth equal", kmin = 2, kmax = 100)
mvnormest.em <- REBMIX(model = "REBMVNORM", Dataset = a.Dataset(mvnorm.simulated), Preprocessing = "histogram", cmax = 20, K = K, Criterion = "BIC", EMcontrol = EM)

REBMIX Version 2.13.1

Dataset = mvnormdataset

summary(mvnormest.em)

        Dataset Preprocessing Criterion c v/k               IC
1 mvnormdataset     histogram       BIC 5  11 4187.41470947874
               logL  M
1 -2003.59553731225 29
Maximum logL = -2003.59553731225 at pos = 1.

plot(mvnormest.em)

mvnormest.em

An object of class "REBMVNORM"
Slot "w":
[1] 0.152279635754161 0.275358384808801 0.199843621310076 0.100173516773880
[5] 0.272344841353083
Slot "Theta":
$pdf1
[1] normal normal

$theta1.1
[1] 2.02428656449457 9.07574734127331

$theta2.1
[1] 0.10981483707693809 0.00352370387753491 0.00352370387753491
[4] 0.46913200788209303

$pdf2
[1] normal normal

$theta1.2
[1] 9.44659813907218 6.63688235904895

$theta2.2
[1]  1.933553716924786 -0.623478526344136 -0.623478526344136  0.450227158988646

$pdf3
[1] normal normal

$theta1.3
[1] 6.269368825668520 0.324380333584456

$theta2.3
[1] 0.467421529302982 0.384510592759515 0.384510592759515 2.219716132871343

$pdf4
[1] normal normal

$theta1.4
[1] 2.68807681282499 3.83486159862061

$theta2.4
[1]  0.9110468749452292 -0.0287615211342961 -0.0287615211342961
[4]  0.4500137858474870

$pdf5
[1] normal normal

$theta1.5
[1] 5.68847779853623 9.13463414245424

$theta2.5
[1]  3.08963942399193 -1.16287620792059 -1.16287620792059  1.24108251562985

Slot "summary":
        Dataset Preprocessing Criterion c v/k               IC
1 mvnormdataset     histogram       BIC 5  11 4187.41470947874
               logL  M
1 -2003.59553731225 29

cat("Iteracions: ", a.summary.EM(mvnormest.em, pos = 1, col.name = "total.iterations.nbr") )

Iteracions:  207

En aquesta estimació sí tenim un model de mixtura de gaussianes amb 5 components. Mostrem una representació gràfica del model i en donem les estimacions dels coeficients de la mixtura (w) i de les mitjanes i les covariàncies (theta). Amb aquesta estimació, la versemblaça val aproximadament -2003.6 Les iteracions necessàries han estat 207.

EXHAUSTIVE

Vegem ara la diferència en els resultats si apliquem l'estratègia exhaustiva.

EM.normal <- new("EM.Control", strategy = "exhaustive", variant = "EM", acceleration = "fixed", acceleration.multiplier = 1, tolerance = 1e-04, maximum.iterations = 1000, K = 0)
mvnormestest.em.normal <- REBMIX(model = "REBMVNORM", Dataset = a.Dataset(mvnorm.simulated), Preprocessing = "histogram", cmax = 15, Criterion = "BIC", EMcontrol = EM.normal)

REBMIX Version 2.13.1

Dataset = mvnormdataset

summary(mvnormestest.em.normal)

        Dataset Preprocessing Criterion c v/k              IC              logL
1 mvnormdataset     histogram       BIC 5  10 4187.3926805214 -2003.58452283358
   M
1 29
Maximum logL = -2003.58452283358 at pos = 1.

plot(mvnormestest.em.normal)

mvnormestest.em.normal

An object of class "REBMVNORM"
Slot "w":
[1] 0.152284985073041 0.274921060140380 0.199843614718997 0.100173515328916
[5] 0.272776824738665
Slot "Theta":
$pdf1
[1] normal normal

$theta1.1
[1] 2.02432486842025 9.07569168839743

$theta2.1
[1] 0.10982941869482202 0.00354384455037328 0.00354384455037328
[4] 0.46900870887927582

$pdf2
[1] normal normal

$theta1.2
[1] 9.44869793968046 6.63550861870828

$theta2.2
[1]  1.931305322008650 -0.622164115054210 -0.622164115054210  0.449014605989509

$pdf3
[1] normal normal

$theta1.3
[1] 6.269368819011667 0.324380140500414

$theta2.3
[1] 0.467421578809030 0.384510553145366 0.384510553145366 2.219715054978578

$pdf4
[1] normal normal

$theta1.4
[1] 2.68807677888187 3.83486155882545

$theta2.4
[1]  0.9110468106082343 -0.0287616236059882 -0.0287616236059882
[4]  0.4500136951479720

$pdf5
[1] normal normal

$theta1.5
[1] 5.69243712082065 9.13204634426902

$theta2.5
[1]  3.09671203987169 -1.16803398563731 -1.16803398563731  1.24419425992394

Slot "summary":
        Dataset Preprocessing Criterion c v/k              IC              logL
1 mvnormdataset     histogram       BIC 5  10 4187.3926805214 -2003.58452283358
   M
1 29

cat("Iteracions: ", a.summary.EM(mvnormestest.em.normal, pos = 1, col.name = "total.iterations.nbr"))

Iteracions:  775

En aquesta estimació hem obtingut pràcticament el mateix valor de versemblança que en el cas de l'estratègia 'single', però l'algorisme ha necessitat 775 iteracions per convergir.

Acceleració de EM

Per abordar la lenta convergència lineal de l'algorisme EM, s'implementen mètodes d'acceleració senzills. Es pot escriure l'increment de l'algorisme EM en cada iteració com $$ \Delta\Theta=\Theta^{i+1}-\Theta^{i} $$

Per reduir el nombre d'iteracions necessàries per a la convergència de l'algorisme EM, aquest increment es pot multiplicar amb algun accelerador. Per tant, l'actualització de cada iteració de EM es converteix en

$$ \Theta^{i+1}=\Theta^{i}+a_{EM}\Delta\Theta $$

El rang desitjable per al multiplicador $a_{EM}$ es troba entre 1.0 i 2.0, on 1.0 dona un increment EM estàndard i 2,0 duplica l'increment EM. Tanmateix, això no vol dir necessàriament que aquesta multiplicació amb un valor de 2,0 duplicarà la velocitat de l'algorisme EM, reduint el nombre necessari d'iteracions per 2. Aquí, 1.5 és un valor segur que accelera majoritàriament l'algorisme, tot conservant bons resultats pels paràmetres estimats. S'implementen els següents mètodes per definir el paràmetre d'acceleració:

• Fixed: el multiplicador $a_{EM}$ s'especifica mitjançant el paràmetre acceleration.multiplier i per a cada iteració de l'algorisme EM, l'increment es multiplica pel valor especificat.

• Line i golden --> es realitza una cerca del valor òptim de $a_{EM}$, per al qual s'aconsegueix el major l'augment de la funció de probabilitat a cada iteració EM. Aquests dos casos no els estudiarem.

#Usem paràmetre d'acceleració fixa, a_{EM}=1.5

EM.fixed1.5 <- new("EM.Control", strategy = "exhaustive", variant = "EM", acceleration = "fixed", acceleration.multiplier = 1.5, tolerance = 1e-04, maximum.iterations = 1000, K = 0)
mvnormest.em.fixed1.5 <- REBMIX(model = "REBMVNORM", Dataset = a.Dataset(mvnorm.simulated), Preprocessing = "histogram", cmax = 15, Criterion = "BIC", EMcontrol = EM.fixed1.5)

REBMIX Version 2.13.1

Dataset = mvnormdataset

summary(mvnormest.em.fixed1.5)

        Dataset Preprocessing Criterion c v/k               IC
1 mvnormdataset     histogram       BIC 5  11 4187.27258636969
               logL  M
1 -2003.52447575772 29
Maximum logL = -2003.52447575772 at pos = 1.

plot(mvnormest.em.fixed1.5)

mvnormest.em.fixed1.5

An object of class "REBMVNORM"
Slot "w":
[1] 0.152168395036592 0.272273714652508 0.199843564849612 0.100173511577744
[5] 0.275540813883544
Slot "Theta":
$pdf1
[1] normal normal

$theta1.1
[1] 2.02444840669804 9.07418964326810

$theta2.1
[1] 0.1098054366810460 0.0037007009177143 0.0037007009177143 0.4667908711983257

$pdf2
[1] normal normal

$theta1.2
[1] 9.46067428686192 6.62718684379548

$theta2.2
[1]  1.919741767875493 -0.614775405528496 -0.614775405528496  0.441646965635795

$pdf3
[1] normal normal

$theta1.3
[1] 6.269368697301642 0.324378681442089

$theta2.3
[1] 0.467421885396044 0.384509915585586 0.384509915585586 2.219706755536705

$pdf4
[1] normal normal

$theta1.4
[1] 2.68807669436544 3.83486141239794

$theta2.4
[1]  0.911046716178501 -0.028762019849457 -0.028762019849457  0.450013406051403

$pdf5
[1] normal normal

$theta1.5
[1] 5.71508350884029 9.11708067669072

$theta2.5
[1]  3.14792563371995 -1.20188916943485 -1.20188916943485  1.26318927352261

Slot "summary":
        Dataset Preprocessing Criterion c v/k               IC
1 mvnormdataset     histogram       BIC 5  11 4187.27258636969
               logL  M
1 -2003.52447575772 29

cat("Iteracions: ", a.summary.EM(mvnormest.em.fixed1.5, pos = 1, col.name = "total.iterations.nbr") )

Iteracions:  560

Veiem com usant el paràmetre fix d'acceleració de 1.5, el nombre d'iteracions s'ha reduït a 560 (quan, amb l'increment estàndard, eren necessàries 775 iteracions). A més, hem obtingut un valor lleugerament més gran de la versemblança: -2003.52