31444abc |
1 | a = sum(beta*c(1,gammai[i,])) |
2 | return(exp(a)/(1+exp(a))) |
3 | }) |
4 | U = runif(length(pydata)) |
5 | ydata = as.numeric(U<pydata) |
6 | plot(sort(pydata),col=ydata[order(pydata)]+2) |
7 | wi = Si %*% beta[-1] |
8 | plot(wi, type='l') |
9 | lines(1:15,rep(0,15)) |
10 | plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green") |
11 | for (ite in 1:30){ |
12 | print(ite) |
13 | Xdata = t(sapply(1:N,FUN = function(i){ |
14 | as.vector(Si%*%gammai[i,])+rnorm(15,0,s2) |
15 | })) |
16 | #beta = c(.1,.5,-.3) |
17 | beta = c(1,-3.5,4.5,-2.5) |
18 | pydata = sapply(1:N,FUN=function(i){ |
19 | a = sum(beta*c(1,gammai[i,])) |
20 | return(exp(a)/(1+exp(a))) |
21 | }) |
22 | U = runif(length(pydata)) |
23 | ydata = as.numeric(U<pydata) |
24 | lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green") |
25 | } |
26 | lines(1:15, rep(0,15), col='red') |
27 | Gamma |
28 | Gamma = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3) |
29 | s2 = 0.01 |
30 | Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3)) |
31 | gammai = t(sapply(1:N,FUN = function(i){ |
32 | rmvnorm(1,mugamma,Gamma) |
33 | })) |
34 | Xdata = t(sapply(1:N,FUN = function(i){ |
35 | as.vector(Si%*%gammai[i,])+rnorm(15,0,s2) |
36 | })) |
37 | matplot(t(Xdata),type="l") |
38 | #beta = c(.1,.5,-.3) |
39 | plot(sort(sapply(1:N,FUN=function(i) 1/(1+exp((-1)*(resEM$betah[1,ITfinal]+sum(resEM$betah[-1,ITfinal]*resEM$gammahat[i,]))))) )) |
40 | plot(wi, type='l') |
41 | wi = Si %*% beta[-1] |
42 | plot(wi, type='l') |
43 | lines(1:15,rep(0,15)) |
44 | plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green") |
45 | for (ite in 1:30){ |
46 | print(ite) |
47 | Xdata = t(sapply(1:N,FUN = function(i){ |
48 | as.vector(Si%*%gammai[i,])+rnorm(15,0,s2) |
49 | })) |
50 | #beta = c(.1,.5,-.3) |
51 | beta = c(1,-3.5,4.5,-2.5) |
52 | pydata = sapply(1:N,FUN=function(i){ |
53 | a = sum(beta*c(1,gammai[i,])) |
54 | return(exp(a)/(1+exp(a))) |
55 | }) |
56 | U = runif(length(pydata)) |
57 | ydata = as.numeric(U<pydata) |
58 | lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green") |
59 | } |
60 | lines(1:15, rep(0,15), col='red') |
61 | N = 200 |
62 | M = 100 |
63 | mugamma = c(1,2,3,2,4,6) |
64 | LGamma = matrix(0,nrow=6,ncol=6) |
65 | LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05) |
66 | LGamma = LGamma + diag(runif(6,0.2,1)) |
67 | Gamma = LGamma%*%t(LGamma) |
68 | s2 = 0.2 |
69 | freq = seq(0,1,length=M) |
70 | intervals = list(12:25,67:89) |
71 | Si = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3)))) |
72 | gammai = t(sapply(1:N,FUN = function(i){ |
73 | rmvnorm(1,mugamma,Gamma) |
74 | })) |
75 | Xdata = t(sapply(1:N,FUN = function(i){ |
76 | as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2) |
77 | })) |
78 | xdata = matrix(0,nrow=N,ncol=M) |
79 | xdata[,intervals[[1]] ] = Xdata[,1:14] |
80 | xdata[,intervals[[2]] ] = Xdata[,15:37] |
81 | beta = c(2,0.5,-2,1,0.2,1,-1) |
82 | pydata = sapply(1:N,FUN=function(i){ |
83 | a = sum(beta*c(1,gammai[i,])) |
84 | return(exp(a)/(1+exp(a))) |
85 | }) |
86 | U = runif(length(pydata)) |
87 | ydata = as.numeric(U<pydata) |
88 | plot(sort(pydata),col=ydata[order(pydata)]+1) |
89 | nBasis = c(3,3) |
90 | L = 20 |
91 | eer = rep(0,L) |
92 | for (ell in 1:L){ |
93 | kapp = sample(1:200,0.8*200,replace=FALSE) |
94 | xdata.app = xdata[kapp,] |
95 | ydata.app = ydata[kapp] |
96 | ktest = setdiff(1:200,kapp) |
97 | xdata.test = xdata[ktest,] |
98 | ydata.test = ydata[ktest] |
99 | resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE) |
100 | Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3)) |
101 | eer[ell] = mean(Yres$Ypred==ydata.test) |
102 | cat("\n", "ell =", ell) |
103 | } |
104 | source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R') |
105 | for (ell in 1:L){ |
106 | kapp = sample(1:200,0.8*200,replace=FALSE) |
107 | xdata.app = xdata[kapp,] |
108 | ydata.app = ydata[kapp] |
109 | ktest = setdiff(1:200,kapp) |
110 | xdata.test = xdata[ktest,] |
111 | ydata.test = ydata[ktest] |
112 | resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE) |
113 | Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3)) |
114 | eer[ell] = mean(Yres$Ypred==ydata.test) |
115 | cat("\n", "ell =", ell) |
116 | } |
117 | N = 3000 |
118 | mugamma = c(1,2,3) |
119 | Gamma = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3) |
120 | s2 = 0.01 |
121 | Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3)) |
122 | gammai = t(sapply(1:N,FUN = function(i){ |
123 | rmvnorm(1,mugamma,Gamma) |
124 | })) |
125 | Xdata = t(sapply(1:N,FUN = function(i){ |
126 | as.vector(Si%*%gammai[i,])+rnorm(15,0,s2) |
127 | })) |
128 | matplot(t(Xdata),type="l") |
129 | #beta = c(.1,.5,-.3) |
130 | beta = c(1,-3.5,4.5,-2.5) |
131 | pydata = sapply(1:N,FUN=function(i){ |
132 | a = sum(beta*c(1,gammai[i,])) |
133 | return(exp(a)/(1+exp(a))) |
134 | }) |
135 | U = runif(length(pydata)) |
136 | ydata = as.numeric(U<pydata) |
137 | plot(sort(pydata),col=ydata[order(pydata)]+2) |
138 | wi = Si %*% beta[-1] |
139 | plot(wi, type='l') |
140 | lines(1:15,rep(0,15)) |
141 | plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green") |
142 | for (ite in 1:30){ |
143 | print(ite) |
144 | Xdata = t(sapply(1:N,FUN = function(i){ |
145 | as.vector(Si%*%gammai[i,])+rnorm(15,0,s2) |
146 | })) |
147 | #beta = c(.1,.5,-.3) |
148 | beta = c(1,-3.5,4.5,-2.5) |
149 | pydata = sapply(1:N,FUN=function(i){ |
150 | a = sum(beta*c(1,gammai[i,])) |
151 | return(exp(a)/(1+exp(a))) |
152 | }) |
153 | U = runif(length(pydata)) |
154 | ydata = as.numeric(U<pydata) |
155 | lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green") |
156 | } |
157 | lines(1:15, rep(0,15), col='red') |
158 | N = 3000 |
159 | mugamma = c(1,2,3) |
160 | Gamma = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,5),3,3) |
161 | s2 = 0.01 |
162 | Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3)) |
163 | gammai = t(sapply(1:N,FUN = function(i){ |
164 | rmvnorm(1,mugamma,Gamma) |
165 | })) |
166 | Xdata = t(sapply(1:N,FUN = function(i){ |
167 | as.vector(Si%*%gammai[i,])+rnorm(15,0,s2) |
168 | })) |
169 | matplot(t(Xdata),type="l") |
170 | #beta = c(.1,.5,-.3) |
171 | beta = c(1,-3.5,4.5,-2.5) |
172 | pydata = sapply(1:N,FUN=function(i){ |
173 | a = sum(beta*c(1,gammai[i,])) |
174 | return(exp(a)/(1+exp(a))) |
175 | }) |
176 | U = runif(length(pydata)) |
177 | ydata = as.numeric(U<pydata) |
178 | plot(sort(pydata),col=ydata[order(pydata)]+2) |
179 | # ydata = ydata.noise |
180 | wi = Si %*% beta[-1] |
181 | plot(wi, type='l') |
182 | lines(1:15,rep(0,15)) |
183 | plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green") |
184 | for (ite in 1:30){ |
185 | print(ite) |
186 | Xdata = t(sapply(1:N,FUN = function(i){ |
187 | as.vector(Si%*%gammai[i,])+rnorm(15,0,s2) |
188 | })) |
189 | #beta = c(.1,.5,-.3) |
190 | beta = c(1,-3.5,4.5,-2.5) |
191 | pydata = sapply(1:N,FUN=function(i){ |
192 | a = sum(beta*c(1,gammai[i,])) |
193 | return(exp(a)/(1+exp(a))) |
194 | }) |
195 | U = runif(length(pydata)) |
196 | ydata = as.numeric(U<pydata) |
197 | lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green") |
198 | } |
199 | lines(1:15, rep(0,15), col='red') |
200 | matplot(t(Xdata[ydata==1,]),typ="l",lty=1,col="red") |
201 | matplot(t(Xdata[ydata==0,]),typ="l",lty=1,col="black",add=TRUE) |
202 | N = 200 |
203 | M = 100 |
204 | mugamma = c(1,2,3,2,4,6) |
205 | LGamma = matrix(0,nrow=6,ncol=6) |
206 | LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05) |
207 | LGamma = LGamma + diag(runif(6,0.2,1)) |
208 | Gamma = LGamma%*%t(LGamma) |
209 | s2 = 0.2 |
210 | freq = seq(0,1,length=M) |
211 | intervals = list(12:25,67:89) |
212 | Si = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3)))) |
213 | gammai = t(sapply(1:N,FUN = function(i){ |
214 | rmvnorm(1,mugamma,Gamma) |
215 | })) |
216 | Xdata = t(sapply(1:N,FUN = function(i){ |
217 | as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2) |
218 | })) |
219 | xdata = matrix(0,nrow=N,ncol=M) |
220 | xdata[,intervals[[1]] ] = Xdata[,1:14] |
221 | xdata[,intervals[[2]] ] = Xdata[,15:37] |
222 | beta = c(2,0.5,-2,1,0.2,1,-1) |
223 | pydata = sapply(1:N,FUN=function(i){ |
224 | a = sum(beta*c(1,gammai[i,])) |
225 | return(exp(a)/(1+exp(a))) |
226 | }) |
227 | U = runif(length(pydata)) |
228 | ydata = as.numeric(U<pydata) |
229 | plot(sort(pydata),col=ydata[order(pydata)]+1) |
230 | nBasis = c(3,3) |
231 | resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE) |
232 | source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R') |
233 | resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE) |
234 | L = 20 |
235 | eer = rep(0,L) |
236 | for (ell in 1:L){ |
237 | kapp = sample(1:200,0.8*200,replace=FALSE) |
238 | xdata.app = xdata[kapp,] |
239 | ydata.app = ydata[kapp] |
240 | ktest = setdiff(1:200,kapp) |
241 | xdata.test = xdata[ktest,] |
242 | ydata.test = ydata[ktest] |
243 | resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE) |
244 | Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3)) |
245 | eer[ell] = mean(Yres$Ypred==ydata.test) |
246 | cat("\n", "ell =", ell) |
247 | } |
248 | source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R') |
249 | for (ell in 1:L){ |
250 | kapp = sample(1:200,0.8*200,replace=FALSE) |
251 | xdata.app = xdata[kapp,] |
252 | ydata.app = ydata[kapp] |
253 | ktest = setdiff(1:200,kapp) |
254 | xdata.test = xdata[ktest,] |
255 | ydata.test = ydata[ktest] |
256 | resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE) |
257 | Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3)) |
258 | eer[ell] = mean(Yres$Ypred==ydata.test) |
259 | cat("\n", "ell =", ell) |
260 | } |
261 | eer |
262 | mean(eer) |
263 | source('~/Dropbox/projet spectres/cluster/EMLiqArtBandeHugues.R') |
264 | library(fda) |
265 | library(Matrix) |
266 | library(mvtnorm) |
267 | source('EMalgorithm.R') |
268 | source('EMLiqArtBandeHugues.R') |
269 | for (nVC in 1:100){ |
270 | EMLiqArtBandeHugues(nVC) |
271 | } |
272 | setwd("~/Dropbox/projet spectres/cluster") |
273 | library(fda) |
274 | library(Matrix) |
275 | library(mvtnorm) |
276 | source('EMalgorithm.R') |
277 | source('EMLiqArtBandeHugues.R') |
278 | for (nVC in 1:100){ |
279 | EMLiqArtBandeHugues(nVC) |
280 | } |
281 | sims= c(63:67,69,71:74,76:85,88:90,93:100) |
282 | #simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300)) |
283 | #for (i in 1:nrow(simulation_settings)){ |
284 | aa=cbind(sims,rep(9, each =length(sims))) |
285 | colnames(aa)=c("isim","model") |
286 | write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",") |
287 | setwd("~/Dropbox/WarpMixedModel/warpingMixedModel/SuperComputer") |
288 | write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",") |
289 | sims= 1:100 |
290 | colnames(sims) = "isim" |
291 | sims= 1:100 |
292 | #simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300)) |
293 | #for (i in 1:nrow(simulation_settings)){ |
294 | aa=cbind(sims,rep(10, each =length(sims))) |
295 | colnames(aa)=c("isim","model") |
296 | write.table(aa,file=paste("SimulationSettings10.txt",sep=""),row.names=FALSE,sep=",") |
297 | load("~/valse/data/data.RData") |
298 | load("~/valse/data/data.RData") |
299 | X |
300 | Y |
301 | library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3") |
302 | library(roxygen2) |
303 | setwd("~/valse") |
304 | document() |
305 | document() |
306 | setwd("~/") |
307 | install("valse") |
308 | c(9.75, |
309 | 20, |
310 | 11.75, |
311 | 11, |
312 | 8.25, |
313 | 12.5, |
314 | 11, |
315 | 18, |
316 | 7, |
317 | 12.75, |
318 | 13, |
319 | 14.75, |
320 | 4.75, |
321 | 11, |
322 | 20, |
323 | 13.5, |
324 | 8, |
325 | 13.25, |
326 | 6, |
327 | 17.5, |
328 | 13.25, |
329 | 8.5, |
330 | 9.5, |
331 | 16, |
332 | 8, |
333 | 9.5, |
334 | 10.25, |
335 | 13.75, |
336 | 9, |
337 | 14.75, |
338 | 12.5, |
339 | 19.5, |
340 | 17.5, |
341 | 7, |
342 | 11.5, |
343 | 4, |
344 | 7.5, |
345 | 13.25, |
346 | 10.5 |
347 | ) |
348 | notes = c(9.75, |
349 | 20, |
350 | 11.75, |
351 | 11, |
352 | 8.25, |
353 | 12.5, |
354 | 11, |
355 | 18, |
356 | 7, |
357 | 12.75, |
358 | 13, |
359 | 14.75, |
360 | 4.75, |
361 | 11, |
362 | 20, |
363 | 13.5, |
364 | 8, |
365 | 13.25, |
366 | 6, |
367 | 17.5, |
368 | 13.25, |
369 | 8.5, |
370 | 9.5, |
371 | 16, |
372 | 8, |
373 | 9.5, |
374 | 10.25, |
375 | 13.75, |
376 | 9, |
377 | 14.75, |
378 | 12.5, |
379 | 19.5, |
380 | 17.5, |
381 | 7, |
382 | 11.5, |
383 | 4, |
384 | 7.5, |
385 | 13.25, |
386 | 10.5 |
387 | ) |
388 | hist(notes) |
389 | hist(notes, nclass = 20) |
390 | notesBis = c(7.375,19.375,10.125,9,6,12.25,11,18.75,10.5,11.5,15.5,13,2.375,12.75,17.75,15.375,8.125,16,8.5,14.5,11.625,11.25,9.625,14, |
391 | 6.5,11.375,11.875,16.25,10.125,12,6.25,9.75,8.75,3.5,0,0,0,5.75,2,3.75,6.625,5.25) |
392 | hist(notesBis) |
393 | notes = c(7, |
394 | 19, |
395 | 10, |
396 | 9, |
397 | 6, |
398 | 12, |
399 | 11, |
400 | 19, |
401 | 11, |
402 | 12, |
403 | 16, |
404 | 13, |
405 | 2, |
406 | 13, |
407 | 18, |
408 | 15, |
409 | 8, |
410 | 16, |
411 | 9, |
412 | 15, |
413 | 12, |
414 | 11, |
415 | 10, |
416 | 14, |
417 | 7, |
418 | 11, |
419 | 12, |
420 | 16, |
421 | 13, |
422 | 15, |
423 | 14, |
424 | 19, |
425 | 14, |
426 | 4, 9, |
427 | 8, |
428 | 7, |
429 | 8, |
430 | 8 |
431 | ) |
432 | hist(notes) |
433 | shiny::runApp('Dropbox/cranview-master') |
434 | packages = 'shock' |
435 | min_date = Sys.Date() - 1 |
436 | for (pkg in packages) |
437 | { |
438 | # api data for package. we want the initial release - the first element of the "timeline" |
439 | pkg_data = httr::GET(paste0("http://crandb.r-pkg.org/", pkg, "/all")) |
440 | pkg_data = httr::content(pkg_data) |
441 | initial_release = pkg_data$timeline[[1]] |
442 | min_date = min(min_date, as.Date(initial_release)) |
443 | } |
444 | min_date |
445 | library("capushe", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3") |
446 | A = matrix(rnorm(25), ncol = 5) |
447 | A = rbind(A, matrix(0, ncol = 5, nrow = 2)) |
448 | A |
449 | A[rowSums(A)==0,]=c() |
450 | A[rowSums(A)==0,]=[] |
451 | A = A[rowSums(A)!=0,] |
452 | A |
453 | source('~/valse/pkg/R/valse.R') |
454 | setwd("~/valse/reports") |
455 | XY = cbind(X,Y) |
456 | X |
457 | p = 10 |
458 | q = 8 |
459 | k = 2 |
460 | D = 20 |
461 | meanX = rep(0,p) |
462 | covX = 0.1*diag(p) |
463 | covY = array(dim = c(q,q,k)) |
464 | covY[,,1] = 0.1*diag(q) |
465 | covY[,,2] = 0.2*diag(q) |
466 | beta = array(dim = c(p,q,2)) |
467 | beta[,,2] = matrix(c(rep(2,(D)),rep(0, p*q-D))) |
468 | beta[,,1] = matrix(c(rep(1,D),rep(0, p*q-D))) |
469 | n = 100 |
470 | pi = c(0.4,0.6) |
471 | data = generateXY(meanX,covX,covY, pi, beta, n) |
472 | source('~/valse/pkg/R/generateSampleInputs.R') |
473 | data = generateXY(meanX,covX,covY, pi, beta, n) |
474 | X = data$X |
475 | Y = data$Y |
476 | XY = cbind(X,Y) |
477 | data |
478 | data$class |
479 | affec = data$class |
480 | XY = cbind(X,Y) |
481 | for (r in 1:K){ |
482 | XY_class[[r]] = XY[affec == r, ] |
483 | } |
484 | K= 2 |
485 | for (r in 1:K){ |
486 | XY_class[[r]] = XY[affec == r, ] |
487 | } |
488 | XY_class= list() |
489 | for (r in 1:K){ |
490 | XY_class[[r]] = XY[affec == r, ] |
491 | } |
492 | for (r in 1:K){ |
493 | XY_class[[r]] = XY[affec == r, ] |
494 | meanPerClass[,r] = apply(XY_class[[r]], 2, mean) |
495 | } |
496 | meanPerClass= matrix() |
497 | for (r in 1:K){ |
498 | XY_class[[r]] = XY[affec == r, ] |
499 | meanPerClass[,r] = apply(XY_class[[r]], 2, mean) |
500 | } |
501 | apply(XY_class[[r]], 2, mean) |
502 | p |
503 | q |
504 | meanPerClass[r,] = apply(XY_class[[r]], 2, mean) |
505 | dim(XY) |
506 | meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2]) |
507 | for (r in 1:K){ |
508 | XY_class[[r]] = XY[affec == r, ] |
509 | meanPerClass[,r] = apply(XY_class[[r]], 2, mean) |
510 | } |
511 | matplot(meanPerClass) |
512 | matplot(meanPerClass, type='l') |