remove extra prints
[valse.git] / reports / essaiPlot.R
CommitLineData
64b28e3e 1### Regression matrices
4c9cc558 2model = Res
64b28e3e 3K = dim(model$phi)[3]
4valMax = max(abs(model$phi))
5
6require(fields)
4c9cc558 7
64b28e3e 8if (K<4){
9 par(mfrow = c(1,K))
4c9cc558 10} else op = par(mfrow = c(2, (K+1)/2))
11
12## Phi
64b28e3e 13
14for (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 18par(mfrow = c(1,K),oma = c(0,0,3,0))
19mtext("Regression matrices in each cluster", side=3, line=4, font=2, cex=2, col='red')
20
21par(mfrow = c(1,2), oma=c(0,0,3,0))
22for (i in 1:4)
23 plot(runif(20), runif(20),
24 main=paste("random plot (",i,")",sep=''))
25par(op)
26mtext("Four plots",
27 side=3, line=4, font=2, cex=2, col='red')
64b28e3e 28
29### Zoom onto two classes we want to compare
30kSel = c(1,2)
31par(mfrow = c(1,3))
32
33for (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}
37image.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
41par(mfrow = c(K, 1))
42for (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
48Gam = matrix(0, ncol = K, nrow = n)
49gam = Gam
50for (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}
57affec = apply(gam, 1,which.max)
58gam2 = matrix(NA, ncol = K, nrow = n)
59for (i in 1:n){
60 gam2[i, affec[i]] = gam[i, affec[i]]
61}
62boxplot(gam2)
63
64### Mean in each cluster
65XY = cbind(X,Y)
66XY_class= list()
67meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2])
68for (r in 1:K){
69 XY_class[[r]] = XY[affec == r, ]
70 meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
71}
72
c0846e18 73matplot(meanPerClass, type='l')