lines(1:15, rep(0,15), col='red')
N       = 300
mugamma = c(1,2,3)
Gamma   = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,1),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')
matplot(Xdata)
matplot(Xdata,type='l')
matplot(t(Xdata),type='l')
matplot(apply(Xdata, 2, mean),type='l')
mean(apply(Xdata, 2, mean))
N       = 3000
mugamma = c(1,2,3)
Gamma   = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,1),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')
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
load("~/valse/data/data.RData")
k = 2
init = initSmallEM(k, X, Y)
setwd("~/valse")
library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
library("roxygen2", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")