.Rhistory

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