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')