reports/bazar Emilie/simulatedData.R

   1 library("mclust")
   2 #library("R.matlab", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
   3 #redrawData = TRUE
   4 #if (redrawData==TRUE){
   5   ###########
   6   ## Model
   7   ###########
   8   K = 2
   9   p = 48
  10   T = seq(0,1.5,length.out = p)
  11   T2 = seq(0,3, length.out = 2*p)
  12   n = 100
  13   x1 = cos(2*pi*T) + 0.2*cos(4*2*pi*T) +2*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7))))
  14   plot(T,x1)
  15   lines(T,x1)
  16
  17   sigmaX = 0.085
  18   sigmaY = 0.1
  19   beta = list()
  20   p1= 0.5
  21   beta[[1]] =diag(c(rep(p1,5),rep(1,5), rep(p1,5), rep(1, p-15)))
  22   p2 = 1
  23   beta[[2]] = diag(c(rep(p2,5),rep(1,5), rep(p2,5), rep(1, p-15)))
  24   ITE = 100
  25   ARI1 = ARI2 = ARI3 = rep(0,ITE)
  26   XY = array(0, dim = c(ITE, 2*p,n))
  27   XYproj = array(0, dim=c(ITE, 96,n))
  28
  29   affec = list()
  30   ###########
  31   ## Iterations
  32   ###########
  33   for (ite in c(1:ITE)){
  34     ###########
  35     ##Sample
  36     ###########
  37     x = x1 + matrix(rnorm(n*p, 0, sigmaX), ncol = n)
  38     affec[[ite]] = sample(c(1,2), n, replace = TRUE)
  39     y  = x
  40     xy = matrix(0,ncol=n, nrow= 2*p)
  41      for (i in c(1:n)){
  42        y[,i] = x[,i] %*% beta[[affec[[ite]][i]]] + rnorm(p, 0, sigmaY)
  43        xy[,i] = c(x[,i],y[,i])
  44        XY[ite,,i] = xy[,i] - mean(xy[,i])
  45     #   Dx = dwt(x[,i], filter='haar')@W
  46     #   Dx = rev(unlist(Dx))
  47     #   Dx = Dx[2:(1+3+6+12+24)]
  48     #   Ax = dwt(x[,i], filter='haar')@V
  49     #   Ax = rev(unlist(Ax))
  50     #   Ax = Ax[2:(1+3)]
  51     #   Dy = dwt(y[,i], filter='haar')@W
  52     #   Dy = rev(unlist(Dy))
  53     #   Dy = Dy[2:(1+3+6+12+24)]
  54     #   Ay = dwt(y[,i], filter='haar')@V
  55     #   Ay = rev(unlist(Ay))
  56     #   Ay = Ay[2:(1+3)]
  57     #   XYproj[ite,,i] = c(Ax,Dx,Ay,Dy)
  58      }
  59     print(ite)
  60     #
  61     #
  62   }
  63   xy[c(7,55),] = NA
  64  # write.table(XY,'data.csv', row.names=FALSE, col.names=FALSE)
  65 matplot(T2,xy[,affec[[ite]]==1],type='l', col='red', lty = 1)
  66 matplot(T2,xy[,affec[[ite]]==2],type='l', col='black', add=TRUE, lty= 1)
  67 abline(v = 1.5)
  68 text(0.75,0,'X', cex = 2 )
  69 text(0.75+1.5,0,'Y', cex = 2 )
  70 #proj =  read.table('dataProj.csv')
  71 #}
  72
  73
  74 #matplot(T,x,type='l', col='black', xlab = '', ylab='', lwd=1.5,lty=1)
  75 #matplot(T,y[,affec[[ite]]==1],type='l', col='red', xlab = '', ylab='', lwd=1.5,lty=1)
  76 #matplot(T,y[,affec[[ite]]==2],type='l', col='black', add=TRUE,lwd=2,lty=1)
  77 # proj2 = array(0,dim=c(ITE,2*p,n))
  78 # for (ite in c(1:ITE)){
  79 #   for (i in c(1:n)){
  80 #     A = proj[ite,(1+(i-1)*96):(i*96)]
  81 #     for (j in 1:96){
  82 #       proj2[ite,j,i] = A[1,j]
  83 #     }
  84 #   }
  85 #   print(ite)
  86 # }
  87 ###########
  88 ## Iterations
  89 ###########
  90 Kmod2 = Kmod1 = rep(0,ITE)
  91 Kmod3 = rep(0,ITE)
  92 for (ite in c(1:ITE)){
  93   print(ite)
  94   ###########
  95   ## k-means 1
  96   ###########
  97   mod1 = Mclust(t(XY[ite,,]),G = 2, mode='VII')
  98   ARI1[ite] = adjustedRandIndex(mod1$classification, affec[[ite]])
  99   Kmod1[ite] = mod1$G
 100   # ###########
 101   # ## k-means 2
 102   # ###########
 103   # #proj2 =
 104   # mod2 = Mclust(t(XYproj[ite,,]),G = 1:8, mode='VII')
 105   # ARI2[ite] = adjustedRandIndex(mod2$classification, affec[[ite]])
 106   # Kmod2[ite] = mod2$G
 107   # ###########
 108   # ## k-means 1
 109   # ###########
 110   # #proj3 =
 111   # mod3 = Mclust(t(XYproj[ite,c(4:12,52:60),]),G = 1:8, mode='VII')
 112   # ARI3[ite] = adjustedRandIndex(mod3$classification, affec[[ite]])
 113   # Kmod3[ite] = mod3$G
 114 }
 115 ARI0 = rep(1,ITE)
 116 par(cex.lab=1.5)
 117 par(cex.axis=1.5)
 118 boxplot(ARI0,ARI1, names = c('LassoMLE','K-means'), lwd=1.3)
 119 table(Kmod1)
 120 table(Kmod2)