Merge branch 'master' of auder.net:valse
[valse.git] / reports / essaiPlot.R
1 ### Regression matrices
2 model = Res
3 K = dim(model$phi)[3]
4 valMax = max(abs(model$phi))
5
6 require(fields)
7
8 if (K<4){
9 par(mfrow = c(1,K))
10 } else op = par(mfrow = c(2, (K+1)/2))
11
12 ## Phi
13
14 for (r in 1:K){
15 image.plot(t(abs(model$phi[,,r])),
16 col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
17 }
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')
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){
34 image.plot(t(abs(model$phi[,,r])),xaxt="n",yaxt="n",
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){
51 for (r in 1:K){
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
73 matplot(meanPerClass, type='l')