[Yuima-commits] r10 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 28 12:08:37 CET 2009
Author: hinohide
Date: 2009-10-28 12:08:28 +0100 (Wed, 28 Oct 2009)
New Revision: 10
Added:
pkg/yuima/R/limiting.gamma.R
pkg/yuima/man/limiting.gamma.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/man/quasi-likelihood.Rd
pkg/yuima/man/yuima-class.Rd
pkg/yuima/man/yuima.model-class.Rd
Log:
add limiting gamma function
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2009-10-28 10:20:11 UTC (rev 9)
+++ pkg/yuima/DESCRIPTION 2009-10-28 11:08:28 UTC (rev 10)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 0.0.67
+Version: 0.0.68
Date: 2009-10-28
Depends: methods, zoo
Author: Yuima project team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2009-10-28 10:20:11 UTC (rev 9)
+++ pkg/yuima/NAMESPACE 2009-10-28 11:08:28 UTC (rev 10)
@@ -15,10 +15,10 @@
exportMethods(
"dim",
"length",
- # "start",
+ ## "start",
"plot",
- # "time",
- # "end",
+ ## "time",
+ ## "end",
"simulate",
"cce",
@@ -28,10 +28,11 @@
"initialize",
- "ql",
- "rql",
- "ml.ql",
- "adaBayes"
+ "ql",
+ "rql",
+ "ml.ql",
+ "adaBayes",
+ "limiting.gamma"
)
## function which we want to expose to the user
@@ -58,4 +59,5 @@
export(ql,rql,ml.ql)
export(adaBayes)
-export(rIG, rNIG, rbgamma, rngamma, rstable) ##:: random number generator for Inverse Gaussian
\ No newline at end of file
+export(rIG, rNIG, rbgamma, rngamma, rstable) ##:: random number generator for Inverse Gaussian
+export(limiting.gamma)
Added: pkg/yuima/R/limiting.gamma.R
===================================================================
--- pkg/yuima/R/limiting.gamma.R (rev 0)
+++ pkg/yuima/R/limiting.gamma.R 2009-10-28 11:08:28 UTC (rev 10)
@@ -0,0 +1,231 @@
+setGeneric("limiting.gamma",
+ function(obj, theta, verbose=FALSE)
+ standardGeneric("limiting.gamma")
+ )
+setMethod("limiting.gamma", "yuima",
+ function(obj, theta, verbose=FALSE){
+ limiting.gamma(obj at model, theta, verbose=verbose)
+ })
+setMethod("limiting.gamma", "yuima.model",
+ function(obj, theta, verbose=FALSE){
+
+ ## error check 1
+ if(missing(obj)){
+ stop("main object is missing.")
+ }
+ if(missing(theta)){
+ stop("theta is missing.")
+ }
+ if(is.list(theta)==FALSE){
+ stop("theta required list.\ntheta <- list(theta1, theta2)")
+ }
+
+ if(verbose){
+ cat("initialize ... ")
+ }
+
+ r.size <- obj at noise.number
+ d.size <- obj at equation.number
+ state <- obj at state.variable
+ THETA.1 <- obj at parameter@diffusion
+ THETA.2 <- obj at parameter@drift
+ mi.size <- c(length(THETA.1), length(THETA.2))
+
+ ## error check 2
+ if(d.size!=1){
+ stop("this program is 1-dimention yuima limitation.")
+ }
+ if(mi.size[1]!=length(theta[[1]])){
+ stop("the length of m1 and theta1 is different.")
+ }
+ if(mi.size[2]!=length(theta[[2]])){
+ stop("the length of m2 and theta2 is different.")
+ }
+
+
+ Differentiation.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)
+ }
+
+ Differentiation.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)
+ }
+
+
+ ## assign theta
+ for(i in 1:mi.size[1]){
+ assign(THETA.1[i], theta[[1]][i])
+ }
+ for(i in 1:mi.size[2]){
+ assign(THETA.2[i], theta[[2]][i])
+ }
+ if(verbose){
+ cat("done\n")
+ }
+ ## p(x)
+ if(verbose){
+ cat("get C ... ")
+ }
+
+ ## a part of p(x) function
+ tmp.y <- function(x.arg){
+ assign(obj at state.variable, x.arg)
+ a.eval <- eval(obj at drift)
+ b.eval <- numeric(d.size)
+ for(i in 1:d.size){
+ b.eval[i] <- eval(obj at diffusion[[1]][i])
+ }
+ return(2*a.eval/as.double(t(b.eval)%*%b.eval))
+ }
+
+ ## a part of p(x) function for integrate
+ integrate.tmp.y <- function(x.arg){
+ tmp <- numeric(length(x.arg))
+ for(i in 1:length(x.arg)){
+ tmp[i]<-tmp.y(x.arg[i])
+ }
+ return(tmp)
+ }
+
+ ## p(x) function without normalize const C
+ p0.x <- function(x.arg){
+ assign(obj at state.variable, x.arg)
+ b.eval <- numeric(d.size)
+ for(i in 1:d.size){
+ b.eval[i] <- eval(obj at diffusion[[1]][i])
+ }
+ return(exp(as.double(integrate(integrate.tmp.y, 0, x.arg)$value))/as.double(t(b.eval)%*%b.eval))
+ }
+
+ ## p0.x function for integrate
+ integrate.p0.x <- function(x.arg){
+ tmp <- numeric(length(x.arg))
+ for(i in 1:length(x.arg)){
+ tmp[i]<-p0.x(x.arg[i])
+ }
+ return(tmp)
+ }
+
+ ## normalize const
+ C <- 1/integrate(integrate.p0.x, -Inf, Inf)$value
+ if(verbose){
+ cat("done\n")
+ }
+
+ ## gamma1
+ if(verbose){
+ cat("get gamma1 ... ")
+ }
+ counter.gamma1 <- 1
+
+ ## gamma1 function
+ get.gamma1 <- function(x.arg){
+ assign(obj at state.variable, x.arg)
+ b.eval <- numeric(d.size)
+ for(i in 1:d.size){
+ b.eval[i] <- eval(obj at diffusion[[1]][i])
+ }
+ Bi <- solve(t(b.eval)%*%b.eval)
+ dTHETA.1.B <- array(0,c(d.size, d.size, mi.size[1]))
+ tmp1.B <- array(0, c(d.size, d.size, mi.size[1]))
+ tmp2.B <- numeric(mi.size[1])
+ for(k in 1:mi.size[1]){
+ dTHETA.1.b <- Differentiation.scalar(obj at diffusion[[1]], THETA.1[k], d.size)
+ dTHETA.1.b.eval <- numeric(d.size)
+ for(i in 1:d.size){
+ dTHETA.1.b.eval[i] <- eval(dTHETA.1.b[i])
+ }
+ tmp <- b.eval%*%t(dTHETA.1.b.eval)
+ dTHETA.1.B[, , k] <- tmp+t(tmp)
+ tmp1.B[, , k] <- dTHETA.1.B[, , k]%*%Bi
+ tmp2.B[k] <- tmp1.B[, , k] #sum(diag(tmp1.B[,,k])) 1-dimention limitation
+ }
+ tmp <- tmp2.B %*% t(tmp2.B) * p0.x(x.arg)
+ return(tmp[((counter.gamma1-1)%%mi.size[1])+1, ((counter.gamma1-1)%/%mi.size[1])+1])
+ }
+
+
+ ## gamma1 function for integrate
+ integrate.get.gamma1 <- function(x.arg){
+ tmp <- numeric(length(x.arg))
+ for(i in 1:length(x.arg)){
+ tmp[i] <- get.gamma1(x.arg[i])
+ }
+ return(tmp*C/2)
+ }
+
+ ## calculating gamma1
+ gamma1 <- matrix(0, mi.size[1], mi.size[1])
+ for(i in 1:(mi.size[1]*mi.size[1])){
+ gamma1[((i-1)%%mi.size[1])+1,((i-1)%/%mi.size[1])+1] <- integrate(integrate.get.gamma1, -Inf, Inf)$value
+ counter.gamma1 <- counter.gamma1+1
+ }
+ if(verbose){
+ cat("done\n")
+ }
+
+ ## gamma2
+ if(verbose){
+ cat("get gamma2 ... ")
+ }
+ counter.gamma2 <- 1
+
+ ## gamma2 function
+ get.gamma2 <- function(x.arg){
+ assign(obj at state.variable, x.arg)
+ if(mi.size[2]==1){
+ dTHETA.2.a <- Differentiation.scalar(obj at drift, THETA.2, d.size)
+ }else{
+ dTHETA.2.a <- Differentiation.vector(obj at drift, THETA.2, d.size, mi.size[2])
+ }
+ dTHETA.2.a.eval <- numeric(mi.size[2])
+ for(i in 1:mi.size[2]){
+ dTHETA.2.a.eval[i] <- eval(dTHETA.2.a[i])
+ }
+ b.eval <- numeric(d.size)
+ for(i in 1:d.size){
+ b.eval[i] <- eval(obj at diffusion[[1]][i])
+ }
+ Bi <- solve(t(b.eval)%*%b.eval)
+ tmp <- dTHETA.2.a.eval %*% Bi %*% t(dTHETA.2.a.eval) * p0.x(x.arg)
+ return(tmp[((counter.gamma2-1)%%mi.size[2])+1, ((counter.gamma2-1)%/%mi.size[2])+1])
+ }
+
+ ## gamma2 function for intefrate
+ integrate.get.gamma2 <- function(x.arg){
+ tmp <- numeric(length(x.arg))
+ for(i in 1:length(x.arg)){
+ tmp[i]<-get.gamma2(x.arg[i])
+ }
+ return(tmp*C)
+ }
+
+ ## calculating gamma2
+ gamma2 <- matrix(0, mi.size[2], mi.size[2])
+ for(i in 1:(mi.size[2]*mi.size[2])){
+ gamma2[((i-1)%%mi.size[2])+1, ((i-1)%/%mi.size[2])+1] <- integrate(integrate.get.gamma2, -Inf, Inf)$value
+ counter.gamma2 <- counter.gamma2+1
+ }
+ if(verbose){
+ cat("done\n")
+ }
+
+ ## make list for return
+ poi1 <- list(gamma1, gamma2)
+ names(poi1) <- c("gamma1", "gamma2")
+ poi2 <- append(gamma1, gamma2)
+ ret <- list(poi1, poi2)
+ names(ret) <- c("list", "vec")
+
+ return(ret)
+ })
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2009-10-28 10:20:11 UTC (rev 9)
+++ pkg/yuima/R/quasi-likelihood.R 2009-10-28 11:08:28 UTC (rev 10)
@@ -1,39 +1,11 @@
-# quasi-likelihood function
+##::quasi-likelihood function
-# extract diffusion term from yuima
-#para: parameter of diffusion term (theta1)
-"calc.diffusion" <- function(yuima,para){
+##::extract drift term from yuima
+##::para: parameter of drift term (theta2)
+"calc.drift" <- function(yuima, para){
r.size <- yuima at model@noise.number
d.size <- yuima at model@equation.number
modelstate <- yuima at model@state.variable
- modelpara <- yuima at model@parameter at diffusion
- DIFFUSION <- yuima at model@diffusion
- division <- length(yuima)[1]
- diff <- array(0,dim=c(d.size,r.size,division))
- X <- as.matrix(onezoo(yuima))
- for(i in 1:length(modelpara)){
- assign(modelpara[i],para[i])
- }
- for(r in 1:r.size){
- for(t in 1:division){
- Xt <- X[t,]
- for(d in 1:d.size){
- assign(modelstate[d],Xt[d])
- }
- for(d in 1:d.size){
- diff[d,r,t] <- eval(DIFFUSION[[d]][r])
- }
- }
- }
- return(diff)
-}
-
-#extract drift term from yuima
-#para: parameter of drift term (theta2)
-"calc.drift" <- function(yuima,para){
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- modelstate <- yuima at model@state.variable
modelpara <- yuima at model@parameter at drift
DRIFT <- yuima at model@drift
division <- length(yuima)[1]
@@ -54,226 +26,376 @@
return(drift)
}
-#calcurate diffusion%*%t(diffusion) matrix
+
+##::extract diffusion term from yuima
+##::para: parameter of diffusion term (theta1)
+"calc.diffusion" <- function(yuima, para){
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ modelstate <- yuima at model@state.variable
+ modelpara <- yuima at model@parameter at diffusion
+ DIFFUSION <- yuima at model@diffusion
+ division <- length(yuima)[1]
+ diff <- array(0, dim=c(d.size, r.size, division))
+ X <- as.matrix(onezoo(yuima))
+ for(k in 1:length(modelpara)){
+ assign(modelpara[k], para[k])
+ }
+ for(r in 1:r.size){
+ for(t in 1:division){
+ Xt <- X[t, ]
+ for(d in 1:d.size){
+ assign(modelstate[d], Xt[d])
+ }
+ for(d in 1:d.size){
+ diff[d, r, t] <- eval(DIFFUSION[[d]][r])
+ }
+ }
+ }
+ return(diff)
+}
+
+##::calcurate diffusion%*%t(diffusion) matrix
calc.B <- function(diff){
d.size <- dim(diff)[1]
division <- dim(diff)[3]
- B <- array(0,dim=c(d.size,d.size,division))
+ B <- array(0, dim=c(d.size, d.size, division))
for(t in 1:division){
- B[,,t] <- diff[,,t]%*%t(diff[,,t])
+ B[, , t] <- diff[, , t]%*%t(diff[, , t])
}
return(B)
}
-#calculate the log quasi-likelihood with parameters (theta2,theta1) and X.
-##yuima : yuima object
-##theta2 : parameter in drift.
-##theta1 : parameter in diffusion.
-##h : time width.
-setGeneric("ql",
- function(yuima,theta2,theta1,h,print=FALSE)
- standardGeneric("ql")
- )
-setMethod("ql", "ANY",
- function(yuima,theta2,theta1,h,print=FALSE){
- if( missing(yuima)){
- cat("\nyuima object is missing.\n")
- return(NULL)
- }
- if( missing(theta2) || missing(theta1)){
- cat("\nparameters of yuima.model is missing.\n")
- return(NULL)
- }
- if( missing(h)){
- cat("\nlength of each time is missing.\n")
- return(NULL)
- }
- if(length(yuima at model@parameter at drift)!=length(theta2)){
- cat("\nlength of drift parameter is strange.\n")
- return(NULL)
- }
- if(length(yuima at model@parameter at diffusion)!=length(theta1)){
- cat("\nlength of diffusion parameter is strange.\n")
- return(NULL)
- }
-
-
+calc.B.grad <- function(yuima, para){
+ r.size <- yuima at model@noise.number
d.size <- yuima at model@equation.number
+ modelstate <- yuima at model@state.variable
+ modelpara <- yuima at model@parameter at diffusion
+ DIFFUSION <- yuima at model@diffusion
division <- length(yuima)[1]
X <- as.matrix(onezoo(yuima))
- deltaX <- matrix(0,division-1,d.size)
- for(t in 1:(division-1)){
- deltaX[t,] <- X[t+1,]-X[t,]
+
+ for(k in 1:length(modelpara)){
+ assign(modelpara[k], para[k])
}
- if(is.nan(theta2) || is.nan(theta1)){
- stop("error: theta is not a namber in parameter matrix")
- }
- drift <- calc.drift(yuima,para=theta2)
- diff <- calc.diffusion(yuima,para=theta1)
- B <- calc.B(diff)
+ ## B <- list(NULL)
+ ## for(d in 1:d.size){
+ ## B[[d]] <- list(NULL)
+ ## for(d2 in 1:d.size){
+ ## if(d2<d){
+ ## B[[d]][[d2]] <- B[[d2]][[d]]
+ ## }else{
+ ## B[[d]][[d2]] <- expression(0)
+ ## for(r in 1:r.size){
+ ## B[[d]][[d2]] <- parse(text=paste(as.character(B[[d]][[d2]]), "+(", as.character(DIFFUSION[[d]][r]), ")*(", as.character(DIFFUSION[[d2]][r]), ")"))
+ ## }
+ ## }
+ ## }
+ ## }
- QL <- 0
- pn <- numeric(division-1)
- for(t in 1:(division-1)){
- if(det(as.matrix(B[,,t]))==0){
- pn[t] <- log(1)
- }else{
- pn[t] <- log(1/((2*pi*h)^(d.size/2)*det(as.matrix(B[,,t]))^(1/2)) *
- exp((-1/(2*h))*t(deltaX[t,]-h*drift[t,])%*%solve(as.matrix(B[,,t]))%*%(deltaX[t,]-h*drift[t,])))
- QL <- QL+pn[t]
- if(pn[t]==-Inf && FALSE){
- cat("t:",t, "\n")
- cat("B[,,t]:",B[,,t], "\n")
- cat("det(B):",det(as.matrix(B[,,t])),"\n")
- cat("deltaX[t,]", deltaX[t,], "\n")
- cat("drift[t,]", drift[t,], "\n")
+ B.grad <- array(0, dim=c(d.size, d.size, division, length(modelpara)))
+ for(k in 1:length(modelpara)){
+ for(t in 1:division){
+ for(d in 1:d.size){
+ assign(modelstate[d], X[t, d])
}
+ for(d in 1:d.size){
+ for(d2 in 1:d.size){
+ if(d2<d){
+ B.grad[d, d2, t, k] <- B.grad[d2, d, t, k]
+ }else{
+ B <- expression(0)
+ for(r in 1:r.size){
+ B <- parse(text=paste(as.character(B), "+(", as.character(DIFFUSION[[d]][r]), ")*(", as.character(DIFFUSION[[d2]][r]), ")"))
+ }
+ B.grad[d, d2, t, k] <- eval(D(B, modelpara[k]))
+ }
+ }
+ }
}
}
- if(QL==-Inf){
- warning("quasi likelihood is too small too calculate.")
- }
- if(print==TRUE){
- print(paste("QL:",QL," theta2:",theta2," theta1:",theta1))
- }
- return(QL)
-})
+ return(B.grad)
+}
+##::calculate the log quasi-likelihood with parameters (theta2, theta1) and X.
+##::yuima : yuima object
+##::theta2 : parameter in drift.
+##::theta1 : parameter in diffusion.
+##::h : time width.
+setGeneric("ql",
+ function(yuima, theta2, theta1, h, print=FALSE)
+ standardGeneric("ql")
+ )
+setMethod("ql", "yuima",
+ function(yuima, theta2, theta1, h, print=FALSE){
+ ##QLG <- ql.grad(yuima, theta2, theta1, h, print=FALSE)
+ ##print(QLG)
+
+ if( missing(yuima)){
+ cat("\nyuima object is missing.\n")
+ return(NULL)
+ }
+ if( missing(theta2) || missing(theta1)){
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
+ }
+ if( missing(h)){
+ cat("\nlength of each time is missing.\n")
+ return(NULL)
+ }
+ if(length(yuima at model@parameter at drift)!=length(theta2)){
+ cat("\nlength of drift parameter is strange.\n")
+ return(NULL)
+ }
+ if(length(yuima at model@parameter at diffusion)!=length(theta1)){
+ cat("\nlength of diffusion parameter is strange.\n")
+ return(NULL)
+ }
+
+
+ d.size <- yuima at model@equation.number
+ division <- length(yuima)[1]
+ X <- as.matrix(onezoo(yuima))
+ deltaX <- matrix(0, division-1, d.size)
+ for(t in 1:(division-1)){
+ deltaX[t, ] <- X[t+1, ]-X[t, ]
+ }
+ if(is.nan(theta2) || is.nan(theta1)){
+ stop("error: theta is not a namber in parameter matrix")
+ }
+ drift <- calc.drift(yuima, para=theta2)
+ diff <- calc.diffusion(yuima, para=theta1)
+ B <- calc.B(diff)
+
+ QL <- 0
+ pn <- numeric(division-1)
+ for(t in 1:(division-1)){
+ if(det(as.matrix(B[, , t]))==0){
+ pn[t] <- log(1)
+ }else{
+ pn[t] <- log( 1/((2*pi*h)^(d.size/2)*det(as.matrix(B[, , t]))^(1/2)) *
+ exp((-1/(2*h))*t(deltaX[t, ]-h*drift[t, ])%*%solve(as.matrix(B[, , t]))%*%(deltaX[t,]-h*drift[t, ])) )
+ QL <- QL+pn[t]
+ if(pn[t]==-Inf && FALSE){
+ cat("t:", t, "\n")
+ cat("B[, , t]:", B[, , t], "\n")
+ cat("det(B):",det(as.matrix(B[, , t])), "\n")
+ cat("deltaX[t, ]", deltaX[t, ], "\n")
+ cat("drift[t, ]", drift[t, ], "\n")
+ }
+ }
+ }
+ if(QL==-Inf){
+ warning("quasi likelihood is too small to calculate.")
+ }
+ if(print==TRUE){
+ print(paste("QL:", QL, " theta2:", theta2, " theta1:", theta1))
+ }
+ return(QL)
+ })
-#calculate the relative log quasi-likelihood with parameters (theta2,theta1) and X.
-##yuima : yuima object
-##theta2 : parameter in drift.
-##theta1 : parameter in diffusion.
-##ptheta2 : parameter in drift in the prevous mcmc step.
-##ptheta1 : parameter in diffusion in the prevous mcmc step.
-##h : time width.
-setGeneric("rql",
- function(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE)
- standardGeneric("rql")
- )
-setMethod("rql", "ANY" ,
- function(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE){
+
+ql.grad <- function(yuima, theta2, theta1, h, print=FALSE){
if( missing(yuima)){
- cat("\nyuima object is missing.\n")
- return(NULL)
+ cat("\nyuima object is missing.\n")
+ return(NULL)
}
if( missing(theta2) || missing(theta1)){
- cat("\nparameters of yuima.model is missing.\n")
- return(NULL)
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
}
if( missing(h)){
- cat("\nlength of each time is missing.\n")
- return(NULL)
+ cat("\nlength of each time is missing.\n")
+ return(NULL)
}
if(length(yuima at model@parameter at drift)!=length(theta2)){
- cat("\nlength of drift parameter is strange.\n")
- return(NULL)
+ cat("\nlength of drift parameter is strange.\n")
+ return(NULL)
}
if(length(yuima at model@parameter at diffusion)!=length(theta1)){
- cat("\nlength of diffusion parameter is strange.\n")
- return(NULL)
+ cat("\nlength of diffusion parameter is strange.\n")
+ return(NULL)
}
-
d.size <- yuima at model@equation.number
division <- length(yuima)[1]
X <- as.matrix(onezoo(yuima))
- deltaX <- matrix(0,division-1,d.size)
+ deltaX <- matrix(0, division-1, d.size)
for(t in 1:(division-1)){
- deltaX[t,] <- X[t+1,]-X[t,]
+ deltaX[t, ] <- X[t+1, ]-X[t, ]
}
if(is.nan(theta2) || is.nan(theta1)){
stop("error: theta is not a namber in parameter matrix")
}
- drift <- calc.drift(yuima,para=theta2)
- diff <- calc.diffusion(yuima,para=theta1)
+ drift <- calc.drift(yuima, para=theta2)
+ diff <- calc.diffusion(yuima, para=theta1)
B <- calc.B(diff)
-
- pdrift <- calc.drift(yuima,para=ptheta2)
- pdiff <- calc.diffusion(yuima,para=ptheta1)
- pB <- calc.B(pdiff)
+ B.grad <- calc.B.grad(yuima, para=theta1)
- rQL <- 0
- pn <- numeric(division-1)
- for(t in 1:(division-1)){
- if(det(as.matrix(B[,,t]))*det(as.matrix(pB[,,t]))==0){
- pn[t] <- log(1)
- }else{
- pn[t] <- -log(det(as.matrix(B[,,t]))^(1/2))-1/(2*h)*t(deltaX[t,]-h*drift[t,])%*%solve(as.matrix(B[,,t]))%*%(deltaX[t,]-h*drift[t,])
- pn[t] <- pn[t]-(-log(det(as.matrix(pB[,,t]))^(1/2))-1/(2*h)*t(deltaX[t,]-h*pdrift[t,])%*%solve(as.matrix(pB[,,t]))%*%(deltaX[t,]-h*pdrift[t,]))
- rQL <- rQL+pn[t]
+ QLG <- numeric(dim(B.grad)[4])
+ for(k in 1:length(QLG)){
+ pg1 <- 0
+ pg2 <- 0
+ for(t in 1:(division-1)){
+ B.tmp <- as.matrix(B[, , t])
+ if(det(B.tmp)!=0){
+ B.grad.tmp <- as.matrix(B.grad[, , t, k])
+ aa <- as.matrix(deltaX[t, ]-h*drift[t, ])
+ pg1 <- pg1 + sum(diag(solve(B.tmp)%*%B.grad.tmp))
+ pg2 <- pg2 + t(aa)%*%solve(B.tmp)%*%B.grad.tmp%*%solve(B.tmp)%*%aa
+ }
}
+ QLG[k] <- (-1/2)*pg1 + 1/(2*h)*pg2
+
+ if(QLG[k]==-Inf){
+ warning(paste("gradient[", k, "] of quasi likelihood is too small too calculate.", sep=""))
+ }
}
- if(print==TRUE){
- print(paste("relative QL:",rQL," theta2:",theta2," theta1:",theta1))
+ QLG <- QLG/sqrt(sum(QLG^2))
+ if(print){
+ print(paste("QLG:", QLG, " theta2:", theta2, " theta1:", theta1))
}
- return(rQL)
-})
+ return(QLG)
+}
-#estimate parameters(theta2,theta1) with a constraint ui%*%theta-ci=0
-##yuima : yuima object
-##theta2 : init parameter in drift term.
-##theta1 : init parameter in diffusion term.
-##h : length between each observation time.
-##theta1.lim, theta2.lim : limit of those parameters.
-###example: 0 <= theta1 <= 1 theta1.lim = matrix(c(0,1),1,2)
-###if theta1, theta2 are matrix, theta.lim can be defined matrix like rbind(c(0,1),c(0,1))
+##::calculate the relative log quasi-likelihood with parameters (theta2, theta1) and X.
+##::yuima : yuima object
+##::theta2 : parameter in drift.
+##::theta1 : parameter in diffusion.
+##::ptheta2 : parameter in drift in the prevous mcmc step.
+##::ptheta1 : parameter in diffusion in the prevous mcmc step.
+##::h : time width.
+setGeneric("rql",
+ function(yuima, theta2, theta1, ptheta2, ptheta1, h, print=FALSE)
+ standardGeneric("rql")
+ )
+setMethod("rql", "yuima",
+ function(yuima, theta2, theta1, ptheta2, ptheta1, h, print=FALSE){
+ if(missing(yuima)){
+ cat("\nyuima object is missing.\n")
+ return(NULL)
+ }
+ if(missing(theta2) || missing(theta1)){
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
+ }
+ if(missing(h)){
+ cat("\nlength of each time is missing.\n")
+ return(NULL)
+ }
+ if(length(yuima at model@parameter at drift)!=length(theta2)){
+ cat("\nlength of drift parameter is strange.\n")
+ return(NULL)
+ }
+ if(length(yuima at model@parameter at diffusion)!=length(theta1)){
+ cat("\nlength of diffusion parameter is strange.\n")
+ return(NULL)
+ }
+
+ d.size <- yuima at model@equation.number
+ division <- length(yuima)[1]
+ X <- as.matrix(onezoo(yuima))
+ deltaX <- matrix(0, division-1, d.size)
+ for(t in 1:(division-1)){
+ deltaX[t, ] <- X[t+1, ]-X[t, ]
+ }
+ if(is.nan(theta2) || is.nan(theta1)){
+ stop("error: theta is not a namber in parameter matrix")
+ }
+ drift <- calc.drift(yuima, para=theta2)
+ diff <- calc.diffusion(yuima, para=theta1)
+ B <- calc.B(diff)
+
+ pdrift <- calc.drift(yuima, para=ptheta2)
+ pdiff <- calc.diffusion(yuima, para=ptheta1)
+ pB <- calc.B(pdiff)
+
+ rQL <- 0
+ pn <- numeric(division-1)
+ for(t in 1:(division-1)){
+ if(det(as.matrix(B[, , t]))*det(as.matrix(pB[, , t]))==0){
+ pn[t] <- log(1)
+ }else{
+ pn[t] <- -log(det(as.matrix(B[, , t]))^(1/2))-1/(2*h)*t(deltaX[t, ]-h*drift[t, ])%*%solve(as.matrix(B[, , t]))%*%(deltaX[t, ]-h*drift[t, ])
+ pn[t] <- pn[t]-(-log(det(as.matrix(pB[, , t]))^(1/2))-1/(2*h)*t(deltaX[t, ]-h*pdrift[t, ])%*%solve(as.matrix(pB[, , t]))%*%(deltaX[t, ]-h*pdrift[t, ]))
+ rQL <- rQL+pn[t]
+ }
+ }
+ if(print==TRUE){
+ print(paste("relative QL:", rQL, " theta2:", theta2, " theta1:", theta1))
+ }
+ return(rQL)
+ })
+
+##::estimate parameters(theta2,theta1) with a constraint ui%*%theta-ci=0
+##::yuima : yuima object
+##::theta2 : init parameter in drift term.
+##::theta1 : init parameter in diffusion term.
+##::h : length between each observation time.
+##::theta1.lim, theta2.lim : limit of those parameters.
+##::example: 0 <= theta1 <= 1 theta1.lim = matrix(c(0,1),1,2)
+##::if theta1, theta2 are matrix, theta.lim can be defined matrix like rbind(c(0,1),c(0,1))
setGeneric("ml.ql",
- function(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE)
+ function(yuima, theta2, theta1, h, theta2.lim=matrix(c(0, 1), 1, 2), theta1.lim=matrix(c(0, 1), 1, 2), print=FALSE, BFGS=FALSE)
standardGeneric("ml.ql")
)
-setMethod("ml.ql", "ANY" ,
- function(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE){
- if( missing(yuima)){
- cat("\nyuima object is missing.\n")
- return(NULL)
- }
- if( missing(theta2) || missing(theta1)){
- cat("\nparameters of yuima.model is missing.\n")
- return(NULL)
- }
- if(length(yuima at model@parameter at drift)!=length(theta2)){
- cat("\nlength of drift parameter is strange.\n")
- return(NULL)
- }
- if(length(yuima at model@parameter at diffusion)!=length(theta1)){
- cat("\nlength of diffusion parameter is strange.\n")
- return(NULL)
- }
- if( missing(h)){
- cat("\nlength of each time is missing.\n")
- return(NULL)
- }
- if(is.matrix(theta2.lim) && is.matrix(theta1.lim)){
- if(ncol(theta2.lim)!=2 || ncol(theta1.lim)!=2){
- cat("\ntheta.lim is not available.\n")
- return(NULL)
- }
- }
+setMethod("ml.ql", "yuima",
+ function(yuima, theta2, theta1, h, theta2.lim=matrix(c(0, 1), 1, 2), theta1.lim=matrix(c(0, 1), 1, 2), print=FALSE, BFGS=FALSE){
+ if( missing(yuima)){
+ cat("\nyuima object is missing.\n")
+ return(NULL)
+ }
+ if( missing(theta2) || missing(theta1)){
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
+ }
+ if(length(yuima at model@parameter at drift)!=length(theta2)){
+ cat("\nlength of drift parameter is strange.\n")
+ return(NULL)
+ }
+ if(length(yuima at model@parameter at diffusion)!=length(theta1)){
+ cat("\nlength of diffusion parameter is strange.\n")
+ return(NULL)
+ }
+ if( missing(h)){
+ cat("\nlength of each time is missing.\n")
+ return(NULL)
+ }
+ if(is.matrix(theta2.lim) && is.matrix(theta1.lim)){
+ if(ncol(theta2.lim)!=2 || ncol(theta1.lim)!=2){
+ cat("\ntheta.lim is not available.\n")
+ return(NULL)
+ }
+ }
+ if( length(theta2)!=1 && length(theta2)!=nrow(theta2.lim)){
+ cat("\nsize of theta2 and theta2.lim are different.\n")
+ return(NULL)
+ }
+ if( length(theta1)!=1 && length(theta1)!=nrow(theta1.lim)){
+ cat("\nsize of theta1 and theta1.lim are different.\n")
+ return(NULL)
+ }
+
+ ui <- rbind(diag(length(theta1)+length(theta2)), (-1)*diag(length(theta1)+length(theta2)))
+ ci <- c(rbind(theta2.lim, theta1.lim))
+ ci[(length(ci)/2+1):length(ci)] <- (-1)*ci[(length(ci)/2+1):length(ci)]
+
+ ql.opt <- function(theta=c(theta2, theta1)){
+ return(ql(yuima, theta2=theta[1:length(theta2)], theta1=theta[(length(theta2)+1):length(theta)], h=h, print=print))
+ }
+ ql.grad.opt <- function(theta=c(theta2, theta1)){
+ return(ql.grad(yuima, theta2=theta[1:length(theta2)], theta1=theta[(length(theta2)+1):length(theta)], h=h, print=print))
+ }
+
+ if(BFGS){
+ opt <- constrOptim(c(theta2, theta1), ql.opt, ql.grad.opt, ui=ui, ci=ci, control=list(fnscale=-1), method="BFGS", outer.iterations=500)
+ }else{
+ opt <- constrOptim(c(theta2, theta1), ql.opt, NULL, ui=ui, ci=ci, control=list(fnscale=-1), outer.iterations=500)
+ }
+
+ if(opt$convergence != 0){
+ print("WARNING:optimization did not converge.")
+ }
+ return(opt)
+ })
- if( length(theta2)!=1 && length(theta2)!=nrow(theta2.lim)){
- cat("\nsize of theta2 and theta2.lim are different.\n")
- return(NULL)
- }
- if( length(theta1)!=1 && length(theta1)!=nrow(theta1.lim)){
- cat("\nsize of theta1 and theta1.lim are different.\n")
- return(NULL)
- }
-
-
-
- ui <- rbind(diag(length(theta1)+length(theta2)),(-1)*diag(length(theta1)+length(theta2)))
- ci <- c(rbind(theta2.lim,theta1.lim))
- ci[(length(ci)/2+1):length(ci)] <- (-1)*ci[(length(ci)/2+1):length(ci)]
-
- ql.opt <- function(theta=c(theta2,theta1)){
- return(ql(yuima,theta2=theta[1:length(theta2)],theta1=theta[(length(theta2)+1):length(theta)],h=h,print=print))
- }
-
- opt <- constrOptim(c(theta2,theta1),ql.opt,NULL,ui=ui,ci=ci,control=list(fnscale=-1),outer.iterations=500)
- if( opt$convergence != 0){
- print("WARNING:optimization did not converge.")
- }
- return(opt)
-})
-
Added: pkg/yuima/man/limiting.gamma.Rd
===================================================================
--- pkg/yuima/man/limiting.gamma.Rd (rev 0)
+++ pkg/yuima/man/limiting.gamma.Rd 2009-10-28 11:08:28 UTC (rev 10)
@@ -0,0 +1,53 @@
+\name{limiting.gamma}
+\alias{limiting.gamma}
+
+\title{calculate the value of limiting covariance matrices : Gamma}
+\description{To confirm assysmptotic normality of theta estimators.}
+\usage{
+limiting.gamma(obj,theta,verbose=FALSE)
+}
+\arguments{
+ \item{obj}{an yuima or yuima.model object.}
+ \item{theta}{true theta}
+ \item{verbose}{an option for display a verbose process.}
+}
+\details{
+ Calculate the value of limiting covariance matrices Gamma.
+ The returned values gamma1 and gamma2 are used to confirm assysmptotic normality of theta estimators.
+ this program is limitted to 1-dimention-sde model for now.
+}
+\value{
+ \item{gamma1}{a theoretical figure for variance of theta1 estimator}
+ \item{gamma2}{a theoretical figure for variance of theta2 estimator}
+}
+\author{YUIMA Project Team}
+\note{
+ we need to fix this routine.
+}
+\examples{
+
+set.seed(123)
+
+## Yuima
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 10
More information about the Yuima-commits
mailing list