update EMGLLF.R and some few details
[valse.git] / .Rhistory
CommitLineData
f227455a 1lines(1:15, rep(0,15), col='red')
2N = 300
3mugamma = c(1,2,3)
4Gamma = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,1),3,3)
5s2 = 0.01
6Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
7gammai = t(sapply(1:N,FUN = function(i){
8rmvnorm(1,mugamma,Gamma)
9}))
10Xdata = t(sapply(1:N,FUN = function(i){
11as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
12}))
13matplot(t(Xdata),type="l")
14#beta = c(.1,.5,-.3)
15beta = c(1,-3.5,4.5,-2.5)
16pydata = sapply(1:N,FUN=function(i){
17a = sum(beta*c(1,gammai[i,]))
18return(exp(a)/(1+exp(a)))
19})
20U = runif(length(pydata))
21ydata = as.numeric(U<pydata)
22plot(sort(pydata),col=ydata[order(pydata)]+2)
23wi = Si %*% beta[-1]
24plot(wi, type='l')
25lines(1:15,rep(0,15))
26plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
27for (ite in 1:30){
28print(ite)
29Xdata = t(sapply(1:N,FUN = function(i){
30as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
31}))
32#beta = c(.1,.5,-.3)
33beta = c(1,-3.5,4.5,-2.5)
34pydata = sapply(1:N,FUN=function(i){
35a = sum(beta*c(1,gammai[i,]))
36return(exp(a)/(1+exp(a)))
37})
38U = runif(length(pydata))
39ydata = as.numeric(U<pydata)
40lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
41}
42lines(1:15, rep(0,15), col='red')
43matplot(Xdata)
44matplot(Xdata,type='l')
45matplot(t(Xdata),type='l')
46matplot(apply(Xdata, 2, mean),type='l')
47mean(apply(Xdata, 2, mean))
48N = 3000
49mugamma = c(1,2,3)
50Gamma = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,1),3,3)
51s2 = 0.01
52Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
53gammai = t(sapply(1:N,FUN = function(i){
54rmvnorm(1,mugamma,Gamma)
55}))
56Xdata = t(sapply(1:N,FUN = function(i){
57as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
58}))
59matplot(t(Xdata),type="l")
60#beta = c(.1,.5,-.3)
61beta = c(1,-3.5,4.5,-2.5)
62pydata = sapply(1:N,FUN=function(i){
63a = sum(beta*c(1,gammai[i,]))
64return(exp(a)/(1+exp(a)))
65})
66U = runif(length(pydata))
67ydata = as.numeric(U<pydata)
68plot(sort(pydata),col=ydata[order(pydata)]+2)
69wi = Si %*% beta[-1]
70plot(wi, type='l')
71lines(1:15,rep(0,15))
72plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
73for (ite in 1:30){
74print(ite)
75Xdata = t(sapply(1:N,FUN = function(i){
76as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
77}))
78#beta = c(.1,.5,-.3)
79beta = c(1,-3.5,4.5,-2.5)
80pydata = sapply(1:N,FUN=function(i){
81a = sum(beta*c(1,gammai[i,]))
82return(exp(a)/(1+exp(a)))
83})
84U = runif(length(pydata))
85ydata = as.numeric(U<pydata)
86lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
87}
88lines(1:15, rep(0,15), col='red')
89Gamma
90Gamma = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
91s2 = 0.01
92Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
93gammai = t(sapply(1:N,FUN = function(i){
94rmvnorm(1,mugamma,Gamma)
95}))
96Xdata = t(sapply(1:N,FUN = function(i){
97as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
98}))
99matplot(t(Xdata),type="l")
100#beta = c(.1,.5,-.3)
101plot(sort(sapply(1:N,FUN=function(i) 1/(1+exp((-1)*(resEM$betah[1,ITfinal]+sum(resEM$betah[-1,ITfinal]*resEM$gammahat[i,]))))) ))
102plot(wi, type='l')
103wi = Si %*% beta[-1]
104plot(wi, type='l')
105lines(1:15,rep(0,15))
106plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
107for (ite in 1:30){
108print(ite)
109Xdata = t(sapply(1:N,FUN = function(i){
110as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
111}))
112#beta = c(.1,.5,-.3)
113beta = c(1,-3.5,4.5,-2.5)
114pydata = sapply(1:N,FUN=function(i){
115a = sum(beta*c(1,gammai[i,]))
116return(exp(a)/(1+exp(a)))
117})
118U = runif(length(pydata))
119ydata = as.numeric(U<pydata)
120lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
121}
122lines(1:15, rep(0,15), col='red')
123N = 200
124M = 100
125mugamma = c(1,2,3,2,4,6)
126LGamma = matrix(0,nrow=6,ncol=6)
127LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
128LGamma = LGamma + diag(runif(6,0.2,1))
129Gamma = LGamma%*%t(LGamma)
130s2 = 0.2
131freq = seq(0,1,length=M)
132intervals = list(12:25,67:89)
133Si = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
134gammai = t(sapply(1:N,FUN = function(i){
135rmvnorm(1,mugamma,Gamma)
136}))
137Xdata = t(sapply(1:N,FUN = function(i){
138as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
139}))
140xdata = matrix(0,nrow=N,ncol=M)
141xdata[,intervals[[1]] ] = Xdata[,1:14]
142xdata[,intervals[[2]] ] = Xdata[,15:37]
143beta = c(2,0.5,-2,1,0.2,1,-1)
144pydata = sapply(1:N,FUN=function(i){
145a = sum(beta*c(1,gammai[i,]))
146return(exp(a)/(1+exp(a)))
147})
148U = runif(length(pydata))
149ydata = as.numeric(U<pydata)
150plot(sort(pydata),col=ydata[order(pydata)]+1)
151nBasis = c(3,3)
152L = 20
153eer = rep(0,L)
154for (ell in 1:L){
155kapp = sample(1:200,0.8*200,replace=FALSE)
156xdata.app = xdata[kapp,]
157ydata.app = ydata[kapp]
158ktest = setdiff(1:200,kapp)
159xdata.test = xdata[ktest,]
160ydata.test = ydata[ktest]
161resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
162Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
163eer[ell] = mean(Yres$Ypred==ydata.test)
164cat("\n", "ell =", ell)
165}
166source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
167for (ell in 1:L){
168kapp = sample(1:200,0.8*200,replace=FALSE)
169xdata.app = xdata[kapp,]
170ydata.app = ydata[kapp]
171ktest = setdiff(1:200,kapp)
172xdata.test = xdata[ktest,]
173ydata.test = ydata[ktest]
174resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
175Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
176eer[ell] = mean(Yres$Ypred==ydata.test)
177cat("\n", "ell =", ell)
178}
179N = 3000
180mugamma = c(1,2,3)
181Gamma = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
182s2 = 0.01
183Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
184gammai = t(sapply(1:N,FUN = function(i){
185rmvnorm(1,mugamma,Gamma)
186}))
187Xdata = t(sapply(1:N,FUN = function(i){
188as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
189}))
190matplot(t(Xdata),type="l")
191#beta = c(.1,.5,-.3)
192beta = c(1,-3.5,4.5,-2.5)
193pydata = sapply(1:N,FUN=function(i){
194a = sum(beta*c(1,gammai[i,]))
195return(exp(a)/(1+exp(a)))
196})
197U = runif(length(pydata))
198ydata = as.numeric(U<pydata)
199plot(sort(pydata),col=ydata[order(pydata)]+2)
200wi = Si %*% beta[-1]
201plot(wi, type='l')
202lines(1:15,rep(0,15))
203plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
204for (ite in 1:30){
205print(ite)
206Xdata = t(sapply(1:N,FUN = function(i){
207as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
208}))
209#beta = c(.1,.5,-.3)
210beta = c(1,-3.5,4.5,-2.5)
211pydata = sapply(1:N,FUN=function(i){
212a = sum(beta*c(1,gammai[i,]))
213return(exp(a)/(1+exp(a)))
214})
215U = runif(length(pydata))
216ydata = as.numeric(U<pydata)
217lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
218}
219lines(1:15, rep(0,15), col='red')
220N = 3000
221mugamma = c(1,2,3)
222Gamma = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,5),3,3)
223s2 = 0.01
224Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
225gammai = t(sapply(1:N,FUN = function(i){
226rmvnorm(1,mugamma,Gamma)
227}))
228Xdata = t(sapply(1:N,FUN = function(i){
229as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
230}))
231matplot(t(Xdata),type="l")
232#beta = c(.1,.5,-.3)
233beta = c(1,-3.5,4.5,-2.5)
234pydata = sapply(1:N,FUN=function(i){
235a = sum(beta*c(1,gammai[i,]))
236return(exp(a)/(1+exp(a)))
237})
238U = runif(length(pydata))
239ydata = as.numeric(U<pydata)
240plot(sort(pydata),col=ydata[order(pydata)]+2)
241# ydata = ydata.noise
242wi = Si %*% beta[-1]
243plot(wi, type='l')
244lines(1:15,rep(0,15))
245plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
246for (ite in 1:30){
247print(ite)
248Xdata = t(sapply(1:N,FUN = function(i){
249as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
250}))
251#beta = c(.1,.5,-.3)
252beta = c(1,-3.5,4.5,-2.5)
253pydata = sapply(1:N,FUN=function(i){
254a = sum(beta*c(1,gammai[i,]))
255return(exp(a)/(1+exp(a)))
256})
257U = runif(length(pydata))
258ydata = as.numeric(U<pydata)
259lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
260}
261lines(1:15, rep(0,15), col='red')
262matplot(t(Xdata[ydata==1,]),typ="l",lty=1,col="red")
263matplot(t(Xdata[ydata==0,]),typ="l",lty=1,col="black",add=TRUE)
264N = 200
265M = 100
266mugamma = c(1,2,3,2,4,6)
267LGamma = matrix(0,nrow=6,ncol=6)
268LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
269LGamma = LGamma + diag(runif(6,0.2,1))
270Gamma = LGamma%*%t(LGamma)
271s2 = 0.2
272freq = seq(0,1,length=M)
273intervals = list(12:25,67:89)
274Si = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
275gammai = t(sapply(1:N,FUN = function(i){
276rmvnorm(1,mugamma,Gamma)
277}))
278Xdata = t(sapply(1:N,FUN = function(i){
279as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
280}))
281xdata = matrix(0,nrow=N,ncol=M)
282xdata[,intervals[[1]] ] = Xdata[,1:14]
283xdata[,intervals[[2]] ] = Xdata[,15:37]
284beta = c(2,0.5,-2,1,0.2,1,-1)
285pydata = sapply(1:N,FUN=function(i){
286a = sum(beta*c(1,gammai[i,]))
287return(exp(a)/(1+exp(a)))
288})
289U = runif(length(pydata))
290ydata = as.numeric(U<pydata)
291plot(sort(pydata),col=ydata[order(pydata)]+1)
292nBasis = c(3,3)
293resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
294source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
295resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
296L = 20
297eer = rep(0,L)
298for (ell in 1:L){
299kapp = sample(1:200,0.8*200,replace=FALSE)
300xdata.app = xdata[kapp,]
301ydata.app = ydata[kapp]
302ktest = setdiff(1:200,kapp)
303xdata.test = xdata[ktest,]
304ydata.test = ydata[ktest]
305resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
306Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
307eer[ell] = mean(Yres$Ypred==ydata.test)
308cat("\n", "ell =", ell)
309}
310source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
311for (ell in 1:L){
312kapp = sample(1:200,0.8*200,replace=FALSE)
313xdata.app = xdata[kapp,]
314ydata.app = ydata[kapp]
315ktest = setdiff(1:200,kapp)
316xdata.test = xdata[ktest,]
317ydata.test = ydata[ktest]
318resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
319Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
320eer[ell] = mean(Yres$Ypred==ydata.test)
321cat("\n", "ell =", ell)
322}
323eer
324mean(eer)
325source('~/Dropbox/projet spectres/cluster/EMLiqArtBandeHugues.R')
326library(fda)
327library(Matrix)
328library(mvtnorm)
329source('EMalgorithm.R')
330source('EMLiqArtBandeHugues.R')
331for (nVC in 1:100){
332EMLiqArtBandeHugues(nVC)
333}
334setwd("~/Dropbox/projet spectres/cluster")
335library(fda)
336library(Matrix)
337library(mvtnorm)
338source('EMalgorithm.R')
339source('EMLiqArtBandeHugues.R')
340for (nVC in 1:100){
341EMLiqArtBandeHugues(nVC)
342}
343sims= 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)){
346aa=cbind(sims,rep(9, each =length(sims)))
347colnames(aa)=c("isim","model")
348write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
349setwd("~/Dropbox/WarpMixedModel/warpingMixedModel/SuperComputer")
350write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
351sims= 1:100
352colnames(sims) = "isim"
353sims= 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)){
356aa=cbind(sims,rep(10, each =length(sims)))
357colnames(aa)=c("isim","model")
358write.table(aa,file=paste("SimulationSettings10.txt",sep=""),row.names=FALSE,sep=",")
359load("~/valse/data/data.RData")
360load("~/valse/data/data.RData")
361X
362Y
363library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
364library(roxygen2)
365setwd("~/valse")
366document()
367document()
368setwd("~/")
369install("valse")
370c(9.75,
37120,
37211.75,
37311,
3748.25,
37512.5,
37611,
37718,
3787,
37912.75,
38013,
38114.75,
3824.75,
38311,
38420,
38513.5,
3868,
38713.25,
3886,
38917.5,
39013.25,
3918.5,
3929.5,
39316,
3948,
3959.5,
39610.25,
39713.75,
3989,
39914.75,
40012.5,
40119.5,
40217.5,
4037,
40411.5,
4054,
4067.5,
40713.25,
40810.5
409)
410notes = c(9.75,
41120,
41211.75,
41311,
4148.25,
41512.5,
41611,
41718,
4187,
41912.75,
42013,
42114.75,
4224.75,
42311,
42420,
42513.5,
4268,
42713.25,
4286,
42917.5,
43013.25,
4318.5,
4329.5,
43316,
4348,
4359.5,
43610.25,
43713.75,
4389,
43914.75,
44012.5,
44119.5,
44217.5,
4437,
44411.5,
4454,
4467.5,
44713.25,
44810.5
449)
450hist(notes)
451hist(notes, nclass = 20)
452notesBis = 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,
4536.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)
454hist(notesBis)
455notes = c(7,
45619,
45710,
4589,
4596,
46012,
46111,
46219,
46311,
46412,
46516,
46613,
4672,
46813,
46918,
47015,
4718,
47216,
4739,
47415,
47512,
47611,
47710,
47814,
4797,
48011,
48112,
48216,
48313,
48415,
48514,
48619,
48714,
4884, 9,
4898,
4907,
4918,
4928
493)
494hist(notes)
495shiny::runApp('Dropbox/cranview-master')
496packages = 'shock'
497min_date = Sys.Date() - 1
498for (pkg in packages)
499{
500# api data for package. we want the initial release - the first element of the "timeline"
501pkg_data = httr::GET(paste0("http://crandb.r-pkg.org/", pkg, "/all"))
502pkg_data = httr::content(pkg_data)
503initial_release = pkg_data$timeline[[1]]
504min_date = min(min_date, as.Date(initial_release))
505}
506min_date
507load("~/valse/data/data.RData")
508k = 2
509init = initSmallEM(k, X, Y)
510setwd("~/valse")
511library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
512library("roxygen2", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")