reports/.Rhistory

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