| 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') |