+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')
projects / valse.git / commitdiff