merge with remote

diff --git a/.git-fat b/.git-fat
new file mode 160000 (submodule)
index 0000000..286bcd3
--- /dev/null
+++ b/.git-fat
@@ -0,0 +1 @@
+Subproject commit 286bcd30cea5f86363a04a9313afaf9e2e7a7b81

diff --git a/.gitmodules b/.gitmodules
new file mode 100644 (file)
index 0000000..16826ed
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,6 @@
+[submodule ".git-fat"]
+       path = .git-fat
+       url = https://github.com/yagu0/git-fat.git
+[submodule ".nbstripout"]
+       path = .nbstripout
+       url = https://github.com/kynan/nbstripout.git

diff --git a/.nbstripout b/.nbstripout
new file mode 160000 (submodule)
index 0000000..10c80cb
--- /dev/null
+++ b/.nbstripout
@@ -0,0 +1 @@
+Subproject commit 10c80cbd5b6356fa38179ca526ee2894d9b3bd20

diff --git a/initialize.sh b/initialize.sh
new file mode 100755 (executable)
index 0000000..48f0ba1
--- /dev/null
+++ b/initialize.sh
@@ -0,0 +1,34 @@
+#!/bin/sh
+
+#initialize submodules, set-up .git/config and .gitattributes, and pre-push hook
+git submodule init && git submodule update --merge
+
+#filter for git-fat
+printf \
+'*.pdf filter=fat
+*.tar.xz filter=fat
+*.png filter=fat
+*.jpg filter=fat
+*.ps filter=fat\n' > .gitattributes
+
+#filter for Jupyter
+python .nbstripout/nbstripout.py --install --attributes .gitattributes
+
+#pre-commit and pre-push hooks: indentation, git fat push, submodules update
+cp hooks/* .git/hooks/
+
+#install formatR
+echo 'if (! "formatR" %in% rownames(installed.packages()))
+       install.packages("formatR",repos="https://cloud.r-project.org")' | R --slave
+
+#.gitfat file with remote on gitfat@auder.net
+printf '[rsync]\nremote = gitfat@auder.net:~/files/valse\n' > .gitfat
+
+#manual git-fat init: with relative path to binary
+#1] remove filter if exists http://stackoverflow.com/a/12179641/4640434
+sed -i '1N;$!N;s/\[filter "fat"\]\n.*\n.*//;P;D' .git/config
+#2] place new filter
+printf \
+'[filter "fat"]
+       clean = ./.git-fat/git-fat filter-clean
+       smudge = ./.git-fat/git-fat filter-smudge\n' >> .git/config

diff --git a/pkg/inst/testdata/TODO.csv b/pkg/inst/testdata/TODO.csv
new file mode 100644 (file)
index 0000000..d679966
--- /dev/null
+++ b/pkg/inst/testdata/TODO.csv
@@ -0,0 +1 @@
+ou alors data_test.RData, possible aussi

diff --git a/pkg/tests/testthat.R b/pkg/tests/testthat.R
new file mode 100644 (file)
index 0000000..88e5631
--- /dev/null
+++ b/pkg/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+library(valse) #ou load_all()
+
+test_check("valse")

diff --git a/pkg/tests/testthat/helper-context1.R b/pkg/tests/testthat/helper-context1.R
new file mode 100644 (file)
index 0000000..b40f358
--- /dev/null
+++ b/pkg/tests/testthat/helper-context1.R
@@ -0,0 +1,5 @@
+# Potential helpers for context 1
+help <- function()
+{
+       #...
+}

diff --git a/pkg/tests/testthat/test-context1.R b/pkg/tests/testthat/test-context1.R
new file mode 100644 (file)
index 0000000..17c633f
--- /dev/null
+++ b/pkg/tests/testthat/test-context1.R
@@ -0,0 +1,11 @@
+context("Context1")
+
+test_that("function 1...",
+{
+       #expect_lte( ..., ...)
+})
+
+test_that("function 2...",
+{
+       #expect_equal(..., ...)
+})

diff --git a/pkg/vignettes/.gitignore b/pkg/vignettes/.gitignore
new file mode 100644 (file)
index 0000000..e6493d4
--- /dev/null
+++ b/pkg/vignettes/.gitignore
@@ -0,0 +1,14 @@
+#ignore jupyter generated file (ipynb, HTML)
+*.html
+*.ipynb
+
+#and various (pdf)LaTeX files, in case of
+*.tex
+*.pdf
+*.aux
+*.dvi
+*.log
+*.out
+*.toc
+*.synctex.gz
+/figure/

diff --git a/reports/.Rhistory b/reports/.Rhistory
new file mode 100644 (file)
index 0000000..dc970cc
--- /dev/null
+++ b/reports/.Rhistory
@@ -0,0 +1,512 @@
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+Gamma
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+plot(sort(sapply(1:N,FUN=function(i) 1/(1+exp((-1)*(resEM$betah[1,ITfinal]+sum(resEM$betah[-1,ITfinal]*resEM$gammahat[i,]))))) ))
+plot(wi, type='l')
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+N       = 200
+M       = 100
+mugamma = c(1,2,3,2,4,6)
+LGamma  = matrix(0,nrow=6,ncol=6)
+LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
+LGamma  = LGamma + diag(runif(6,0.2,1))
+Gamma   = LGamma%*%t(LGamma)
+s2      = 0.2
+freq      = seq(0,1,length=M)
+intervals = list(12:25,67:89)
+Si        = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
+}))
+xdata = matrix(0,nrow=N,ncol=M)
+xdata[,intervals[[1]] ] = Xdata[,1:14]
+xdata[,intervals[[2]] ] = Xdata[,15:37]
+beta = c(2,0.5,-2,1,0.2,1,-1)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+1)
+nBasis    = c(3,3)
+L   = 20
+eer = rep(0,L)
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+N       = 3000
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+N       = 3000
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+# ydata = ydata.noise
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+matplot(t(Xdata[ydata==1,]),typ="l",lty=1,col="red")
+matplot(t(Xdata[ydata==0,]),typ="l",lty=1,col="black",add=TRUE)
+N       = 200
+M       = 100
+mugamma = c(1,2,3,2,4,6)
+LGamma  = matrix(0,nrow=6,ncol=6)
+LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
+LGamma  = LGamma + diag(runif(6,0.2,1))
+Gamma   = LGamma%*%t(LGamma)
+s2      = 0.2
+freq      = seq(0,1,length=M)
+intervals = list(12:25,67:89)
+Si        = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
+}))
+xdata = matrix(0,nrow=N,ncol=M)
+xdata[,intervals[[1]] ] = Xdata[,1:14]
+xdata[,intervals[[2]] ] = Xdata[,15:37]
+beta = c(2,0.5,-2,1,0.2,1,-1)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+1)
+nBasis    = c(3,3)
+resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
+L   = 20
+eer = rep(0,L)
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+eer
+mean(eer)
+source('~/Dropbox/projet spectres/cluster/EMLiqArtBandeHugues.R')
+library(fda)
+library(Matrix)
+library(mvtnorm)
+source('EMalgorithm.R')
+source('EMLiqArtBandeHugues.R')
+for (nVC in 1:100){
+EMLiqArtBandeHugues(nVC)
+}
+setwd("~/Dropbox/projet spectres/cluster")
+library(fda)
+library(Matrix)
+library(mvtnorm)
+source('EMalgorithm.R')
+source('EMLiqArtBandeHugues.R')
+for (nVC in 1:100){
+EMLiqArtBandeHugues(nVC)
+}
+sims= c(63:67,69,71:74,76:85,88:90,93:100)
+#simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300))
+#for (i in 1:nrow(simulation_settings)){
+aa=cbind(sims,rep(9, each =length(sims)))
+colnames(aa)=c("isim","model")
+write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
+setwd("~/Dropbox/WarpMixedModel/warpingMixedModel/SuperComputer")
+write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
+sims= 1:100
+colnames(sims) = "isim"
+sims= 1:100
+#simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300))
+#for (i in 1:nrow(simulation_settings)){
+aa=cbind(sims,rep(10, each =length(sims)))
+colnames(aa)=c("isim","model")
+write.table(aa,file=paste("SimulationSettings10.txt",sep=""),row.names=FALSE,sep=",")
+load("~/valse/data/data.RData")
+load("~/valse/data/data.RData")
+X
+Y
+library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
+library(roxygen2)
+setwd("~/valse")
+document()
+document()
+setwd("~/")
+install("valse")
+c(9.75,
+20,
+11.75,
+11,
+8.25,
+12.5,
+11,
+18,
+7,
+12.75,
+13,
+14.75,
+4.75,
+11,
+20,
+13.5,
+8,
+13.25,
+6,
+17.5,
+13.25,
+8.5,
+9.5,
+16,
+8,
+9.5,
+10.25,
+13.75,
+9,
+14.75,
+12.5,
+19.5,
+17.5,
+7,
+11.5,
+4,
+7.5,
+13.25,
+10.5
+)
+notes = c(9.75,
+20,
+11.75,
+11,
+8.25,
+12.5,
+11,
+18,
+7,
+12.75,
+13,
+14.75,
+4.75,
+11,
+20,
+13.5,
+8,
+13.25,
+6,
+17.5,
+13.25,
+8.5,
+9.5,
+16,
+8,
+9.5,
+10.25,
+13.75,
+9,
+14.75,
+12.5,
+19.5,
+17.5,
+7,
+11.5,
+4,
+7.5,
+13.25,
+10.5
+)
+hist(notes)
+hist(notes, nclass = 20)
+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,
+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)
+hist(notesBis)
+notes = c(7,
+19,
+10,
+9,
+6,
+12,
+11,
+19,
+11,
+12,
+16,
+13,
+2,
+13,
+18,
+15,
+8,
+16,
+9,
+15,
+12,
+11,
+10,
+14,
+7,
+11,
+12,
+16,
+13,
+15,
+14,
+19,
+14,
+4,  9,
+8,
+7,
+8,
+8
+)
+hist(notes)
+shiny::runApp('Dropbox/cranview-master')
+packages = 'shock'
+min_date = Sys.Date() - 1
+for (pkg in packages)
+{
+# api data for package. we want the initial release - the first element of the "timeline"
+pkg_data = httr::GET(paste0("http://crandb.r-pkg.org/", pkg, "/all"))
+pkg_data = httr::content(pkg_data)
+initial_release = pkg_data$timeline[[1]]
+min_date = min(min_date, as.Date(initial_release))
+}
+min_date
+library("capushe", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
+A = matrix(rnorm(25), ncol = 5)
+A = rbind(A, matrix(0, ncol = 5, nrow = 2))
+A
+A[rowSums(A)==0,]=c()
+A[rowSums(A)==0,]=[]
+A = A[rowSums(A)!=0,]
+A
+source('~/valse/pkg/R/valse.R')
+setwd("~/valse/reports")
+XY = cbind(X,Y)
+X
+p = 10
+q = 8
+k = 2
+D = 20
+meanX = rep(0,p)
+covX = 0.1*diag(p)
+covY = array(dim = c(q,q,k))
+covY[,,1] = 0.1*diag(q)
+covY[,,2] = 0.2*diag(q)
+beta = array(dim = c(p,q,2))
+beta[,,2] = matrix(c(rep(2,(D)),rep(0, p*q-D)))
+beta[,,1] = matrix(c(rep(1,D),rep(0, p*q-D)))
+n = 100
+pi = c(0.4,0.6)
+data = generateXY(meanX,covX,covY, pi, beta, n)
+source('~/valse/pkg/R/generateSampleInputs.R')
+data = generateXY(meanX,covX,covY, pi, beta, n)
+X = data$X
+Y = data$Y
+XY = cbind(X,Y)
+data
+data$class
+affec = data$class
+XY = cbind(X,Y)
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+}
+K= 2
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+}
+XY_class= list()
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+}
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
+}
+meanPerClass= matrix()
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
+}
+apply(XY_class[[r]], 2, mean)
+p
+q
+meanPerClass[r,] = apply(XY_class[[r]], 2, mean)
+dim(XY)
+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)
+matplot(meanPerClass, type='l')