remove also NY21 as in R code
[valse.git] / reports / essaiPlot.R
... / ...
CommitLineData
1### Regression matrices
2model = res_valse
3K = dim(model$phi)[3]
4valMax = max(abs(model$phi))
5
6require(fields)
7if (K<4){
8 par(mfrow = c(1,K))
9} else par(mfrow = c(2, (K+1)/2))
10
11for (r in 1:K){
12 image.plot(t(abs(model$phi[,,r])),
13 col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
14}
15
16### Zoom onto two classes we want to compare
17kSel = c(1,2)
18par(mfrow = c(1,3))
19
20for (r in kSel){
21 image.plot(t(abs(model$phi[,,r])),xaxt="n",yaxt="n",
22 col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
23}
24image.plot(t(abs(model$phi[,,kSel[1]]-model$phi[,,kSel[2]])),
25 col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
26
27### Covariance matrices
28par(mfrow = c(K, 1))
29for (r in 1:K){
30 image.plot(matrix(diag(model$rho[,,r]), ncol= 1),
31 col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
32}
33
34### proportions
35Gam = matrix(0, ncol = K, nrow = n)
36gam = Gam
37for (i in 1:n){
38 for (r in 1:k){
39 sqNorm2 = sum( (Y[i,]%*%model$rho[,,r]-X[i,]%*%model$phi[,,r])^2 )
40 Gam[i,r] = model$pi[r] * exp(-0.5*sqNorm2)* det(model$rho[,,r])
41 }
42 gam[i,] = Gam[i,] / sum(Gam[i,])
43}
44affec = apply(gam, 1,which.max)
45gam2 = matrix(NA, ncol = K, nrow = n)
46for (i in 1:n){
47 gam2[i, affec[i]] = gam[i, affec[i]]
48}
49boxplot(gam2)
50
51### Mean in each cluster
52XY = cbind(X,Y)
53XY_class= list()
54meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2])
55for (r in 1:K){
56 XY_class[[r]] = XY[affec == r, ]
57 meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
58}
59
60matplot(meanPerClass, type='l')