reports/essaiPlot.R

   1 ### Regression matrices
   2 model = res_valse
   3 K = dim(model$phi)[3]
   4 valMax = max(abs(model$phi))
   5
   6 require(fields)
   7 if (K<4){
   8   par(mfrow = c(1,K))
   9 } else par(mfrow = c(2, (K+1)/2))
  10
  11 for (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
  17 kSel = c(1,2)
  18 par(mfrow = c(1,3))
  19
  20 for (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 }
  24 image.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
  28 par(mfrow = c(K, 1))
  29 for (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
  35 Gam = matrix(0, ncol = K, nrow = n)
  36 gam  = Gam
  37 for (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 }
  44 affec = apply(gam, 1,which.max)
  45 gam2 = matrix(NA, ncol = K, nrow = n)
  46 for (i in 1:n){
  47   gam2[i, affec[i]] = gam[i, affec[i]]
  48 }
  49 boxplot(gam2)
  50
  51 ### Mean in each cluster
  52 XY = cbind(X,Y)
  53 XY_class= list()
  54 meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2])
  55 for (r in 1:K){
  56   XY_class[[r]] = XY[affec == r, ]
  57   meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
  58 }
  59
  60 matplot(meanPerClass, type='l')