X-Git-Url: https://git.auder.net/?p=valse.git;a=blobdiff_plain;f=reports%2FessaiPlot.R;h=dea7e00ebee921f85e8a62957755947baeb2b18f;hp=e69de29bb2d1d6434b8b29ae775ad8c2e48c5391;hb=64b28e3edeef11b4442b6014ec89246810ebc1cf;hpb=3f62d540d32b70f42b9090087f72426c18cb219e

diff --git a/reports/essaiPlot.R b/reports/essaiPlot.R
index e69de29..dea7e00 100644
--- a/reports/essaiPlot.R
+++ b/reports/essaiPlot.R
@@ -0,0 +1,60 @@
+### Regression matrices
+model = res_valse
+K = dim(model$phi)[3]
+valMax = max(abs(model$phi))
+
+require(fields)
+if (K<4){
+  par(mfrow = c(1,K))
+} else par(mfrow = c(2, (K+1)/2))
+
+for (r in 1:K){
+  image.plot(t(abs(model$phi[,,r])), 
+             col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
+}
+
+### Zoom onto two classes we want to compare 
+kSel = c(1,2)
+par(mfrow = c(1,3))
+
+for (r in kSel){
+  image.plot(t(abs(model$phi[,,r])),title="hat{beta}",xaxt="n",yaxt="n", 
+             col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
+}
+image.plot(t(abs(model$phi[,,kSel[1]]-model$phi[,,kSel[2]])), 
+           col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
+
+### Covariance matrices
+par(mfrow = c(K, 1))
+for (r in 1:K){
+  image.plot(matrix(diag(model$rho[,,r]), ncol= 1), 
+             col=gray(rev(seq(0,64,length.out=65))/65),breaks=seq(0,valMax,length.out=66))
+}
+
+### proportions
+Gam = matrix(0, ncol = K, nrow = n)
+gam  = Gam
+for (i in 1:n){
+  for (r in 1:k){
+    sqNorm2 = sum( (Y[i,]%*%model$rho[,,r]-X[i,]%*%model$phi[,,r])^2 )
+    Gam[i,r] = model$pi[r] * exp(-0.5*sqNorm2)* det(model$rho[,,r])
+  }
+  gam[i,] = Gam[i,] / sum(Gam[i,])
+}
+affec = apply(gam, 1,which.max)
+gam2 = matrix(NA, ncol = K, nrow = n)
+for (i in 1:n){
+  gam2[i, affec[i]] = gam[i, affec[i]]
+}
+boxplot(gam2)
+
+### Mean in each cluster
+XY = cbind(X,Y)
+XY_class= list()
+meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2])
+for (r in 1:K){
+  XY_class[[r]] = XY[affec == r, ]
+  meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
+}
+
+matplot(meanPerClass, type='l')
\ No newline at end of file