[Yuima-commits] r257 - in pkg/yuima: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Oct 27 17:47:46 CET 2013
Author: kyuta
Date: 2013-10-27 17:47:46 +0100 (Sun, 27 Oct 2013)
New Revision: 257
Added:
pkg/yuima/src/
pkg/yuima/src/cce_functions.c
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/NEWS
pkg/yuima/R/cce.R
pkg/yuima/R/llag.R
pkg/yuima/R/sim.euler.R
pkg/yuima/man/bns.test.Rd
pkg/yuima/man/cce.Rd
pkg/yuima/man/llag.Rd
pkg/yuima/man/mpv.Rd
pkg/yuima/man/noisy.sampling.Rd
Log:
add cce_functions.c
modify cce.R, llag.R, sim.euler.R, bns.test.Rd, cce.Rd, llag.Rd, mpv.Rd, noisy.sampling.Rd
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2013-10-11 13:59:30 UTC (rev 256)
+++ pkg/yuima/DESCRIPTION 2013-10-27 16:47:46 UTC (rev 257)
@@ -1,14 +1,14 @@
-Package: yuima
-Type: Package
-Title: The YUIMA Project package (unstable version)
-Version: 0.1.211
-Date: 2013-10-10
-Depends: methods, zoo, stats4, utils
-Suggests: cubature, mvtnorm
-Author: YUIMA Project Team.
-Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
-Description: The YUIMA Project for Simulation and Inference
- for Stochastic Differential Equations.
-License: GPL-2
-URL: http://R-Forge.R-project.org/projects/yuima/
-
+Package: yuima
+Type: Package
+Title: The YUIMA Project package (unstable version)
+Version: 0.1.212
+Date: 2013-10-28
+Depends: methods, zoo, stats4, utils
+Suggests: cubature, mvtnorm
+Author: YUIMA Project Team.
+Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
+Description: The YUIMA Project for Simulation and Inference
+ for Stochastic Differential Equations.
+License: GPL-2
+URL: http://R-Forge.R-project.org/projects/yuima/
+
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2013-10-11 13:59:30 UTC (rev 256)
+++ pkg/yuima/NAMESPACE 2013-10-27 16:47:46 UTC (rev 257)
@@ -1,122 +1,123 @@
-import("methods")
-
-##importFrom("stats", "end", "time", "start")
-importFrom("graphics", "plot")
-importFrom("zoo")
-importFrom(stats, confint)
-
-
-importFrom(utils, toLatex)
-
-
-
-
-
-
-exportClasses("yuima",
- "yuima.data",
- "yuima.sampling",
- "yuima.characteristic",
- "yuima.model",
- "model.parameter"
- )
-
-exportMethods(
- "dim",
- "length",
- ## "start",
- "plot",
- ## "time",
- ## "end",
- "simulate",
- "cce",
- "llag",
- "poisson.random.sampling",
- "get.zoo.data",
- "initialize",
-## "ql",
-## "rql",
-## "ml.ql",
- "adaBayes",
- "limiting.gamma",
- "getF",
- "getf",
- "getxinit",
- "gete",
- "simFunctional",
- "F0",
- "Fnorm",
- "asymptotic_term",
- "cbind.yuima"
- )
-
-## function which we want to expose to the user
-## all other functions are used internally by the
-## package
-export(setYuima)
-export(setModel) ## builds sde model
-export(setData)
-export(setSampling)
-export(setCharacteristic)
-
-export(dim)
-export(length)
-#export(start)
-export(plot)
-#export(time)
-#export(end)
-
-export(simulate) # simulates couple of processes
-export(subsampling)
-export(cce)
-export(llag)
-export(poisson.random.sampling)
-export(noisy.sampling)
-export(mpv)
-export(bns.test)
-
-export(get.zoo.data)
-
-##export(ql,rql,ml.ql)
-##export(rql)
-export(adaBayes)
-export(rIG, rNIG, rbgamma, rngamma, rstable) ##:: random number generator for Inverse Gaussian
-export(limiting.gamma)
-
-export(setFunctional)
-export(getF)
-export(getf)
-export(getxinit)
-export(gete)
-
-export(simFunctional)
-export(F0)
-export(Fnorm)
-export(asymptotic_term)
-
-##export(LSE)
-export(lse)
-
-export(qmle)
-export(quasilogl)
-export(phi.test)
-export(lasso)
-export(CPoint)
-export(qmleR)
-export(qmleL)
-
-export(qgv)
-export(mmfrac)
-
-
-
-export(cbind.yuima)
-
-S3method(print, phitest)
-S3method(print, qgv)
-S3method(print, mmfrac)
-
-S3method(toLatex, yuima)
-S3method(toLatex, yuima.model)
-
-
+import("methods")
+
+##importFrom("stats", "end", "time", "start")
+importFrom("graphics", "plot")
+importFrom("zoo")
+importFrom(stats, confint)
+
+
+importFrom(utils, toLatex)
+
+
+
+
+
+
+exportClasses("yuima",
+ "yuima.data",
+ "yuima.sampling",
+ "yuima.characteristic",
+ "yuima.model",
+ "model.parameter"
+ )
+
+exportMethods(
+ "dim",
+ "length",
+ ## "start",
+ "plot",
+ ## "time",
+ ## "end",
+ "simulate",
+ "cce",
+ "llag",
+ "poisson.random.sampling",
+ "get.zoo.data",
+ "initialize",
+## "ql",
+## "rql",
+## "ml.ql",
+ "adaBayes",
+ "limiting.gamma",
+ "getF",
+ "getf",
+ "getxinit",
+ "gete",
+ "simFunctional",
+ "F0",
+ "Fnorm",
+ "asymptotic_term",
+ "cbind.yuima"
+ )
+
+## function which we want to expose to the user
+## all other functions are used internally by the
+## package
+export(setYuima)
+export(setModel) ## builds sde model
+export(setData)
+export(setSampling)
+export(setCharacteristic)
+
+export(dim)
+export(length)
+#export(start)
+export(plot)
+#export(time)
+#export(end)
+
+export(simulate) # simulates couple of processes
+export(subsampling)
+export(cce)
+export(llag)
+export(poisson.random.sampling)
+export(noisy.sampling)
+export(mpv)
+export(bns.test)
+
+export(get.zoo.data)
+
+##export(ql,rql,ml.ql)
+##export(rql)
+export(adaBayes)
+export(rIG, rNIG, rbgamma, rngamma, rstable) ##:: random number generator for Inverse Gaussian
+export(limiting.gamma)
+
+export(setFunctional)
+export(getF)
+export(getf)
+export(getxinit)
+export(gete)
+
+export(simFunctional)
+export(F0)
+export(Fnorm)
+export(asymptotic_term)
+
+##export(LSE)
+export(lse)
+
+export(qmle)
+export(quasilogl)
+export(phi.test)
+export(lasso)
+export(CPoint)
+export(qmleR)
+export(qmleL)
+
+export(qgv)
+export(mmfrac)
+
+
+
+export(cbind.yuima)
+
+S3method(print, phitest)
+S3method(print, qgv)
+S3method(print, mmfrac)
+
+S3method(toLatex, yuima)
+S3method(toLatex, yuima.model)
+
+useDynLib(yuima)
+
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2013-10-11 13:59:30 UTC (rev 256)
+++ pkg/yuima/NEWS 2013-10-27 16:47:46 UTC (rev 257)
@@ -11,3 +11,5 @@
2013/04/14: modify qmle.R
2013/04/14: modify adaBayes.R
2013/04/14: modify bns.test.Rd, mpv.Rd
+2013/10/28: add cce_functions.c
+ modify cce.R, llag.R, sim.euler.R, bns.test.Rd, cce.Rd, llag.Rd, mpv.Rd, noisy.sampling.Rd
\ No newline at end of file
Modified: pkg/yuima/R/cce.R
===================================================================
--- pkg/yuima/R/cce.R 2013-10-11 13:59:30 UTC (rev 256)
+++ pkg/yuima/R/cce.R 2013-10-27 16:47:46 UTC (rev 257)
@@ -25,64 +25,82 @@
ser.times <- lapply(lapply(data,"time"),"as.numeric")
ser.lengths <- sapply(data,"length")
ser.samplings <- vector(d.size, mode="list")
- refresh.times <- c()
+ #refresh.times <- c()
- tmp <- sapply(ser.times,"[",1)
- refresh.times[1] <- max(tmp)
-
- I <- rep(1,d.size)
-
- for(d in 1:d.size){
- #ser.samplings[[d]][1] <- max(which(ser.times[[d]]<=refresh.times[1]))
- while(!(ser.times[[d]][I[d]+1]>refresh.times[1])){
- I[d] <- I[d]+1
- if((I[d]+1)>ser.lengths[d]){
- break
- }
- }
- ser.samplings[[d]][1] <- I[d]
- }
-
- i <- 1
-
- #while(all(sapply(ser.samplings,"[",i)<ser.lengths)){
- while(all(I<ser.lengths)){
+ if(0){
+ tmp <- sapply(ser.times,"[",1)
+ refresh.times[1] <- max(tmp)
- J <- I
- tmp <- double(d.size)
+ I <- rep(1,d.size)
for(d in 1:d.size){
- #tmp[d] <- ser.times[[d]][min(which(ser.times[[d]]>refresh.times[i]))
- repeat{
- J[d] <- J[d] + 1
- tmp[d] <- ser.times[[d]][J[d]]
- if((!(J[d]<ser.lengths))||(tmp[d]>refresh.times[i])){
+ #ser.samplings[[d]][1] <- max(which(ser.times[[d]]<=refresh.times[1]))
+ while(!(ser.times[[d]][I[d]+1]>refresh.times[1])){
+ I[d] <- I[d]+1
+ if((I[d]+1)>ser.lengths[d]){
break
}
}
+ ser.samplings[[d]][1] <- I[d]
}
- i <- i+1
+ i <- 1
- refresh.times[i] <- max(tmp)
-
- for(d in 1:d.size){
- #ser.samplings[[d]][i] <- max(which(ser.times[[d]]<=refresh.times[i]))
- while(!(ser.times[[d]][I[d]+1]>refresh.times[i])){
- I[d] <- I[d]+1
- if((I[d]+1)>ser.lengths[d]){
- break
+ #while(all(sapply(ser.samplings,"[",i)<ser.lengths)){
+ while(all(I<ser.lengths)){
+
+ J <- I
+ tmp <- double(d.size)
+
+ for(d in 1:d.size){
+ #tmp[d] <- ser.times[[d]][min(which(ser.times[[d]]>refresh.times[i]))
+ repeat{
+ J[d] <- J[d] + 1
+ tmp[d] <- ser.times[[d]][J[d]]
+ if((!(J[d]<ser.lengths))||(tmp[d]>refresh.times[i])){
+ break
+ }
}
}
- ser.samplings[[d]][i] <- I[d]
+
+ i <- i+1
+
+ refresh.times[i] <- max(tmp)
+
+ for(d in 1:d.size){
+ #ser.samplings[[d]][i] <- max(which(ser.times[[d]]<=refresh.times[i]))
+ while(!(ser.times[[d]][I[d]+1]>refresh.times[i])){
+ I[d] <- I[d]+1
+ if((I[d]+1)>ser.lengths[d]){
+ break
+ }
+ }
+ ser.samplings[[d]][i] <- I[d]
+ }
+
}
-
}
+ MinL <- min(ser.lengths)
+
+ refresh.times <- double(MinL)
+ refresh.times[1] <- max(sapply(ser.times,"[",1))
+
+ idx <- matrix(.C("refreshsampling",
+ as.integer(d.size),
+ integer(d.size),
+ as.double(unlist(ser.times)),
+ as.double(refresh.times),
+ as.integer(append(ser.lengths,0,after=0)),
+ min(sapply(ser.times,FUN="tail",n=1)),
+ as.integer(MinL),
+ result=integer(d.size*MinL))$result,ncol=d.size)
+
result <- vector(d.size, mode="list")
for(d in 1:d.size){
- result[[d]] <- data[[d]][ser.samplings[[d]]]
+ #result[[d]] <- data[[d]][ser.samplings[[d]]]
+ result[[d]] <- data[[d]][idx[,d]]
}
return(result)
@@ -104,48 +122,70 @@
ser.times <- lapply(lapply(data,"time"),"as.numeric")
ser.lengths <- sapply(data,"length")
- refresh.times <- max(sapply(ser.times,"[",1))
+ #refresh.times <- max(sapply(ser.times,"[",1))
ser.samplings <- vector(d.size,mode="list")
- for(d in 1:d.size){
- ser.samplings[[d]][1] <- 1
- }
-
- I <- rep(1,d.size)
- i <- 1
-
- while(all(I<ser.lengths)){
-
- flag <- integer(d.size)
-
+ if(0){
for(d in 1:d.size){
- while(I[d]<ser.lengths[d]){
- I[d] <- I[d]+1
- if(ser.times[[d]][I[d]]>refresh.times[i]){
- flag[d] <- 1
- ser.samplings[[d]][i+1] <- I[d]
- break
- }
- }
+ ser.samplings[[d]][1] <- 1
}
- if(any(flag<rep(1,d.size))){
- break
- }else{
- i <- i+1
- tmp <- double(d.size)
+ I <- rep(1,d.size)
+ i <- 1
+
+ while(all(I<ser.lengths)){
+
+ flag <- integer(d.size)
+
for(d in 1:d.size){
- tmp[d] <- ser.times[[d]][ser.samplings[[d]][i]]
+ while(I[d]<ser.lengths[d]){
+ I[d] <- I[d]+1
+ if(ser.times[[d]][I[d]]>refresh.times[i]){
+ flag[d] <- 1
+ ser.samplings[[d]][i+1] <- I[d]
+ break
+ }
+ }
}
- refresh.times <- append(refresh.times,max(tmp))
+
+ if(any(flag<rep(1,d.size))){
+ break
+ }else{
+ i <- i+1
+ tmp <- double(d.size)
+ for(d in 1:d.size){
+ tmp[d] <- ser.times[[d]][ser.samplings[[d]][i]]
+ }
+ refresh.times <- append(refresh.times,max(tmp))
+ }
+
}
-
}
+ MinL <- min(ser.lengths)
+
+ refresh.times <- double(MinL)
+ refresh.times[1] <- max(sapply(ser.times,"[",1))
+
+ obj <- .C("refreshsamplingphy",
+ as.integer(d.size),
+ integer(d.size),
+ as.double(unlist(ser.times)),
+ rtimes=as.double(refresh.times),
+ as.integer(append(ser.lengths,0,after=0)),
+ min(sapply(ser.times,FUN="tail",n=1)),
+ as.integer(MinL),
+ Samplings=integer(d.size*(MinL+1)),
+ rNum=integer(1))
+
+ refresh.times <- obj$rtimes[1:obj$rNum]
+ idx <- matrix(obj$Samplings,ncol=d.size)
+
result <- vector(d.size, mode="list")
for(d in 1:d.size){
- result[[d]] <- data[[d]][ser.samplings[[d]]]
+ #result[[d]] <- data[[d]][ser.samplings[[d]]]
+ result[[d]] <- data[[d]][idx[,d]]
}
return(list(data=result,refresh.times=refresh.times))
@@ -166,116 +206,147 @@
xlength <- length(xtime)
ylength <- length(ytime)
+ N.max <- max(xlength,ylength)
- mu <- c()
- w <- c()
- q <- c()
- r <- c()
-
- if(xtime[1]<ytime[1]){
- I <- 1
- while(xtime[I]<ytime[1]){
- I <- I+1
- if(!(I<xlength)){
- break
- }
- }
- #mu0 <- min(which(ytime[1]<=xtime))
- mu0 <- I
- if(ytime[1]<xtime[mu0]){
- q[1] <- mu0-1
- }else{
- q[1] <- mu0
- }
- r[1] <- 1
- }else if(xtime[1]>ytime[1]){
- I <- 1
- while(xtime[I]<ytime[1]){
- I <- I+1
- if(!(I<xlength)){
- break
- }
- }
- #w0 <- min(which(xtime[1]<=ytime))
- w0 <- I
- q[1] <- 1
- if(xtime[1]<ytime[w0]){
- r[1] <- w0-1
- }else{
- r[1] <- w0
- }
- }else{
- q[1] <- 1
- r[1] <- 1
- }
-
- i <- 1
-
- #repeat{
- while((q[i]<xlength)&&(r[i]<ylength)){
- if(xtime[q[i]]<ytime[r[i]]){
- #tmp <- which(ytime[r[i]]<=xtime)
- #if(identical(all.equal(tmp,integer(0)),TRUE)){
- # break
- #}
- #mu[i] <- min(tmp)
- if(xtime[xlength]<ytime[r[i]]){
- break
- }
- I <- q[i]
- while(xtime[I]<ytime[r[i]]){
+ if(0){
+ mu <- integer(N.max)
+ w <- integer(N.max)
+ q <- integer(N.max)
+ r <- integer(N.max)
+
+ if(xtime[1]<ytime[1]){
+ I <- 1
+ while(xtime[I]<ytime[1]){
I <- I+1
+ if(!(I<xlength)){
+ break
+ }
}
- mu[i] <- I
- w[i] <- r[i]
- if(ytime[r[i]]<xtime[mu[i]]){
- q[i+1] <- mu[i]-1
+ #mu0 <- min(which(ytime[1]<=xtime))
+ mu0 <- I
+ if(ytime[1]<xtime[mu0]){
+ #q[1] <- mu0-1
+ q[1] <- mu0
}else{
- q[i+1] <- mu[i]
+ #q[1] <- mu0
+ q[1] <- mu0+1
}
- r[i+1] <- r[i]+1
- }else if(xtime[q[i]]>ytime[r[i]]){
- #tmp <- which(xtime[q[i]]<=ytime)
- #if(identical(all.equal(tmp,integer(0)),TRUE)){
- # break
- #}
- if(xtime[q[i]]>ytime[ylength]){
- break
- }
- mu[i] <- q[i]
- #w[i] <- min(tmp)
- I <- r[i]
- while(xtime[q[i]]>ytime[I]){
+ r[1] <- 2
+ }else if(xtime[1]>ytime[1]){
+ I <- 1
+ while(xtime[I]<ytime[1]){
I <- I+1
+ if(!(I<xlength)){
+ break
+ }
}
- w[i] <- I
- q[i+1] <- q[i]+1
- if(xtime[q[i]]<ytime[w[i]]){
- r[i+1] <- w[i]-1
+ #w0 <- min(which(xtime[1]<=ytime))
+ w0 <- I
+ q[1] <- 2
+ if(xtime[1]<ytime[w0]){
+ #r[1] <- w0-1
+ r[1] <- w0
}else{
- r[i+1] <- w[i]
+ #r[1] <- w0
+ r[1] <- w0+1
}
}else{
- mu[i] <- q[i]
- w[i] <- r[i]
- q[i+1] <- q[i]+1
- r[i+1] <- r[i]+1
+ q[1] <- 2
+ r[1] <- 2
}
- i <- i+1
+ i <- 1
- #if((q[i]>xlength)||(r[i]>ylength)){
- # break
- #}
+ repeat{
+ #while((q[i]<xlength)&&(r[i]<ylength)){
+ if(xtime[q[i]]<ytime[r[i]]){
+ #tmp <- which(ytime[r[i]]<=xtime)
+ #if(identical(all.equal(tmp,integer(0)),TRUE)){
+ # break
+ #}
+ #mu[i] <- min(tmp)
+ if(xtime[xlength]<ytime[r[i]]){
+ break
+ }
+ I <- q[i]
+ while(xtime[I]<ytime[r[i]]){
+ I <- I+1
+ }
+ mu[i] <- I
+ w[i] <- r[i]
+ if(ytime[r[i]]<xtime[mu[i]]){
+ #q[i+1] <- mu[i]-1
+ q[i+1] <- mu[i]
+ }else{
+ #q[i+1] <- mu[i]
+ q[i+1] <- mu[i]+1
+ }
+ r[i+1] <- r[i]+1
+ }else if(xtime[q[i]]>ytime[r[i]]){
+ #tmp <- which(xtime[q[i]]<=ytime)
+ #if(identical(all.equal(tmp,integer(0)),TRUE)){
+ # break
+ #}
+ if(xtime[q[i]]>ytime[ylength]){
+ break
+ }
+ mu[i] <- q[i]
+ #w[i] <- min(tmp)
+ I <- r[i]
+ while(xtime[q[i]]>ytime[I]){
+ I <- I+1
+ }
+ w[i] <- I
+ q[i+1] <- q[i]+1
+ if(xtime[q[i]]<ytime[w[i]]){
+ #r[i+1] <- w[i]-1
+ r[i+1] <- w[i]
+ }else{
+ #r[i+1] <- w[i]
+ r[i+1] <- w[i]+1
+ }
+ }else{
+ mu[i] <- q[i]
+ w[i] <- r[i]
+ q[i+1] <- q[i]+1
+ r[i+1] <- r[i]+1
+ }
+
+ i <- i+1
+
+ if((q[i]>=xlength)||(r[i]>=ylength)){
+ break
+ }
+
+ }
+ mu[i] <- q[i]
+ w[i] <- r[i]
}
- return(list(xg=as.numeric(x)[mu[1:(i-1)]],
- xl=as.numeric(x)[q[1:(i-1)]],
- ygamma=as.numeric(y)[w[1:(i-1)]],
- ylambda=as.numeric(y)[r[1:(i-1)]],
- num.data=i-1))
+ sdata <- .C("bibsynchro",
+ as.double(xtime),
+ as.double(ytime),
+ as.integer(xlength),
+ as.integer(ylength),
+ mu=integer(N.max),
+ w=integer(N.max),
+ q=integer(N.max),
+ r=integer(N.max),
+ Num=integer(1))
+ Num <- sdata$Num
+
+ #return(list(xg=as.numeric(x)[mu[1:i]],
+ # xl=as.numeric(x)[q[1:i]-1],
+ # ygamma=as.numeric(y)[w[1:i]],
+ # ylambda=as.numeric(y)[r[1:i]-1],
+ # num.data=i))
+ return(list(xg=as.numeric(x)[sdata$mu[1:Num]+1],
+ xl=as.numeric(x)[sdata$q[1:Num]],
+ ygamma=as.numeric(y)[sdata$w[1:Num]+1],
+ ylambda=as.numeric(y)[sdata$r[1:Num]],
+ num.data=Num))
}
@@ -295,28 +366,41 @@
grid <- seq(ztime[1],end,by=frequency)
n.sparse <- length(grid)
- result <- double(frequency)
+ #result <- double(frequency)
#result <- 0
- I <- rep(1,n.sparse)
+ #I <- rep(1,n.sparse)
- for(t in 1:frequency){
- for(i in 1:n.sparse){
- while((ztime[I[i]+1]<=grid[i])&&(I[i]<n.size)){
- I[i] <- I[i]+1
- }
- }
- result[t] <- sum(diff(znum[I])^2)
+ #for(t in 1:frequency){
+ # for(i in 1:n.sparse){
+ # while((ztime[I[i]+1]<=grid[i])&&(I[i]<n.size)){
+ # I[i] <- I[i]+1
+ # }
+ # }
+ # result[t] <- sum(diff(znum[I])^2)
#result <- result+sum(diff(znum[I])^2)
- grid <- grid+rep(1,n.sparse)
- if(grid[n.sparse]>end){
- grid <- grid[-n.sparse]
- I <- I[-n.sparse]
- n.sparse <- n.sparse-1
- }
- }
+ # grid <- grid+rep(1,n.sparse)
+ # if(grid[n.sparse]>end){
+ # grid <- grid[-n.sparse]
+ # I <- I[-n.sparse]
+ # n.sparse <- n.sparse-1
+ # }
+ #}
+ K <- floor(end-grid[n.sparse]) + 1
+
+ zmat <- matrix(.C("ctsubsampling",as.double(znum),as.double(ztime),
+ as.integer(frequency),as.integer(n.sparse),
+ as.integer(n.size),as.double(grid),
+ result=double(frequency*n.sparse))$result,
+ n.sparse,frequency)
+
+ result <- double(frequency)
+ result[1:K] <- colSums(diff(zmat[,1:K])^2)
+ result[-(1:K)] <- colSums(diff(zmat[-n.sparse,-(1:K)])^2)
+
return(mean(result))
#return(result/frequency)
+ #return(znum[I])
}
Omega_BNHLS <- function(zdata,sec=120,utime){
@@ -530,33 +614,37 @@
cmat <- matrix(0, n.series, n.series) # cov
for(i in 1:n.series){
for(j in i:n.series){
- I <- rep(1,n.series)
+ #I <- rep(1,n.series)
#Checking Starting Point
#repeat{
- while((I[i]<length(ser.times[[i]])) && (I[j]<length(ser.times[[j]]))){
- if(ser.times[[i]][I[i]] >= ser.times[[j]][I[j]+1]){
- I[j] <- I[j]+1
- }else if(ser.times[[i]][I[i]+1] <= ser.times[[j]][I[j]]){
- I[i] <- I[i]+1
- }else{
- break
- }
- }
+ #while((I[i]<length(ser.times[[i]])) && (I[j]<length(ser.times[[j]]))){
+ # if(ser.times[[i]][I[i]] >= ser.times[[j]][I[j]+1]){
+ # I[j] <- I[j]+1
+ # }else if(ser.times[[i]][I[i]+1] <= ser.times[[j]][I[j]]){
+ # I[i] <- I[i]+1
+ # }else{
+ # break
+ # }
+ #}
#Main Component
if(i!=j){
- while((I[i]<length(ser.times[[i]])) && (I[j]<length(ser.times[[j]]))) {
- cmat[j,i] <- cmat[j,i] + (ser.diffX[[i]])[I[i]]*(ser.diffX[[j]])[I[j]]
- if(ser.times[[i]][I[i]+1]>ser.times[[j]][I[j]+1]){
- I[j] <- I[j] + 1
- }else if(ser.times[[i]][I[i]+1]<ser.times[[j]][I[j]+1]){
- I[i] <- I[i] +1
- }else{
- I[i] <- I[i]+1
- I[j] <- I[j]+1
- }
- }
+ #while((I[i]<length(ser.times[[i]])) && (I[j]<length(ser.times[[j]]))) {
+ # cmat[j,i] <- cmat[j,i] + (ser.diffX[[i]])[I[i]]*(ser.diffX[[j]])[I[j]]
+ # if(ser.times[[i]][I[i]+1]>ser.times[[j]][I[j]+1]){
+ # I[j] <- I[j] + 1
+ # }else if(ser.times[[i]][I[i]+1]<ser.times[[j]][I[j]+1]){
+ # I[i] <- I[i] +1
+ # }else{
+ # I[i] <- I[i]+1
+ # I[j] <- I[j]+1
+ # }
+ #}
+ cmat[j,i] <- .C("HayashiYoshida",as.integer(length(ser.times[[i]])),
+ as.integer(length(ser.times[[j]])),as.double(ser.times[[i]]),
+ as.double(ser.times[[j]]),as.double(ser.diffX[[i]]),
+ as.double(ser.diffX[[j]]),value=double(1))$value
}else{
cmat[i,j] <- sum(ser.diffX[[i]]^2)
}
@@ -588,10 +676,12 @@
cmat <- matrix(0, n.series, n.series) # cov
# synchronize the data and take the differences
- diffX <- lapply(lapply(refresh_sampling(data),"as.numeric"),"diff")
- diffX <- do.call("rbind",diffX) # transform to matrix
+ #diffX <- lapply(lapply(refresh_sampling(data),"as.numeric"),"diff")
+ #diffX <- do.call("rbind",diffX) # transform to matrix
+ diffX <- diff(do.call("cbind",lapply(refresh_sampling(data),"as.numeric")))
- Num <- ncol(diffX)
+ #Num <- ncol(diffX)
+ Num <- nrow(diffX)
if(missing(kn)) kn <- max(floor(mean(theta)*Num^(1/2+delta)),2)
@@ -599,15 +689,17 @@
psi2.kn <- sum(weight^2)
# pre-averaging
- myfunc <- function(dx)rollapplyr(dx,width=kn-1,FUN="%*%",weight)
- barX <- apply(diffX,1,FUN=myfunc)
+ #myfunc <- function(dx)rollapplyr(dx,width=kn-1,FUN="%*%",weight)
+ #barX <- apply(diffX,1,FUN=myfunc)
+ barX <- filter(diffX,filter=rev(weight),method="c",sides=1)[(kn-1):Num,]
cmat <- (Num/(Num-kn+2))*t(barX)%*%barX/psi2.kn
if(delta==0){
psi1.kn <- kn^2*sum(diff(c(0,weight,0))^2)
scale <- psi1.kn/(theta^2*psi2.kn*2*Num)
- cmat <- cmat-scale*diffX%*%t(diffX)
+ #cmat <- cmat-scale*diffX%*%t(diffX)
+ cmat <- cmat-scale*t(diffX)%*%diffX
if(adj) cmat <- cmat/(1-scale)
}
@@ -661,23 +753,36 @@
weight <- sapply((1:(kn-1))/kn,g)
psi.kn <- sum(weight)
- ser.barX <- list(rollapplyr(ser.diffX[[1]],width=kn-1,FUN="%*%",weight),
- rollapplyr(ser.diffX[[2]],width=kn-1,FUN="%*%",weight))
+ #ser.barX <- list(rollapplyr(ser.diffX[[1]],width=kn-1,FUN="%*%",weight),
+ # rollapplyr(ser.diffX[[2]],width=kn-1,FUN="%*%",weight))
+ ser.barX <- list(filter(ser.diffX[[1]],rev(weight),method="c",
+ sides=1)[(kn-1):length(ser.diffX[[1]])],
+ filter(ser.diffX[[2]],rev(weight),method="c",
+ sides=1)[(kn-1):length(ser.diffX[[2]])])
ser.num.barX <- sapply(ser.barX,"length")-1
# core part of cce
- start <- kn+1
- end <- 1
- for(k in 1:ser.num.barX[1]){
- while(!(ser.times[[1]][k]<ser.times[[2]][start])&&((start-kn)<ser.num.barX[2])){
- start <- start + 1
- }
- while((ser.times[[1]][k+kn]>ser.times[[2]][end+1])&&(end<ser.num.barX[2])){
- end <- end + 1
- }
- cmat[i,j] <- cmat[i,j] + ser.barX[[1]][k]*sum(ser.barX[[2]][(start-kn):end])
- }
+ #start <- kn+1
+ #end <- 1
+ #for(k in 1:ser.num.barX[1]){
+ # while(!(ser.times[[1]][k]<ser.times[[2]][start])&&((start-kn)<ser.num.barX[2])){
+ # start <- start + 1
+ # }
+ # while((ser.times[[1]][k+kn]>ser.times[[2]][end+1])&&(end<ser.num.barX[2])){
+ # end <- end + 1
+ # }
+ # cmat[i,j] <- cmat[i,j] + ser.barX[[1]][k]*sum(ser.barX[[2]][(start-kn):end])
+ #}
+ cmat[i,j] <- .C("pHayashiYoshida",
+ as.integer(kn),
+ as.integer(ser.num.barX[1]),
+ as.integer(ser.num.barX[2]),
+ as.double(ser.times[[1]]),
+ as.double(ser.times[[2]]),
+ as.double(ser.barX[[1]]),
+ as.double(ser.barX[[2]]),
+ value=double(1))$value
cmat[i,j] <- cmat[i,j]/(psi.kn^2)
cmat[j,i] <- cmat[i,j]
@@ -691,9 +796,13 @@
weight <- sapply((1:(kn-1))/kn,g)
psi.kn <- sum(weight)
- barX <- rollapplyr(diffX,width=kn-1,FUN="%*%",weight)
+ #barX <- rollapplyr(diffX,width=kn-1,FUN="%*%",weight)
+ barX <- filter(diffX,rev(weight),method="c",
+ sides=1)[(kn-1):length(diffX)]
tmp <- barX[-length(barX)]
- cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kn-1)+1,FUN="sum",partial=TRUE)/(psi.kn)^2
+ #cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kn-1)+1,FUN="sum",partial=TRUE)/(psi.kn)^2
+ cmat[i,j] <- tmp%*%rollsum(c(double(kn-1),tmp,double(kn-1)),
+ k=2*(kn-1)+1)/(psi.kn)^2
}
}
@@ -723,23 +832,36 @@
weight <- sapply((1:(knij-1))/knij,g)
psi.kn <- sum(weight)
- ser.barX <- list(rollapplyr(ser.diffX[[1]],width=knij-1,FUN="%*%",weight),
- rollapplyr(ser.diffX[[2]],width=knij-1,FUN="%*%",weight))
+ #ser.barX <- list(rollapplyr(ser.diffX[[1]],width=knij-1,FUN="%*%",weight),
+ # rollapplyr(ser.diffX[[2]],width=knij-1,FUN="%*%",weight))
+ ser.barX <- list(filter(ser.diffX[[1]],rev(weight),method="c",
+ sides=1)[(knij-1):length(ser.diffX[[1]])],
+ filter(ser.diffX[[2]],rev(weight),method="c",
+ sides=1)[(knij-1):length(ser.diffX[[2]])])
ser.num.barX <- sapply(ser.barX,"length")-1
# core part of cce
- start <- knij+1
- end <- 1
- for(k in 1:ser.num.barX[1]){
- while(!(ser.times[[1]][k]<ser.times[[2]][start])&&((start-knij)<ser.num.barX[2])){
- start <- start + 1
- }
- while((ser.times[[1]][k+knij]>ser.times[[2]][end+1])&&(end<ser.num.barX[2])){
- end <- end + 1
- }
- cmat[i,j] <- cmat[i,j] + ser.barX[[1]][k]*sum(ser.barX[[2]][(start-knij):end])
- }
+ #start <- knij+1
+ #end <- 1
+ #for(k in 1:ser.num.barX[1]){
+ # while(!(ser.times[[1]][k]<ser.times[[2]][start])&&((start-knij)<ser.num.barX[2])){
+ # start <- start + 1
+ # }
+ # while((ser.times[[1]][k+knij]>ser.times[[2]][end+1])&&(end<ser.num.barX[2])){
+ # end <- end + 1
+ # }
+ # cmat[i,j] <- cmat[i,j] + ser.barX[[1]][k]*sum(ser.barX[[2]][(start-knij):end])
+ #}
+ cmat[i,j] <- .C("pHayashiYoshida",
+ as.integer(knij),
+ as.integer(ser.num.barX[1]),
+ as.integer(ser.num.barX[2]),
+ as.double(ser.times[[1]]),
+ as.double(ser.times[[2]]),
+ as.double(ser.barX[[1]]),
+ as.double(ser.barX[[2]]),
+ value=double(1))$value
cmat[i,j] <- cmat[i,j]/(psi.kn^2)
cmat[j,i] <- cmat[i,j]
@@ -752,12 +874,15 @@
weight <- sapply((1:(kni-1))/kni,g)
psi.kn <- sum(weight)
- barX <- rollapplyr(diffX,width=kni-1,FUN="%*%",weight)
+ #barX <- rollapplyr(diffX,width=kni-1,FUN="%*%",weight)
+ barX <- filter(diffX,rev(weight),method="c",
+ sides=1)[(kni-1):length(diffX)]
tmp <- barX[-length(barX)]
# core part of cce
- cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kni-1)+1,FUN="sum",partial=TRUE)/(psi.kn)^2
-
+ #cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kni-1)+1,FUN="sum",partial=TRUE)/(psi.kn)^2
+ cmat[i,j] <- tmp%*%rollsum(c(double(kni-1),tmp,double(kni-1)),
+ k=2*(kni-1)+1)/(psi.kn)^2
}
}
}
@@ -793,7 +918,9 @@
ser.diffX[[i]] <- diff( ser.X[[i]] )
# pre-averaging
- ser.barX[[i]] <- rollapplyr(ser.diffX[[i]],width=kn-1,FUN="%*%",weight)
+ #ser.barX[[i]] <- rollapplyr(ser.diffX[[i]],width=kn-1,FUN="%*%",weight)
+ ser.barX[[i]] <- filter(ser.diffX[[i]],rev(weight),method="c",
+ sides=1)[(kn-1):length(ser.diffX[[i]])]
ser.num.barX[i] <- length(ser.barX[[i]])-1
}
@@ -802,21 +929,32 @@
for(i in 1:n.series){
for(j in i:n.series){
if(i!=j){
- start <- kn+1
- end <- 1
- for(k in 1:ser.num.barX[i]){
- while(!(ser.times[[i]][k]<ser.times[[j]][start])&&((start-kn)<ser.num.barX[j])){
- start <- start + 1
- }
- while((ser.times[[i]][k+kn]>ser.times[[j]][end+1])&&(end<ser.num.barX[j])){
- end <- end + 1
- }
- cmat[i,j] <- cmat[i,j] + ser.barX[[i]][k]*sum(ser.barX[[j]][(start-kn):end])
- }
+ #start <- kn+1
+ #end <- 1
+ #for(k in 1:ser.num.barX[i]){
+ # while(!(ser.times[[i]][k]<ser.times[[j]][start])&&((start-kn)<ser.num.barX[j])){
+ # start <- start + 1
+ # }
+ # while((ser.times[[i]][k+kn]>ser.times[[j]][end+1])&&(end<ser.num.barX[j])){
+ # end <- end + 1
+ # }
+ # cmat[i,j] <- cmat[i,j] + ser.barX[[i]][k]*sum(ser.barX[[j]][(start-kn):end])
+ #}
+ cmat[i,j] <- .C("pHayashiYoshida",
+ as.integer(kn),
+ as.integer(ser.num.barX[i]),
+ as.integer(ser.num.barX[j]),
+ as.double(ser.times[[i]]),
+ as.double(ser.times[[j]]),
+ as.double(ser.barX[[i]]),
+ as.double(ser.barX[[j]]),
+ value=double(1))$value
cmat[j,i] <- cmat[i,j]
}else{
tmp <- ser.barX[[i]][1:ser.num.barX[i]]
- cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kn-1)+1,FUN="sum",partial=TRUE)
+ #cmat[i,j] <- tmp%*%rollapply(tmp,width=2*(kn-1)+1,FUN="sum",partial=TRUE)
+ cmat[i,j] <- tmp%*%rollsum(c(double(kn-1),tmp,double(kn-1)),
+ k=2*(kn-1)+1)
}
}
}
@@ -856,8 +994,12 @@
# pre-averaging
knij <- min(kn[i,j],sapply(ser.X,"length")) # kn must be less than the numbers of the observations
weight <- sapply((1:(knij-1))/knij,g)
- ser.barX <- list(rollapplyr(ser.diffX[[1]],width=knij-1,FUN="%*%",weight),
- rollapplyr(ser.diffX[[2]],width=knij-1,FUN="%*%",weight))
+ #ser.barX <- list(rollapplyr(ser.diffX[[1]],width=knij-1,FUN="%*%",weight),
+ # rollapplyr(ser.diffX[[2]],width=knij-1,FUN="%*%",weight))
+ ser.barX <- list(filter(ser.diffX[[1]],rev(weight),method="c",
+ sides=1)[(knij-1):length(ser.diffX[[1]])],
+ filter(ser.diffX[[2]],rev(weight),method="c",
+ sides=1)[(knij-1):length(ser.diffX[[2]])])
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 257
More information about the Yuima-commits
mailing list