Adjustments for CRAN upload
[valse.git] / reports / essaiPlot.R
... / ...
CommitLineData
1### Regression matrices
2model = Res
3K = dim(model$phi)[3]
4valMax = max(abs(model$phi))
5
6require(fields)
7
8if (K<4){
9 par(mfrow = c(1,K))
10} else op = par(mfrow = c(2, (K+1)/2))
11
12## Phi
13
14for (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}
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')
28
29### Zoom onto two classes we want to compare
30kSel = c(1,2)
31par(mfrow = c(1,3))
32
33for (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}
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){
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}
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
73matplot(meanPerClass, type='l')