basic test OK, but problem with 'fast' EMGLLF output in last EMilie test
[valse.git] / 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*base::pi*T) + 0.2*cos(4*2*base::pi*T) + 0.3*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7))))+1
14 plot(T,x1)
15 lines(T,x1)
16
17 sigmaX = 0.12
18 sigmaY = 0.12
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 = 1: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)