64b28e3e |
1 | ### Regression matrices |
4c9cc558 |
2 | model = Res |
64b28e3e |
3 | K = dim(model$phi)[3] |
4 | valMax = max(abs(model$phi)) |
5 | |
6 | require(fields) |
4c9cc558 |
7 | |
64b28e3e |
8 | if (K<4){ |
9 | par(mfrow = c(1,K)) |
4c9cc558 |
10 | } else op = par(mfrow = c(2, (K+1)/2)) |
11 | |
12 | ## Phi |
64b28e3e |
13 | |
14 | for (r in 1:K){ |
4c9cc558 |
15 | image.plot(t(abs(model$phi[,,r])), |
64b28e3e |
16 | col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66)) |
17 | } |
4c9cc558 |
18 | par(mfrow = c(1,K),oma = c(0,0,3,0)) |
19 | mtext("Regression matrices in each cluster", side=3, line=4, font=2, cex=2, col='red') |
20 | |
21 | par(mfrow = c(1,2), oma=c(0,0,3,0)) |
22 | for (i in 1:4) |
23 | plot(runif(20), runif(20), |
24 | main=paste("random plot (",i,")",sep='')) |
25 | par(op) |
26 | mtext("Four plots", |
27 | side=3, line=4, font=2, cex=2, col='red') |
64b28e3e |
28 | |
29 | ### Zoom onto two classes we want to compare |
30 | kSel = c(1,2) |
31 | par(mfrow = c(1,3)) |
32 | |
33 | for (r in kSel){ |
93285a2d |
34 | image.plot(t(abs(model$phi[,,r])),xaxt="n",yaxt="n", |
64b28e3e |
35 | col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66)) |
36 | } |
37 | image.plot(t(abs(model$phi[,,kSel[1]]-model$phi[,,kSel[2]])), |
38 | col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66)) |
39 | |
40 | ### Covariance matrices |
41 | par(mfrow = c(K, 1)) |
42 | for (r in 1:K){ |
43 | image.plot(matrix(diag(model$rho[,,r]), ncol= 1), |
44 | col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66)) |
45 | } |
46 | |
47 | ### proportions |
48 | Gam = matrix(0, ncol = K, nrow = n) |
49 | gam = Gam |
50 | for (i in 1:n){ |
4c9cc558 |
51 | for (r in 1:K){ |
64b28e3e |
52 | sqNorm2 = sum( (Y[i,]%*%model$rho[,,r]-X[i,]%*%model$phi[,,r])^2 ) |
53 | Gam[i,r] = model$pi[r] * exp(-0.5*sqNorm2)* det(model$rho[,,r]) |
54 | } |
55 | gam[i,] = Gam[i,] / sum(Gam[i,]) |
56 | } |
57 | affec = apply(gam, 1,which.max) |
58 | gam2 = matrix(NA, ncol = K, nrow = n) |
59 | for (i in 1:n){ |
60 | gam2[i, affec[i]] = gam[i, affec[i]] |
61 | } |
62 | boxplot(gam2) |
63 | |
64 | ### Mean in each cluster |
65 | XY = cbind(X,Y) |
66 | XY_class= list() |
67 | meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2]) |
68 | for (r in 1:K){ |
69 | XY_class[[r]] = XY[affec == r, ] |
70 | meanPerClass[,r] = apply(XY_class[[r]], 2, mean) |
71 | } |
72 | |
c0846e18 |
73 | matplot(meanPerClass, type='l') |