[Yuima-commits] r52 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Dec 26 09:34:16 CET 2009
Author: hinohide
Date: 2009-12-26 09:34:15 +0100 (Sat, 26 Dec 2009)
New Revision: 52
Added:
pkg/yuima/R/asymptotic_term.R
pkg/yuima/R/simFunctional.R
pkg/yuima/R/yuima.functional.R
pkg/yuima/man/asymptotic_term.Rd
pkg/yuima/man/setFunctional.Rd
pkg/yuima/man/simFunctional.Rd
pkg/yuima/man/yuima.functional-class.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/AllClasses.R
pkg/yuima/R/yuima.R
pkg/yuima/man/setYuima.Rd
pkg/yuima/man/yuima-class.Rd
Log:
add asymptotic expansion
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2009-12-24 07:38:56 UTC (rev 51)
+++ pkg/yuima/DESCRIPTION 2009-12-26 08:34:15 UTC (rev 52)
@@ -1,9 +1,9 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 0.0.81
-Date: 2009-12-24
-Depends: methods, zoo
+Version: 0.0.82
+Date: 2010-01-01
+Depends: methods, zoo, adapt
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Description: The YUIMA Project for Simulation and Inference
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2009-12-24 07:38:56 UTC (rev 51)
+++ pkg/yuima/NAMESPACE 2009-12-26 08:34:15 UTC (rev 52)
@@ -28,7 +28,17 @@
"rql",
"ml.ql",
"adaBayes",
- "limiting.gamma"
+ "limiting.gamma",
+
+
+ "getF",
+ "getf",
+ "getxinit",
+ "gete",
+ "simFunctional",
+ "F0",
+ "Fnorm",
+ "asymptotic_term"
)
## function which we want to expose to the user
@@ -58,3 +68,15 @@
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)
+
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2009-12-24 07:38:56 UTC (rev 51)
+++ pkg/yuima/R/AllClasses.R 2009-12-26 08:34:15 UTC (rev 52)
@@ -66,6 +66,18 @@
)
)
+# Class 'yuima.functional'
+
+# functional model used in 'asymptotic term' procedure
+
+setClass("yuima.functional", representation(F = "ANY",
+ f = "list",
+ xinit = "numeric",
+ e = "numeric"
+ )
+ )
+
+
# Class 'yuima'
# this is the principal class of yuima project. It may contain up to
@@ -80,6 +92,7 @@
setClass("yuima", representation(data = "yuima.data",
model = "yuima.model",
sampling = "yuima.sampling",
- characteristic = "yuima.characteristic"
+ characteristic = "yuima.characteristic",
+ functional = "yuima.functional"
)
)
Added: pkg/yuima/R/asymptotic_term.R
===================================================================
--- pkg/yuima/R/asymptotic_term.R (rev 0)
+++ pkg/yuima/R/asymptotic_term.R 2009-12-26 08:34:15 UTC (rev 52)
@@ -0,0 +1,1476 @@
+# in this source we note formulae like latex
+
+
+setGeneric("asymptotic_term",
+ function(yuima,block=100, rho, g)
+ standardGeneric("asymptotic_term")
+ )
+
+setMethod("asymptotic_term",signature(yuima="yuima"), function(yuima,block=100, rho, g){
+
+ 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)
+ e <- gete(yuima at functional)
+ Terminal <- yuima at sampling@Terminal
+ division <- yuima at sampling@n
+ 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)
+
+ X.t0 <- Get.X.t0(yuima)
+ delta <- deltat(X.t0)
+ pars <- 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
+ 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])
+ }
+ ## 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 thesis 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 derived f0 by x
+ dxF[k] <- deriv(F[k],state) #expression of derived F by x
+ 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
+ }
+ if(any(eigen(sigma)$value<=0.0001)){
+ # Singularity check
+ cat("Eigen value of covariance matrix in very small.\n")
+ }
+ 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)
+ }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)
+ tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi0)
+ }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 ){
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 52
More information about the Yuima-commits
mailing list