Fix test/ folder: tests of C code run fine now
[valse.git] / reports / .Rhistory
CommitLineData
31444abc 1a = sum(beta*c(1,gammai[i,]))
2return(exp(a)/(1+exp(a)))
3})
4U = runif(length(pydata))
5ydata = as.numeric(U<pydata)
6plot(sort(pydata),col=ydata[order(pydata)]+2)
7wi = Si %*% beta[-1]
8plot(wi, type='l')
9lines(1:15,rep(0,15))
10plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
11for (ite in 1:30){
12print(ite)
13Xdata = t(sapply(1:N,FUN = function(i){
14as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
15}))
16#beta = c(.1,.5,-.3)
17beta = c(1,-3.5,4.5,-2.5)
18pydata = sapply(1:N,FUN=function(i){
19a = sum(beta*c(1,gammai[i,]))
20return(exp(a)/(1+exp(a)))
21})
22U = runif(length(pydata))
23ydata = as.numeric(U<pydata)
24lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
25}
26lines(1:15, rep(0,15), col='red')
27Gamma
28Gamma = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
29s2 = 0.01
30Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
31gammai = t(sapply(1:N,FUN = function(i){
32rmvnorm(1,mugamma,Gamma)
33}))
34Xdata = t(sapply(1:N,FUN = function(i){
35as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
36}))
37matplot(t(Xdata),type="l")
38#beta = c(.1,.5,-.3)
39plot(sort(sapply(1:N,FUN=function(i) 1/(1+exp((-1)*(resEM$betah[1,ITfinal]+sum(resEM$betah[-1,ITfinal]*resEM$gammahat[i,]))))) ))
40plot(wi, type='l')
41wi = Si %*% beta[-1]
42plot(wi, type='l')
43lines(1:15,rep(0,15))
44plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
45for (ite in 1:30){
46print(ite)
47Xdata = t(sapply(1:N,FUN = function(i){
48as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
49}))
50#beta = c(.1,.5,-.3)
51beta = c(1,-3.5,4.5,-2.5)
52pydata = sapply(1:N,FUN=function(i){
53a = sum(beta*c(1,gammai[i,]))
54return(exp(a)/(1+exp(a)))
55})
56U = runif(length(pydata))
57ydata = as.numeric(U<pydata)
58lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
59}
60lines(1:15, rep(0,15), col='red')
61N = 200
62M = 100
63mugamma = c(1,2,3,2,4,6)
64LGamma = matrix(0,nrow=6,ncol=6)
65LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
66LGamma = LGamma + diag(runif(6,0.2,1))
67Gamma = LGamma%*%t(LGamma)
68s2 = 0.2
69freq = seq(0,1,length=M)
70intervals = list(12:25,67:89)
71Si = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
72gammai = t(sapply(1:N,FUN = function(i){
73rmvnorm(1,mugamma,Gamma)
74}))
75Xdata = t(sapply(1:N,FUN = function(i){
76as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
77}))
78xdata = matrix(0,nrow=N,ncol=M)
79xdata[,intervals[[1]] ] = Xdata[,1:14]
80xdata[,intervals[[2]] ] = Xdata[,15:37]
81beta = c(2,0.5,-2,1,0.2,1,-1)
82pydata = sapply(1:N,FUN=function(i){
83a = sum(beta*c(1,gammai[i,]))
84return(exp(a)/(1+exp(a)))
85})
86U = runif(length(pydata))
87ydata = as.numeric(U<pydata)
88plot(sort(pydata),col=ydata[order(pydata)]+1)
89nBasis = c(3,3)
90L = 20
91eer = rep(0,L)
92for (ell in 1:L){
93kapp = sample(1:200,0.8*200,replace=FALSE)
94xdata.app = xdata[kapp,]
95ydata.app = ydata[kapp]
96ktest = setdiff(1:200,kapp)
97xdata.test = xdata[ktest,]
98ydata.test = ydata[ktest]
99resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
100Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
101eer[ell] = mean(Yres$Ypred==ydata.test)
102cat("\n", "ell =", ell)
103}
104source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
105for (ell in 1:L){
106kapp = sample(1:200,0.8*200,replace=FALSE)
107xdata.app = xdata[kapp,]
108ydata.app = ydata[kapp]
109ktest = setdiff(1:200,kapp)
110xdata.test = xdata[ktest,]
111ydata.test = ydata[ktest]
112resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
113Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
114eer[ell] = mean(Yres$Ypred==ydata.test)
115cat("\n", "ell =", ell)
116}
117N = 3000
118mugamma = c(1,2,3)
119Gamma = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
120s2 = 0.01
121Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
122gammai = t(sapply(1:N,FUN = function(i){
123rmvnorm(1,mugamma,Gamma)
124}))
125Xdata = t(sapply(1:N,FUN = function(i){
126as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
127}))
128matplot(t(Xdata),type="l")
129#beta = c(.1,.5,-.3)
130beta = c(1,-3.5,4.5,-2.5)
131pydata = sapply(1:N,FUN=function(i){
132a = sum(beta*c(1,gammai[i,]))
133return(exp(a)/(1+exp(a)))
134})
135U = runif(length(pydata))
136ydata = as.numeric(U<pydata)
137plot(sort(pydata),col=ydata[order(pydata)]+2)
138wi = Si %*% beta[-1]
139plot(wi, type='l')
140lines(1:15,rep(0,15))
141plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
142for (ite in 1:30){
143print(ite)
144Xdata = t(sapply(1:N,FUN = function(i){
145as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
146}))
147#beta = c(.1,.5,-.3)
148beta = c(1,-3.5,4.5,-2.5)
149pydata = sapply(1:N,FUN=function(i){
150a = sum(beta*c(1,gammai[i,]))
151return(exp(a)/(1+exp(a)))
152})
153U = runif(length(pydata))
154ydata = as.numeric(U<pydata)
155lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
156}
157lines(1:15, rep(0,15), col='red')
158N = 3000
159mugamma = c(1,2,3)
160Gamma = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,5),3,3)
161s2 = 0.01
162Si = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
163gammai = t(sapply(1:N,FUN = function(i){
164rmvnorm(1,mugamma,Gamma)
165}))
166Xdata = t(sapply(1:N,FUN = function(i){
167as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
168}))
169matplot(t(Xdata),type="l")
170#beta = c(.1,.5,-.3)
171beta = c(1,-3.5,4.5,-2.5)
172pydata = sapply(1:N,FUN=function(i){
173a = sum(beta*c(1,gammai[i,]))
174return(exp(a)/(1+exp(a)))
175})
176U = runif(length(pydata))
177ydata = as.numeric(U<pydata)
178plot(sort(pydata),col=ydata[order(pydata)]+2)
179# ydata = ydata.noise
180wi = Si %*% beta[-1]
181plot(wi, type='l')
182lines(1:15,rep(0,15))
183plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
184for (ite in 1:30){
185print(ite)
186Xdata = t(sapply(1:N,FUN = function(i){
187as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
188}))
189#beta = c(.1,.5,-.3)
190beta = c(1,-3.5,4.5,-2.5)
191pydata = sapply(1:N,FUN=function(i){
192a = sum(beta*c(1,gammai[i,]))
193return(exp(a)/(1+exp(a)))
194})
195U = runif(length(pydata))
196ydata = as.numeric(U<pydata)
197lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
198}
199lines(1:15, rep(0,15), col='red')
200matplot(t(Xdata[ydata==1,]),typ="l",lty=1,col="red")
201matplot(t(Xdata[ydata==0,]),typ="l",lty=1,col="black",add=TRUE)
202N = 200
203M = 100
204mugamma = c(1,2,3,2,4,6)
205LGamma = matrix(0,nrow=6,ncol=6)
206LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
207LGamma = LGamma + diag(runif(6,0.2,1))
208Gamma = LGamma%*%t(LGamma)
209s2 = 0.2
210freq = seq(0,1,length=M)
211intervals = list(12:25,67:89)
212Si = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
213gammai = t(sapply(1:N,FUN = function(i){
214rmvnorm(1,mugamma,Gamma)
215}))
216Xdata = t(sapply(1:N,FUN = function(i){
217as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
218}))
219xdata = matrix(0,nrow=N,ncol=M)
220xdata[,intervals[[1]] ] = Xdata[,1:14]
221xdata[,intervals[[2]] ] = Xdata[,15:37]
222beta = c(2,0.5,-2,1,0.2,1,-1)
223pydata = sapply(1:N,FUN=function(i){
224a = sum(beta*c(1,gammai[i,]))
225return(exp(a)/(1+exp(a)))
226})
227U = runif(length(pydata))
228ydata = as.numeric(U<pydata)
229plot(sort(pydata),col=ydata[order(pydata)]+1)
230nBasis = c(3,3)
231resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
232source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
233resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
234L = 20
235eer = rep(0,L)
236for (ell in 1:L){
237kapp = sample(1:200,0.8*200,replace=FALSE)
238xdata.app = xdata[kapp,]
239ydata.app = ydata[kapp]
240ktest = setdiff(1:200,kapp)
241xdata.test = xdata[ktest,]
242ydata.test = ydata[ktest]
243resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
244Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
245eer[ell] = mean(Yres$Ypred==ydata.test)
246cat("\n", "ell =", ell)
247}
248source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
249for (ell in 1:L){
250kapp = sample(1:200,0.8*200,replace=FALSE)
251xdata.app = xdata[kapp,]
252ydata.app = ydata[kapp]
253ktest = setdiff(1:200,kapp)
254xdata.test = xdata[ktest,]
255ydata.test = ydata[ktest]
256resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
257Yres = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
258eer[ell] = mean(Yres$Ypred==ydata.test)
259cat("\n", "ell =", ell)
260}
261eer
262mean(eer)
263source('~/Dropbox/projet spectres/cluster/EMLiqArtBandeHugues.R')
264library(fda)
265library(Matrix)
266library(mvtnorm)
267source('EMalgorithm.R')
268source('EMLiqArtBandeHugues.R')
269for (nVC in 1:100){
270EMLiqArtBandeHugues(nVC)
271}
272setwd("~/Dropbox/projet spectres/cluster")
273library(fda)
274library(Matrix)
275library(mvtnorm)
276source('EMalgorithm.R')
277source('EMLiqArtBandeHugues.R')
278for (nVC in 1:100){
279EMLiqArtBandeHugues(nVC)
280}
281sims= 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)){
284aa=cbind(sims,rep(9, each =length(sims)))
285colnames(aa)=c("isim","model")
286write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
287setwd("~/Dropbox/WarpMixedModel/warpingMixedModel/SuperComputer")
288write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
289sims= 1:100
290colnames(sims) = "isim"
291sims= 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)){
294aa=cbind(sims,rep(10, each =length(sims)))
295colnames(aa)=c("isim","model")
296write.table(aa,file=paste("SimulationSettings10.txt",sep=""),row.names=FALSE,sep=",")
297load("~/valse/data/data.RData")
298load("~/valse/data/data.RData")
299X
300Y
301library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
302library(roxygen2)
303setwd("~/valse")
304document()
305document()
306setwd("~/")
307install("valse")
308c(9.75,
30920,
31011.75,
31111,
3128.25,
31312.5,
31411,
31518,
3167,
31712.75,
31813,
31914.75,
3204.75,
32111,
32220,
32313.5,
3248,
32513.25,
3266,
32717.5,
32813.25,
3298.5,
3309.5,
33116,
3328,
3339.5,
33410.25,
33513.75,
3369,
33714.75,
33812.5,
33919.5,
34017.5,
3417,
34211.5,
3434,
3447.5,
34513.25,
34610.5
347)
348notes = c(9.75,
34920,
35011.75,
35111,
3528.25,
35312.5,
35411,
35518,
3567,
35712.75,
35813,
35914.75,
3604.75,
36111,
36220,
36313.5,
3648,
36513.25,
3666,
36717.5,
36813.25,
3698.5,
3709.5,
37116,
3728,
3739.5,
37410.25,
37513.75,
3769,
37714.75,
37812.5,
37919.5,
38017.5,
3817,
38211.5,
3834,
3847.5,
38513.25,
38610.5
387)
388hist(notes)
389hist(notes, nclass = 20)
390notesBis = 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,
3916.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)
392hist(notesBis)
393notes = c(7,
39419,
39510,
3969,
3976,
39812,
39911,
40019,
40111,
40212,
40316,
40413,
4052,
40613,
40718,
40815,
4098,
41016,
4119,
41215,
41312,
41411,
41510,
41614,
4177,
41811,
41912,
42016,
42113,
42215,
42314,
42419,
42514,
4264, 9,
4278,
4287,
4298,
4308
431)
432hist(notes)
433shiny::runApp('Dropbox/cranview-master')
434packages = 'shock'
435min_date = Sys.Date() - 1
436for (pkg in packages)
437{
438# api data for package. we want the initial release - the first element of the "timeline"
439pkg_data = httr::GET(paste0("http://crandb.r-pkg.org/", pkg, "/all"))
440pkg_data = httr::content(pkg_data)
441initial_release = pkg_data$timeline[[1]]
442min_date = min(min_date, as.Date(initial_release))
443}
444min_date
445library("capushe", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
446A = matrix(rnorm(25), ncol = 5)
447A = rbind(A, matrix(0, ncol = 5, nrow = 2))
448A
449A[rowSums(A)==0,]=c()
450A[rowSums(A)==0,]=[]
451A = A[rowSums(A)!=0,]
452A
453source('~/valse/pkg/R/valse.R')
454setwd("~/valse/reports")
455XY = cbind(X,Y)
456X
457p = 10
458q = 8
459k = 2
460D = 20
461meanX = rep(0,p)
462covX = 0.1*diag(p)
463covY = array(dim = c(q,q,k))
464covY[,,1] = 0.1*diag(q)
465covY[,,2] = 0.2*diag(q)
466beta = array(dim = c(p,q,2))
467beta[,,2] = matrix(c(rep(2,(D)),rep(0, p*q-D)))
468beta[,,1] = matrix(c(rep(1,D),rep(0, p*q-D)))
469n = 100
470pi = c(0.4,0.6)
471data = generateXY(meanX,covX,covY, pi, beta, n)
472source('~/valse/pkg/R/generateSampleInputs.R')
473data = generateXY(meanX,covX,covY, pi, beta, n)
474X = data$X
475Y = data$Y
476XY = cbind(X,Y)
477data
478data$class
479affec = data$class
480XY = cbind(X,Y)
481for (r in 1:K){
482XY_class[[r]] = XY[affec == r, ]
483}
484K= 2
485for (r in 1:K){
486XY_class[[r]] = XY[affec == r, ]
487}
488XY_class= list()
489for (r in 1:K){
490XY_class[[r]] = XY[affec == r, ]
491}
492for (r in 1:K){
493XY_class[[r]] = XY[affec == r, ]
494meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
495}
496meanPerClass= matrix()
497for (r in 1:K){
498XY_class[[r]] = XY[affec == r, ]
499meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
500}
501apply(XY_class[[r]], 2, mean)
502p
503q
504meanPerClass[r,] = apply(XY_class[[r]], 2, mean)
505dim(XY)
506meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2])
507for (r in 1:K){
508XY_class[[r]] = XY[affec == r, ]
509meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
510}
511matplot(meanPerClass)
512matplot(meanPerClass, type='l')