[Yuima-commits] r155 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 3 13:18:21 CEST 2011
Author: hinohide
Date: 2011-06-03 13:18:21 +0200 (Fri, 03 Jun 2011)
New Revision: 155
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/asymptotic_term.R
Log:
bugs in asymptotic_term is fixed
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2011-05-19 12:11:48 UTC (rev 154)
+++ pkg/yuima/DESCRIPTION 2011-06-03 11:18:21 UTC (rev 155)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.190
-Date: 2011-05-19
+Version: 0.1.191
+Date: 2011-06-03
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/R/asymptotic_term.R
===================================================================
--- pkg/yuima/R/asymptotic_term.R 2011-05-19 12:11:48 UTC (rev 154)
+++ pkg/yuima/R/asymptotic_term.R 2011-06-03 11:18:21 UTC (rev 155)
@@ -1,1531 +1,1712 @@
-# in this source we note formulae like latex
-
-
-setGeneric("asymptotic_term",
- function(yuima,block=100, rho, g, expand.var="e")
- standardGeneric("asymptotic_term")
- )
-
-setMethod("asymptotic_term",signature(yuima="yuima"), function(yuima,block=100, rho, g, expand.var="e"){
-
- if(is.null(yuima at model)) stop("model object is missing!")
- if(is.null(yuima at sampling)) stop("sampling object is missing!")
- if(is.null(yuima at functional)) stop("functional object is missing!")
-
- f <- getf(yuima at functional)
- F <- getF(yuima at functional)
-
- ##:: fix bug 07/23
- #e <- gete(yuima at functional)
- assign(expand.var, gete(yuima at functional))
-
-## Terminal <- yuima at sampling@Terminal
- Terminal <- yuima at sampling@Terminal[1]
-## division <- yuima at sampling@n
- division <- yuima at sampling@n[1]
- xinit <- getxinit(yuima at functional)
- state <- yuima at model@state.variable
- V0 <- yuima at model@drift
- V <- yuima at model@diffusion
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- k.size <- length(F)
-
- print("compute X.t0")
- X.t0 <- Get.X.t0(yuima, expand.var=expand.var)
- delta <- deltat(X.t0)
-
- ##:: fix bug 07/23
- pars <- expand.var #yuima at model@parameter at all[1] #epsilon
-
- # function to return symbolic derivatives of myfunc by mystate(multi-state)
- Derivation.vector <- function(myfunc,mystate,dim1,dim2){
- tmp <- vector(dim1*dim2,mode="expression")
- for(i in 1:dim1){
- for(j in 1:dim2){
- tmp[(i-1)*dim2+j] <- parse(text=deparse(D(myfunc[i],mystate[j])))
- }
- }
- return(tmp)
- }
-
- # function to return symbolic derivatives of myfunc by mystate(single state)
- Derivation.scalar <- function(myfunc,mystate,dim){
- tmp <- vector(dim,mode="expression")
- for(i in 1:dim){
- tmp[i] <- parse(text=deparse(D(myfunc[i],mystate)))
- }
- return(tmp)
- }
-
- # function to solve Y_{t} (between (13.9) and (13.10)) using runge kutta method. Y_{t} is GL(d) valued (matrices)
- Get.Y <- function(){
- ## init
- dt <- Terminal/division
- assign(pars,0) ## epsilon=0
- Yinit <- as.vector(diag(d.size))
- Yt <- Yinit
- Y <- Yinit
- k <- numeric(d.size*d.size)
- k1 <- numeric(d.size*d.size)
- k2 <- numeric(d.size*d.size)
- k3 <- numeric(d.size*d.size)
- k4 <- numeric(d.size*d.size)
- Ystate <- paste("y",1:(d.size*d.size),sep="")
- F <- NULL
- F.n <- vector(d.size,mode="expression")
- for(n in 1:d.size){
- for(i in 1:d.size){
- F.tmp <- dx.drift[((i-1)*d.size+1):(i*d.size)]
- F.n[i] <- parse(text=paste(paste(F.tmp,"*",Ystate[((n-1)*d.size+1):(n*d.size)],sep=""),collapse="+"))
- }
- F <- append(F,F.n)
- }
- ## runge kutta
- for(t in 1:division){
- ## Xt
- for(i in 1:d.size){
- assign(state[i],X.t0[t,i]) ## state[i] is x_i, for example.
- }
- ## k1
- for(i in 1:(d.size*d.size)){
- assign(Ystate[i],Yt[i])
- }
-
- for(i in 1:(d.size*d.size)){
- k1[i] <- dt*eval(F[i])
- }
- ## k2
- for(i in 1:(d.size*d.size)){
- assign(Ystate[i],Yt[i]+k1[i]/2)
- }
- for(i in 1:(d.size*d.size)){
- k2[i] <- dt*eval(F[i])
- }
- ## k3
- for(i in 1:(d.size*d.size)){
- assign(Ystate[i],Yt[i]+k2[i]/2)
- }
- for(i in 1:(d.size*d.size)){
- k3[i] <- dt*eval(F[i])
- }
- ## k4
- for(i in 1:(d.size*d.size)){
- assign(Ystate[i],Yt[i]+k3[i])
- }
- for(i in 1:(d.size*d.size)){
- k4[i] <- dt*eval(F[i])
- }
- ## F(Y(t+dt))
- k <- (k1+k2+k2+k3+k3+k4)/6
- Yt <- Yt+k
- Y <- rbind(Y,Yt)
- }
- ## return matrix : (division+1)*(d.size*d.size)
- rownames(Y) <- NULL
- colnames(Y) <- Ystate
- return(ts(Y,deltat=dt,start=0))
- }
-
- # function to calculate Y_{t}Y_{s}^{-1}
- Get.YtYis <- function(t,s,range){
- yt <- matrix(Y[(range[t]-1)*delta/deltat(Y)+1,],d.size,d.size)
- yis <- solve(matrix(Y[(range[s]-1)*delta/deltat(Y)+1,],d.size,d.size))
- return(yt%*%yis)
- }
-
- # function to calculate lambda_{t,s}
- ## require: de.diffusion, ytyis
- Get.lambda.ts <- function(t,s,range){
- tmp <- matrix(0,d.size,r.size)
- assign(pars,0) ## epsilon=0
- #ytyis <- Get.YtYis((range[t]-1)*delta,(range[s]-1)*delta)
- for(i in 1:d.size){
- assign(state[i],X.t0[(range[s]-1)*delta/deltat(X.t0)+1,i])
- }
- for(i in 1:d.size){
- for(j in 1:r.size){
- tmp[i,j] <- eval(de.diffusion[[i]][j]) # dV/de
- }
- }
- return(ytyis[t,s,,]%*%tmp)
- }
-
- # function to calculate lambda_{t,s,0}
- ## require: de.drift, ytyis
- Get.lambda.ts0 <- function(t,s,range){
- tmp <- seq(0,0,length=d.size)
- assign(pars,0) ## epsilon=0
- #ytyis <- Get.YtYis((range[t]-1)*delta,(range[s]-1)*delta)
- for(i in 1:d.size){
- assign(state[i],X.t0[(range[s]-1)*delta/deltat(X.t0)+1,i])
- }
- for(i in 1:d.size){
- tmp[i] <- eval(de.drift[i]) # dV0/de
- }
- return(ytyis[t,s,,]%*%tmp)
- }
-
- # function to calculate mu_{i,t,s}
- ## require: ytyis
- Get.mu.its <- function(i.state,t,s,range){
- tmp <- matrix(0,d.size,r.size)
- assign(pars,0) ## epsilon=0
- #ytyis <- Get.YtYis((range[t]-1)*delta,(range[s]-1)*delta)
- for(i in 1:d.size){
- assign(state[i],X.t0[(range[s]-1)*delta/deltat(X.t0)+1,i])
- }
- # make expression of dV/di
- diV <- as.list(NULL)
- for(i in 1:d.size){
- diV[i] <- list(Derivation.scalar(V[[i]],state[i.state],r.size))
- }
- # make expression of d(dV/di)/de
- dideV <- as.list(NULL)
- for(i in 1:d.size){
- dideV[i] <- list(Derivation.scalar(diV[[i]],pars,r.size))
- }
- # evaluate expression
- for(i in 1:d.size){
- for(j in 1:r.size){
- tmp[i,j] <- eval(dideV[[i]][j])
- }
- }
- return(ytyis[t,s,,]%*%tmp)
- }
-
- # function to calculate mu_{i,t,s,0}
- ## require: ytyis
- Get.mu.its0 <- function(i.state,t,s,range){
- tmp <- seq(0,0,length=d.size)
- assign(pars,0) ## epsilon=0
- #ytyis <- Get.YtYis((range[t]-1)*delta,(range[s]-1)*delta)
- for(i in 1:d.size){
- assign(state[i],X.t0[(range[t]-1)*delta/deltat(X.t0)+1,i])
- }
- diV0 <- Derivation.scalar(V0,state[i.state],d.size) # dV0/di
- dideV0 <- Derivation.scalar(V0,pars,d.size) #d(dV0/di)/de
- for(i in 1:d.size){
- tmp[i] <- eval(dideV0[i])
- }
- return(ytyis[t,s,,]%*%tmp)
- }
-
- # function to calculate nu_{i,j,t,s}
- ## require: ytyis
- Get.nu.ijts <- function(i.state,j.state,t,s,range){
- tmp <- seq(0,0,length=d.size)
- assign(pars,0) ## epsilon=0
- #ytyis <- Get.YtYis((range[t]-1)*delta,(range[s]-1)*delta)
- for(i in 1:d.size){
- assign(state[i],X.t0[(range[s]-1)*delta/deltat(X.t0)+1,i])
- }
- diV0 <- Derivation.scalar(V0,state[i.state],d.size) #dV0/di
- didjV0 <- Derivation.scalar(diV0,state[j.state],d.size) #d(dV0/di)/dj
- for(i in 1:d.size){
- tmp[i] <- eval(didjV0[i])
- }
- return(ytyis[t,s,,]%*%tmp)
- }
-
- # function to calculate nu_{t,s}
- ## require: dede.diffusion, ytyis
- Get.nu.ts <- function(t,s,range){
- tmp <- matrix(0,d.size,r.size)
- assign(pars,0) ## epsilon=0
- #ytyis <- Get.YtYis((range[t]-1)*delta,(range[s]-1)*delta)
- for(i in 1:d.size){
- assign(state[i],X.t0[(range[s]-1)*delta/deltat(X.t0)+1,i])
- }
- for(i in 1:d.size){
- for(j in 1:r.size){
- tmp[i,j] <- eval(dede.diffusion[[i]][j])
- }
- }
- return(ytyis[t,s,,]%*%tmp)
- }
-
- # function to calculate nu_{t,s,0}
- ## require: dede.drift, ytyis
- Get.nu.ts0 <- function(t,s,range){
- tmp <- seq(0,0,length=d.size)
- assign(pars,0) ## epsilon=0
- #ytyis <- Get.YtYis((range[t]-1)*delta,(range[s]-1)*delta)
- for(i in 1:d.size){
- assign(state[i],X.t0[(range[s]-1)*delta/deltat(X.t0)+1,i])
- }
- for(i in 1:d.size){
- tmp[i] <- eval(dede.drift[i])
- }
- return(ytyis[t,s,,]%*%tmp)
- }
-
- # function to calculate mu in thesis p5
- funcmu <- function(e=0){
- division <- nrow(X.t0)
- XT <- X.t0[division,] #data X0 observated last. size:vector[d.size]
-
- ## calculate derived F by e with XT and e=0. deF(XT,0)
- deF0 <- c() #size:vector[k.size]
- for(d in 1:d.size){
- assign(state[d],XT[d]) #input XT in state to use eval function
- }
- for(k in 1:k.size){
- deF <- deriv(F[k],"e") #expression of derived F by e
- deF0[k] <- attr(eval(deF),"gradient") #calculate deF(derived F by e) with XT
- }
-
- ## calculate derived f0 by e with Xt and e=0. def0(Xt,0)
- def0 <- matrix(0,k.size,division) #size:matrix[k.size,division]
- for(k in 1:k.size){
- def <- deriv(f[[1]][k],"e") #expression of derived f0 by e
- for(t in 1:division){
- X0t <- X.t0[t,] #data X0 observated on time t
- for(d in 1:d.size){
- assign(state[d],X0t[d]) #input X0t in state
- }
- def0[k,t] <- attr(eval(def),"gradient") #calculate def(derived f0 by e) with X0t
- }
- }
-
- # integrate def0 (just sum it)
- def0 <- apply(def0,1,sum) #sum of def0. size:vector[k.size]
- def0 <- def0*(1/(division-1))
- mu <- def0+deF0 #size:vector[k.size]
- return(mu)
- }
-
- # function to calculate a_{s}^{alpha} in bookchapter p5
- funca <- function(e=0){
- #init
- division <- nrow(X.t0)
- XT <- X.t0[division,] #data X0 observated last. size:vector[d.size]
- defa <- array(0,dim=c(k.size,r.size,division)) #size:array[k.size,r.size,division]
- deva <- array(0,dim=c(d.size,r.size,division)) #size:array[d.size,r.size,division]
- dxF0 <- matrix(0,k.size,d.size) #size:matrix[k.size,d.size]
- dxf0 <- array(0,dim=c(k.size,d.size,division)) #size:array[k.size,d.size,division]
- dxf <- c()
- dxF <- c()
- def <- list()
- dev <- list()
-
- # prepare expression of derivatives
- for(k in 1:k.size){
- dxf[k] <- deriv(f[[1]][k],state) #expression of d f0/dx
- dxF[k] <- deriv(F[k],state) #expression of d F/dx
- def[[k]] <- list()
- for(r in 2:(r.size+1)){
- def[[k]][[r-1]] <- deriv(f[[r]][k],"e") #expression of derived fa by e
- }
- }
- for(r in 1:r.size){
- dev[[r]] <- list()
- for(d in 1:d.size){
- dev[[r]][[d]] <- deriv(V[[d]][r],"e") #expression of derived Vr by e
- }
- }
-
- # evaluate derivative expressions
- for(t in 1:division){
- X0t <- X.t0[t,]
- for(d in 1:d.size){
- assign(state[d],X0t[d]) #input X0t in state to use eval function to V
- }
- for(k in 1:k.size){
- ##calculate derived f0 by x with Xt and e=0. dxf0(Xt,0)
- dxf0[k,,t] <- attr(eval(dxf[k]),"gradient") #calculate dxf(derived f0 by x) with X0t
- ##calculate derived F by x with XT and e=0. dxF(XT,0)
- dxF0[k,] <- attr(eval(dxF[k]),"gradient") #calculate dxF(derived F by e) with XT
- for(r in 2:(r.size+1)){
- ##calculate derived fa by e with Xt and e=0. defa(Xt,0)
- defa[k,r-1,t] <- attr(eval(def[[k]][[r-1]]),"gradient") #calculate def(derived fa by e) with X0t
- }
- }
- for(r in 1:r.size){
- for(d in 1:d.size){
- ##calculate derived Va by e with Xt and e=0. deVa(Xt,0)
- deva[d,r,t] <- attr(eval(dev[[r]][[d]]),"gradient") #calculate dev(derived Vr by e) with X0t
- }
- }
- }
-
- # prepare Y and Y^{-1}
- arrayY <- array(0,dim=c(d.size,d.size,division))
- invY <- array(0,dim=c(d.size,d.size,division))
- for(t in 1:division){
- arrayY[,,t] <- matrix(Y[t,],d.size,d.size)
- invY[,,t] <- solve(arrayY[,,t])
- }
-
- # calculate dxF*Y^{T}*Y^{-1}*deV_{a}
- second <- array(0,dim=c(k.size,r.size,division))
- temp <- dxF0 %*% arrayY[,,division]
- for(t in 1:division) {
- second[,,t] <- temp %*% invY[,,t] %*% deva[,,t]
- }
-
- #calculate integral
- fIntegral <- array(0,dim=c(k.size,r.size,division))
- third <- array(0,dim=c(k.size,r.size,division))
- dt <- Terminal / division
- third[,,division] <- dxf0[,,division] %*% arrayY[,,division] %*% invY[,,division] %*% deva[,,division] * dt
- for(s in (division-1):1) {
- third[,,s] <- third[,,s+1] + dxf0[,,s] %*% arrayY[,,s] %*% invY[,,s] %*% deva[,,s] * dt
- }
-
- # defa <- aperm(defa,c(1,3,2))*1.0
- # deva <- aperm(deva,c(1,3,2))*1.0
- # dxF0 <- dxF0*1.0
- # dxf0 <- dxf0*1.0
-
- ##use C source
- # dyn.load("yuima.so")
- # a <- .Call("get_a",defa,dxF0,arrayY,invY,deva,dxf0,
- # 1.0*dim(defa),1.0*dim(arrayY),1.0*dim(invY),1.0*dim(deva),1.0*dim(dxf0),
- # 1.0*a,1.0*dim(a))
-
- return(defa + second + third) #size:array[k.size,r.size,division]
- }
-
- # function to calculate sigma in thesis p5
- # require: aMat
- funcsigma <- function(e=0){
-
- division <- nrow(X.t0)
- sigma <- matrix(0,k.size,k.size) #size:matrix[k.size,k.size]
- for(t in 1:division){
-## sigma <- sigma+(aMat[,,t]%*%t(aMat[,,t])) /(division-1) #calculate sigma: old code, modified 20110519
- sigma <- sigma+(t(aMat[,,t])%*%(aMat[,,t])) /(division-1) #calculate sigma
- }
- if(any(eigen(sigma)$value<=0.0001)){
- # Singularity check
- yuima.warn("Eigen value of covariance matrix in very small.")
- }
- return(sigma)
- }
-
- ## integrate start:1 end:t number to integrate:block
- # because integration at all 0:T takes too much time,
- # we sample only 'block' points and use trapezium rule to integrate
- make.range.for.trapezium.fomula <- function(t,block){
- if(t/block <= 1){ # block >= t : just use all points
- range <- c(1:t)
- }else{ # make array which includes points to use
- range <- as.integer( (c(0:block) * (t/block))+1)
- if( range[block+1] < t){
- range[block+2] <- t
- }else if( range[block+1] > t){
- range[block+1] <- t
- }
- }
- return(range)
- }
-
- # function to return expressions of df0/dxi
- deriv.f0.for.state<- function(f0){
- tmp_deriv_f0 <- function(i,f){
- d_xi_f <- c()
- for(j in 1:k.size){
- d_xi_f[j] <- deriv(f[j],state[i])
- }
- return(d_xi_f)
- }
- tmp <- apply(as.matrix(1:length(state)),1,tmp_deriv_f0,f0)
- ##return list of (df0/dxi)
- return(tmp)
- }
-
- ## This function return value of expr(X0[t])
- ## expr:list of derived for state x1,x2,...
- ## ex) expr[[1]] : f0 derived for x1
- ## t: time index (t=1,2,...,division+1)
- input.deriv <- function(t,expr,l=1,X0){
- df <- c()
- ##input x1,x2,...,
- for(i in 1:d.size){
- df[i] <- expr[[i]][l]
- assign(state[i],X0[t,i])
- }
- ##epsilon = 0
- assign(pars[1],0)
- tmp <- c()
- for(i in 1:d.size){
- tmp[i] <- attr(eval(df[i]),"gradient")
- }
- return(tmp)
- }
-
- ## get hessian for(state,"e")
- hessian.f0.di.de<- function(f0){
- tmp.hessian.f0 <- function(i,f){
- d.xi.de.f <- c()
- for(j in 1:k.size){
- d.xi.de.f[j] <- deriv(f[j],c(state[i],"e"),hessian=T)
- }
- return(d.xi.de.f)
- }
- tmp <- apply(as.matrix(1:length(state)),1,tmp.hessian.f0,f0)
- ##return list of (d^2 f0/dxi de)
- return(tmp)
- }
-
- hessian.f.dxi.de<- function(f){
- list_l_dxi_de <- list(NULL)
- list_dxi_de <- list(NULL)
- de <- c()
- for(k in 1:k.size){
- for(i in 1:d.size){
- for(r in 1:r.size){
- de[r] <- deriv(f[[r+1]][k],c(state[i],"e"),hessian=T)
- }
- list_dxi_de[[i]] <- de
- }
- list_l_dxi_de[[k]] <- list_dxi_de
- }
- ##return list of (d^2 f0/dxi de)
- return(list_l_dxi_de)
- }
-
- ## This function return list deriv f0 for (xi,xj)
- ## list_k_dxi_dxj[[ k ]][[ i ]][ j ] is expression f0[k] derived for (state[i],state[j])
- ## f0:1,...,k.size expression
- hessian.f0.di.dj<- function(f0){
- list_k_dxi_dxj <- as.list(NULL)
- list_dxi_dxj<- as.list(NULL)
- dxj <- c()
- for(k in 1:k.size){
- for(i in 1:d.size){
- for(j in 1:d.size){
- dxj[j] <- deriv(f0[k],c(state[i],state[j]),hessian=T)
- }
- list_dxi_dxj[[i]] <- dxj
- }
- list_k_dxi_dxj[[k]] <- list_dxi_dxj
- }
- return(list_k_dxi_dxj)
- }
-
- # following funcs (named 'input.hessian~') solve expression of hessian
- ## This function return (d x 1 vector)
- input.hessian <- function(t,expr,n=1,m=2,l=1,X0){
- h_f <- c()
- ##input x1,x2,...,
- for(i in 1:d.size){
- h_f[i] <- expr[[i]][l]
- assign(state[i],X0[t,i])
- }
- ##epsilon = 0
- assign(pars[1],0)
- tmp <- c()
- for(i in 1:d.size){
- tmp[i] <- attr(eval(h_f[i]),"hessian")[1,n,m]
- }
- return(tmp)
- }
-
- ## This function returns (d x d matrix)
- ## t: time index expr:derived expression
- input.hessian.dxi.dxj <- function(t,expr,n=1,m=2,l=1,X0){
- h_f <- list(NULL)
- ##input x1,x2,...,
- for(i in 1:d.size){
- h_f[[i]] <- expr[[l]][[i]]
- assign(state[i],X0[t,i])
- }
- ##epsilon = 0
- assign(pars[1],0)
-
- tmp <- matrix(0,d.size,d.size)
- for(i in 1:d.size){
- for(j in 1:d.size){
- if(i == j){
- tmp[i,j] <- attr(eval(h_f[[i]][j]),"hessian")[1,1,1]
- }else{
- tmp[i,j] <- attr(eval(h_f[[i]][j]),"hessian")[1,i,j]
- }
- }
- }
- dim(tmp) <- c()
- return(tmp)
- }
-
- ## This function returns (d x r matrix)
- ## t:index , expr:hessian for f, l=1,...,k.size
- input.hessian.dxi.de.f<- function(t,expr,n=1,m=2,l=1,X0){
- h.f <- list(NULL)
- h.f.r <- c()
- for(i in 1:d.size){
- for(r in 1:r.size){
- h.f.r[r] <- expr[[l]][[i]][r]
- }
- assign(state[i],X0[t,i])
- h.f[[i]] <- h.f.r
- }
- assign(pars[1],0)
- tmp <- matrix(0,d.size,r.size)
- for(d in 1:d.size){
- for(r in 1:r.size){
- tmp[d,r] <- attr(eval(h.f[[d]][r]),"hessian")[1,n,m]
- }
- }
- return(tmp)
- }
-
- ## multi dimension gausian distribusion
- normal <- function(x,mu,Sigma){
- if(length(x)!=length(mu)){
- print("Error:length of x != one of mu")
- }
- dimension <- length(x)
- tmp <- 1/((sqrt(2*pi))^dimension * sqrt(det(Sigma))) * exp(-1/2 * t((x-mu)) %*% solve(Sigma) %*% (x-mu) )
- return(tmp)
- }
-
- ## get d0
- ## required library(adapt)
- get.d0.term <- function(){
- lambda.max<- max(eigen(Sigma)$values)
- ## get g(z)*pi0(z)
- gz_pi0 <- function(z){
- return( g(z) * H0 *normal(z,mu=mu,Sigma=Sigma))
- }
- gz_pi02 <- function(z){
- return( g(z) * H0 *dnorm(z,mean=mu,sd=sqrt(Sigma)))
- }
-
- ## integrate
- if( k.size ==1){ # use 'integrate' if k=1
- tmp <- integrate(gz_pi02,-Inf,Inf)$value
- }else if( 2 <= k.size || k.size <= 20 ){ # use library 'adapt' to solve multi-dimentional integration
- max <- 10 * lambda.max
- min <- -10 * lambda.max
- L <- (max - min)
-# if(require(adapt)){
- if(require(cubature)){
-#tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi0)$value
- tmp <- adaptIntegrate(gz_pi0, lower=rep(min,k.size),upper=rep(max,k.size))$integral
- } else {
- tmp <- NA
- }
- }else{
- stop("length k is too big.")
- }
- return(tmp)
- }
-
- ###############################################################################
- # following funcs are part of d1 term
- # because they are finally integrated at 'get.d1.term()',
- # these funcs are called over and over again.
- # so, we use trapezium rule for integration to save time and memory.
-
- # these funcs almost calculate each formulas including trapezium integration.
- # see each formulas in thesis to know what these funcs do.
-
- # some funcs do alternative calculation at k=1.
- # it depends on 'integrate()' function
- ###############################################################################
-
- # p.9 Lemma2 (a)
-
- ## This function returns First term of di.bar (d.size x block martix) and
- ## part of second term(Second.tmp)
- ## Second.tmp is ((d x block) x k) matrix
- ## aMat.tmp is k x (d x block) matrix
- get.di.bar.init <- function(){
- # trapezium rule
- tmp.mat <- rep(1,block)
- tmp.mat2 <- rep(Diff,d.size)
- dim(tmp.mat2) <- c((block*block),d.size)
- tmp.mat2 <- t(tmp.mat2)
- dim(tmp.mat2) <- c((d.size*block),block)
- First <- (lambda.ts0 * tmp.mat2) %*% tmp.mat * delta
- dim(First) <- c(d.size,block)
-
- tmp.mat3 <- tmp.mat2
- for(i in 1:r.size){
- if(i != 1){
- tmp.mat3 <- rbind(tmp.mat3,tmp.mat2)
- }
- }
- dim(tmp.mat3) <- c(d.size*block,r.size*block)
- Second.tmp <- (lambda.ts * tmp.mat3) %*% t(aMat.tmp) * delta
- tmp <- list(First=First,
- Second.tmp=Second.tmp)
- return(tmp)
- }
-
- ## dependency:dat.di.bar
- get.di.bar <- function(x){
- if(k.size ==1){
- First <- dat.di.bar$First
-
- tmp.di.bar <- dat.di.bar$Second.tmp
- dim(tmp.di.bar) <- c(d.size,block)
- Second.tmp <- tmp.di.bar
- for(i in 1:(length(x))){
- if(i!=1){
- First <- rbind(First,dat.di.bar$First)
- }
- }
- for(i in 1:length(x)){
- if(i!=1){
- Second.tmp <- rbind(Second.tmp,tmp.di.bar)
- }
- }
- tmp.x <- x
- for(i in 1:d.size){
- if(i !=1){
- tmp.x <- rbind(tmp.x,x)
- }
- }
- dim(tmp.x) <- c()
- tmp.x <- rep(tmp.x,block)
- dim(tmp.x) <- c((d.size*length(x)),block)
- Second <- tmp.x * Second.tmp / as.double(Sigma)
- tmp <- First + Second
- }else{
- Second <- dat.di.bar$Second.tmp %*% solve(Sigma) %*% x
- dim(Second) <- c(d.size,block)
- tmp <- dat.di.bar$First + Second
- }
- return(tmp)
- }
-
- ## h is (d x block matrix)
- ## x is k dimension vector
- get.Di.bar <- function(h,x){
- if(k.size==1){
- tmp.Diff <- Diff[block,]
- for(i in 1:d.size){
- if( i != 1 ){
- tmp.Diff <- rbind( tmp.Diff,Diff[block,] )
- }
- }
- tmp.h <- h * tmp.Diff
- for(i in 1:length(x)){
- if(i !=1){
- tmp.h <- rbind( tmp.h , (h*tmp.Diff) )
- }
- }
- tmp <- tmp.h * get.di.bar(x) *delta
- tmp <- tmp %*% rep(1,block)
- dim(tmp) <- c(d.size,length(x))
- }else{
- tmp.Diff <- Diff[block,]
- for(i in 1:d.size){
- if( i != 1 ){
- tmp.Diff <- rbind( tmp.Diff,Diff[block,] )
- }
- }
- tmp <- h * get.di.bar(x) * tmp.Diff * delta
- tmp <- as.vector(tmp %*% rep(1,block))
- }
- return(tmp)
- }
-
- # p.10 (b)
-
- ## this function returns first term and part of second term of Di
- ## h: d x (r x block) matrix
- ## x: k dimension vector
- ## dependency:dat.di.bar
- get.Di.init <- function(h){
- ## First term
- tmp1 <- dat.di.bar$First
- for(i in 1:r.size){
- if(i !=1){
- tmp1 <- rbind( tmp1 , dat.di.bar$First )
- }
- }
- dim(tmp1) <- c( d.size , r.size * block)
- tmp.Diff <- Diff[block,]
- for(i in 1:(k.size * r.size)){
- if(i != 1) tmp.Diff <- rbind(tmp.Diff,Diff[block,])
- }
- dim(tmp.Diff) <- c(k.size,(r.size*block))
- First <- (tmp1 * h) %*% t(aMat.tmp * tmp.Diff) %*% solve(Sigma) * delta
- ## End of First term
-
- tmp.mat1 <- t(dat.di.bar$Second.tmp)
- dim(tmp.mat1) <- c((k.size * d.size),block)
- Second.tmp <- array(0,dim=c(k.size,k.size,d.size))
- for( i in 1:d.size){
- tmp.mat2 <- tmp.mat1[((i-1)*k.size + 1):(i*k.size),]
- dim(tmp.mat2) <- c()
- tmp.mat3 <- tmp.mat2
- tmp2 <- h[i,]
- dim(tmp2) <- c(r.size,block)
- tmp.mat4 <- tmp2
- for(j in 1:r.size){
- if(j != 1) tmp.mat3 <- rbind( tmp.mat3 , tmp.mat2 )
- }
- for(j in 1:k.size){
- if(j !=1) tmp.mat4 <- rbind( tmp.mat4 , tmp2 )
- }
- dim(tmp.mat4) <- c(r.size,(k.size * block))
- tmp3 <- (tmp.mat4 * tmp.mat3)
- tmp4 <- matrix(0,(r.size*block),k.size)
- for(k in 1:k.size){
- tmp4[,k] <- tmp3[,(c(1:block)*k.size + ( k - k.size))]
- }
- tmp.Diff <- Diff[block,]
- for( j in 1:(k.size*r.size)){
- if(j !=1){
- tmp.Diff <- rbind(tmp.Diff,Diff[block,])
- }
- }
- dim(tmp.Diff)<-c(k.size,(r.size*block))
- Second.tmp[,,i] <- (aMat.tmp * tmp.Diff) %*% tmp4 * delta
- }
- tmp <- list(First=First,
- Second.tmp=Second.tmp
- )
- }
-
- ## dependency: dat.Di.init include First and Second.tmp
- get.Di <- function(dat.Di.init,x){
- if(k.size==1){
- tmp <- numeric(d.size)
- for(i in 1:d.size){
- tmp[i] <- dat.Di.init$Second.tmp[,,i]
- }
- ##dat.Di.init$First is d x k matrix
- First <- dat.Di.init$First %*% x / as.double(Sigma)
- Second <- (tmp/(as.double(Sigma)^2) ) %*% t(x^2 - as.double(Sigma))
- }else{
- First <- dat.Di.init$First %*% solve(Sigma) %*% x
- Second <- numeric(d.size)
- for(i in 1:d.size){
- tmp <- dat.Di.init$Second.tmp[,,i]
- Second[i] <- sum( diag( solve(Sigma) %*% tmp %*% solve(Sigma) %*% (x%*%t(x) - Sigma) ) )
- }
- }
- Di <- First + Second
- return(Di)
- }
-
- # p.10 (c)
-
- ## preparation to calculate d^{ij}(x)_{t}
- ## dependency: dat.di.bar
- get.dij.init <- function(){
- ## First term
- tmp1.1<- dat.di.bar$First
- tmp1.2 <- rep(dat.di.bar$First,d.size)
- dim(tmp1.2) <- c((d.size*block),d.size)
- tmp1.2 <- t(tmp1.2)
- dim(tmp1.2) <- c((d.size^2),block)
- for(i in 1:d.size){
- if(i !=1){
- tmp1.1 <- rbind( tmp1.1 , dat.di.bar$First )
- }
- }
- First <- tmp1.1 * tmp1.2 ## dim(First) = c( d.size^2 , block)
-
- ## Second term
- tmp2.1 <- tmp1.1
- dim(tmp2.1) <- c()
- tmp2.1 <- rep(tmp2.1,k.size)
- dim(tmp2.1) <- c((d.size^2)*block,k.size)
- tmp2.2 <- rep(dat.di.bar$Second.tmp,d.size)
- dim(tmp2.2) <- c(d.size*k.size*block,d.size)
- tmp2.2 <- t(tmp2.2)
- dim(tmp2.2) <- c(d.size^2*block,k.size)
- Second.tmp <- (tmp2.1 * tmp2.2) %*% solve(Sigma) ## dim(Second.tmp) = c(d.size^2 * block ,k.size)
-
- ## Third term
- tmp3.1 <- rep(dat.di.bar$First,d.size)
- dim(tmp3.1) <- c(d.size*block,d.size)
- tmp3.1 <- as.vector(t(tmp3.1))
- tmp3.1 <- rep(tmp3.1,k.size)
- dim(tmp3.1) <- c((d.size^2)*block,k.size)
-
- tmp3 <- dat.di.bar$Second.tmp
- dim(tmp3) <- c(d.size,k.size*block)
- tmp3.2 <- tmp3
- for( i in 1:d.size){
- if(i != 1){
- tmp3.2 <- rbind(tmp3.2,tmp3)
- }
- }
- dim(tmp3.2) <- c((d.size^2)*block,k.size)
- Third.tmp <- (tmp3.1 * tmp3.2) %*% solve(Sigma) ## dim(Third.t,p) = c(d.size^2 * block ,k.size)
-
- ##Fourth term
- tmp4 <- t(dat.di.bar$Second.tmp)
- dim(tmp4) <- c(k.size*d.size,block)
- tmp4.1 <- tmp4
- for(i in 1:d.size){
- if(i != 1) tmp4.1 <- rbind(tmp4.1,tmp4)
- }
- # rm(tmp4)
- dim(tmp4.1) <- c(k.size,(d.size^2)*block)
- tmp4.2 <- tmp4.1
- for(k in 1:k.size){
- if(k != 1) tmp4.2 <- rbind(tmp4.2,tmp4.1)
- }
- # rm(tmp4.1)
- dim(tmp4.2) <- c(k.size,k.size*(d.size^2)*block)
-
- tmp4.3 <- t(dat.di.bar$Second.tmp)
- for(i in 1:d.size){
- if(i != 1) tmp4.3 <- rbind(tmp4.3,t(dat.di.bar$Second.tmp))
- }
- dim(tmp4.3) <- c()
- tmp4.4 <- tmp4.3
- for( k in 1:k.size){
- if( k != 1) tmp4.4 <- rbind(tmp4.4,tmp4.3)
- }
- # rm(tmp4.3)
- dim(tmp4.4) <- c(k.size,k.size*(d.size^2)*block)
- tmp4.5 <- tmp4.2 * tmp4.4
- tmp4.6 <- matrix(0,k.size*(d.size^2)*block,k.size)
- for(k in 1:k.size){
- tmp4.6[,k] <- tmp4.5[,(1:(d.size^2 * block)-1)*k.size + k ]
- }
- tmp4.6 <- t(tmp4.6)
- tmp4.7 <- solve(Sigma) %*% (tmp4.5 + tmp4.6)
- # rm(tmp4.5)
- # rm(tmp4.6)
- tmp4.8 <- matrix(0,k.size*(d.size^2)*block,k.size)
- for(k in 1:k.size){
- tmp4.8[,k] <- tmp4.7[,(1:(d.size^2 * block)-1)*k.size + k ]
- }
- Fourth.tmp <- tmp4.8 %*% solve(Sigma) ## dim(Fourth.tmp) =c( k.size * d.size^2 * block , k.size)
-
- ## Fifth term
- tmp.Diff <- rep(Diff,d.size)
- dim(tmp.Diff) <- c(block^2,d.size)
- tmp.Diff <- t(tmp.Diff)
- dim(tmp.Diff) <- c(d.size*block,block)
- tmp.Diff2 <- tmp.Diff
- for(r in 1:r.size){
- if(r != 1) tmp.Diff <- rbind(tmp.Diff,tmp.Diff2)
- }
- dim(tmp.Diff) <- c(d.size*block,r.size*block)
- # rm(tmp.Diff2)
- Fifth <- matrix(0,d.size^2,block)
-
- for(t in 1:block){
- start <- (t-1) * d.size + 1
- end <- (t-1) * d.size +d.size
- if( d.size == 1 ){
- Fifth[,t] <- (lambda.ts[start:end,] * tmp.Diff[start:end,]) %*% lambda.ts[start:end,] *delta
- }else{
- Fifth[,t] <- (lambda.ts[start:end,] * tmp.Diff[start:end,]) %*% t(lambda.ts[start:end,]) *delta
- }
- }
- ## dim(Fifth) = c( d.size^2 ,block )
- tmp <- list(First=First,
- Second.tmp=Second.tmp,
- Third.tmp=Third.tmp,
- Fourth.tmp=Fourth.tmp,
- Fifth=Fifth
- )
- return(tmp)
- }
-
- ## dependency: dat.dij
- get.dij <- function(x){
- if( k.size == 1 ){
- First <- dat.dij$First
- First.tmp <- First
-
- Fifth <- dat.dij$Fifth
- Fifth.tmp <- Fifth
-
- Second <- dat.dij$Second.tmp
- dim(Second) <- c(d.size^2,block)
- Second.tmp <- Second
- for( i in 1:length(x)){
- if(i != 1){
- First <- rbind(First,First.tmp)
- Second <- rbind(Second,Second.tmp)
- Fifth <- rbind(Fifth,Fifth.tmp)
- }
- }
- # rm(First.tmp)
- # rm(Second.tmp)
- # rm(Fifth.tmp)
-
- tmp.x <- rep(x,d.size^2 * block)
- dim(tmp.x) <- c(length(x),d.size^2 * block)
- tmp.x <- t(tmp.x)
- dim(tmp.x) <- c(block,d.size^2*length(x))
- tmp.x <- t(tmp.x)
- Second <- Second * tmp.x
-
- Third <- dat.dij$Third.tmp
- dim(Third) <- c(d.size^2,block)
- Third.tmp <- Third
- for( i in 1:length(x)){
- if(i != 1) Third <- rbind(Third,Third.tmp)
- }
- # rm(Third.tmp)
- Third <- Third * tmp.x
-
- Fourth <- dat.dij$Fourth.tmp
- dim(Fourth) <- c(d.size^2,block)
- Fourth.tmp <- Fourth
- for( i in 1:length(x)){
- if(i != 1) Fourth <- rbind(Fourth,Fourth.tmp)
- }
- # rm(Fourth.tmp)
- tmp.x <- rep((x^2 - Sigma),d.size^2 * block)
- dim(tmp.x) <- c(length(x),d.size^2 * block)
- tmp.x <- t(tmp.x)
- dim(tmp.x) <- c(block,d.size^2*length(x))
- tmp.x <- t(tmp.x)
- Fourth <- Fourth * tmp.x /2
-
- tmp <- First + Second + Third + Fourth + Fifth
- }else{
- Second <- dat.dij$Second.tmp %*% x
- dim(Second) <- c(d.size^2,block)
-
- Third <- dat.dij$Third.tmp %*% x
- dim(Third) <- c(d.size^2,block)
-
- tmp.x <- rep((x %*% t(x)) - Sigma, d.size^2 * block )
- dim(tmp.x) <- c(k.size , k.size * d.size^2 * block )
-
- Fourth <- rep(1,k.size) %*% t( dat.dij$Fourth.tmp * t(tmp.x) )
- dim(Fourth) <- c(k.size ,d.size^2 * block)
- Fourth <- ( rep(1,k.size) %*% Fourth ) /2
- dim(Fourth) <- c(d.size^2,block)
- tmp <- dat.dij$First + Second + Third + Fourth + dat.dij$Fifth
- }
- return( tmp )
- }
-
- get.Dij <- function(h,x){
- if(k.size == 1){
- tmp.Diff <- Diff[block,]
- for(i in 1:(d.size^2)){
- if( i != 1 ){
- tmp.Diff <- rbind( tmp.Diff,Diff[block,] )
- }
- }
- tmp.h <- h * tmp.Diff
- for(i in 1:length(x)){
- if(i !=1){
- tmp.h <- rbind( tmp.h , (h*tmp.Diff) )
- }
- }
- tmp <- tmp.h * get.dij(x) *delta
- tmp <- tmp %*% rep(1,block)
- dim(tmp) <- c(d.size^2,length(x))
- }else{
- tmp.Diff <- Diff[block,]
- for(i in 1:(d.size^2)){
- if( i != 1 ){
- tmp.Diff <- rbind( tmp.Diff,Diff[block,] )
- }
- }
- tmp <- h * get.dij(x) * tmp.Diff * delta
- tmp <- as.vector(tmp %*% rep(1,block))
- }
- return(tmp)
- }
-
- # p.11 (d)
-
- get.el.init <- function(){
- tmp.mat <- rep(1,block)
- tmp.mat2 <- rep(Diff,d.size)
- dim(tmp.mat2) <- c((block*block),d.size)
- tmp.mat2 <- t(tmp.mat2)
- dim(tmp.mat2) <- c((d.size*block),block)
- First <- (nu.ts0 * tmp.mat2) %*% tmp.mat * delta
-
- dim(First) <- c(d.size,block)
- tmp.mat3 <- tmp.mat2
- for(i in 1:r.size){
- if(i != 1){
- tmp.mat3 <- rbind(tmp.mat3,tmp.mat2)
- }
- }
- dim(tmp.mat3) <- c(d.size*block,r.size*block)
- Second.tmp <- (nu.ts * tmp.mat3) %*% t(aMat.tmp) * delta
- tmp <- list(Third=First,
- Fifth.tmp=Second.tmp)
- return(tmp)
- }
-
-
- ## x : k.size dimension vector
- get.el <- function(x){
- if(k.size==1){
- First <- matrix(0,d.size*length(x),block)
- Second <- matrix(0,d.size*length(x),block)
- Third <- matrix(0,d.size*length(x),block)
- Fourth <- matrix(0,d.size*length(x),block)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 155
More information about the Yuima-commits
mailing list