[Yuima-commits] r141 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 3 04:22:06 CET 2011
Author: hinohide
Date: 2011-01-03 04:22:00 +0100 (Mon, 03 Jan 2011)
New Revision: 141
Removed:
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/man/quasi-likelihood.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
Log:
ml.ql and ql are replaced with qmle and quasilogl
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-11-02 13:12:00 UTC (rev 140)
+++ pkg/yuima/DESCRIPTION 2011-01-03 03:22:00 UTC (rev 141)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.180
-Date: 2010-11-02
+Version: 0.1.181
+Date: 2011-01-03
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2010-11-02 13:12:00 UTC (rev 140)
+++ pkg/yuima/NAMESPACE 2011-01-03 03:22:00 UTC (rev 141)
@@ -32,9 +32,9 @@
"poisson.random.sampling",
"get.zoo.data",
"initialize",
- "ql",
- "rql",
- "ml.ql",
+## "ql",
+## "rql",
+## "ml.ql",
"adaBayes",
"limiting.gamma",
"getF",
@@ -72,7 +72,8 @@
export(get.zoo.data)
-export(ql,rql,ml.ql)
+##export(ql,rql,ml.ql)
+##export(rql)
export(adaBayes)
export(rIG, rNIG, rbgamma, rngamma, rstable) ##:: random number generator for Inverse Gaussian
export(limiting.gamma)
@@ -88,7 +89,7 @@
export(Fnorm)
export(asymptotic_term)
-export(LSE)
+##export(LSE)
export(lse)
export(qmle)
Deleted: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2010-11-02 13:12:00 UTC (rev 140)
+++ pkg/yuima/R/quasi-likelihood.R 2011-01-03 03:22:00 UTC (rev 141)
@@ -1,1121 +0,0 @@
-##::quasi-likelihood function
-
-##::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
- n <- length(yuima)[1]
- drift <- matrix(0,n,d.size)
- X <- as.matrix(onezoo(yuima))
- for(i in 1:length(modelpara)){
- assign(modelpara[i],para[i])
- }
- for(t in 1:n){
- Xt <- X[t,]
- for(d in 1:d.size){
- assign(modelstate[d],Xt[d])
- }
- for(d in 1:d.size){
- drift[t,d] <- eval(DRIFT[d])
- }
- }
- return(drift)
-}
-
-
-##::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
- n <- length(yuima)[1]
- diff <- array(0, dim=c(d.size, r.size, n))
- 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:n){
- 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]
- n <- dim(diff)[3]
- B <- array(0, dim=c(d.size, d.size, n))
- for(t in 1:n){
- B[, , t] <- diff[, , t]%*%t(diff[, , t])
- }
- return(B)
-}
-
-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
- n <- length(yuima)[1]
- X <- as.matrix(onezoo(yuima))
-
- for(k in 1:length(modelpara)){
- assign(modelpara[k], para[k])
- }
- ## 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]), ")"))
- ## }
- ## }
- ## }
- ## }
-
- B.grad <- array(0, dim=c(d.size, d.size, n, length(modelpara)))
- for(k in 1:length(modelpara)){
- for(t in 1:n){
- 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]))
- }
- }
- }
- }
- }
- 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, param)
- standardGeneric("ql")
- )
-setMethod("ql", "yuima",
- function(yuima, theta2, theta1, h, print=FALSE, param){
- ##QLG <- ql.grad(yuima, theta2, theta1, h, print=FALSE)
- ##print(QLG)
- if( missing(yuima)){
- yuima.warn("yuima object is missing.")
- return(NULL)
- }
-
- ## param handling
- if( missing(param) ){
- if( missing(theta2) || missing(theta1) ){
- yuima.warn("Parameters of yuima.model are missing.")
- return(NULL)
- }
- }else{
- if( missing(theta2) && missing(theta1) ){
- if( !is.list(param) ){
- yuima.warn("param must be list.")
- return(NULL)
- }
-
- if( length(param)!=2){
- yuima.warn("length of param is strange.")
- return(NULL)
- }
-
- ## get theta2 and theta1 from param
- if( is.null(names(param)) ){
- theta2 <- as.vector(param[[1]])
- theta1 <- as.vector(param[[2]])
- }
- else if( sum( names(param)==c("theta2", "theta1") ) == 2 ){
- theta2 <- as.vector(param[[1]])
- theta1 <- as.vector(param[[2]])
- }
- else if( sum( names(param)==c("theta1", "theta2") ) == 2 ){
- theta2 <- as.vector(param[[2]])
- theta1 <- as.vector(param[[1]])
- }
- else{
- yuima.warn("names of param are strange.")
- return(NULL)
- }
- }else{
- yuima.warn("Conflict in parameter specification method.")
- return(NULL)
- }
- }
- ## END param handling
-
- if( missing(h)){
- yuima.warn("length of each time is missing.")
- return(NULL)
- }
-
- if(length(yuima at model@parameter at drift)!=length(theta2)){
- yuima.warn("length of drift parameter is strange.")
- return(NULL)
- }
-
- if(length(yuima at model@parameter at diffusion)!=length(theta1)){
- yuima.warn("length of diffusion parameter is strange.")
- return(NULL)
- }
-
- d.size <- yuima at model@equation.number
- n <- length(yuima)[1]
- X <- as.matrix(onezoo(yuima))
- deltaX <- matrix(0, n-1, d.size)
- for(t in 1:(n-1))
- deltaX[t, ] <- X[t+1, ]-X[t, ]
-
-#
-# print(system.time(for(t in 1:(n-1))
-# deltaX[t, ] <- X[t+1, ]-X[t, ]
-# ))
-# print(system.time(deltaX1 <- diff(X)))
-# print(deltaX)
-# print(deltaX1)
-
- 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(n-1)
- for(t in 1:(n-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)
- })
-
-
-ql.grad <- function(yuima, theta2, theta1, h, print=FALSE){
- if( missing(yuima)){
- yuima.warn("yuima object is missing.")
- return(NULL)
- }
- if( missing(theta2) || missing(theta1)){
- yuima.warn("parameters of yuima.model are missing.")
- return(NULL)
- }
- if( missing(h)){
- yuima.warn("length of each time is missing.")
- return(NULL)
- }
- if(length(yuima at model@parameter at drift)!=length(theta2)){
- yuima.warn("length of drift parameter is strange.")
- return(NULL)
- }
- if(length(yuima at model@parameter at diffusion)!=length(theta1)){
- yuima.warn("length of diffusion parameter is strange.")
- return(NULL)
- }
-
- d.size <- yuima at model@equation.number
- n <- length(yuima)[1]
- X <- as.matrix(onezoo(yuima))
- deltaX <- matrix(0, n-1, d.size)
- for(t in 1:(n-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)
- B.grad <- calc.B.grad(yuima, para=theta1)
-
- QLG <- numeric(dim(B.grad)[4])
- for(k in 1:length(QLG)){
- pg1 <- 0
- pg2 <- 0
- for(t in 1:(n-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=""))
- }
- }
- QLG <- QLG/sqrt(sum(QLG^2))
- if(print){
- print(paste("QLG:", QLG, " theta2:", theta2, " theta1:", theta1))
- }
- return(QLG)
-}
-
-##::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, param, prevparam)
- standardGeneric("rql")
- )
-setMethod("rql", "yuima",
- function(yuima, theta2, theta1, ptheta2, ptheta1, h, print=FALSE, param, prevparam){
- if(missing(yuima)){
- yuima.warn("yuima object is missing.")
- return(NULL)
- }
-
- ## param handling
- if( missing(param) ){
- if( missing(theta2) || missing(theta1) ){
- yuima.warn("parameters of yuima.model is missing.")
- return(NULL)
- }
- }else{
- if( missing(theta2) && missing(theta1) ){
- if( !is.list(param) ){
- yuima.warn("param must be list.")
- return(NULL)
- }
-
- if( length(param)!=2){
- yuima.warn("length of param is strange.")
- return(NULL)
- }
-
- ## get theta2 and theta1 from param
- if( is.null(names(param)) ){
- theta2 <- as.vector(param[[1]])
- theta1 <- as.vector(param[[2]])
- }
- else if( sum( names(param)==c("theta2", "theta1") ) == 2 ){
- theta2 <- as.vector(param[[1]])
- theta1 <- as.vector(param[[2]])
- }
- else if( sum( names(param)==c("theta1", "theta2") ) == 2 ){
- theta2 <- as.vector(param[[2]])
- theta1 <- as.vector(param[[1]])
- }
- else{
- yuima.warn("names of param are strange.")
- return(NULL)
- }
- }else{
- yuima.warn("Conflict in parameter specification method.")
- return(NULL)
- }
- }
- ## END param handling
-
- ## prevparam handling
- if( missing(prevparam) ){
- if( missing(ptheta2) || missing(ptheta1) ){
- yuima.warn("parameters of yuima.model is missing.")
- return(NULL)
- }
- }else{
- if( missing(ptheta2) && missing(ptheta1) ){
- if( !is.list(prevparam) ){
- yuima.warn("param must be list.")
- return(NULL)
- }
- if( length(prevparam)!=2){
- yuima.warn("length of param is strange.")
- return(NULL)
- }
-
- ## get theta2 and theta1 from param
- if( is.null(names(prevparam)) ){
- ptheta2 <- as.vector(prevparam[[1]])
- ptheta1 <- as.vector(prevparam[[2]])
- }
- else if( sum( names(prevparam)==c("ptheta2", "ptheta1") ) == 2 ){
- ptheta2 <- as.vector(prevparam[[1]])
- ptheta1 <- as.vector(prevparam[[2]])
- }
- else if( sum( names(prevparam)==c("ptheta1", "ptheta2") ) == 2 ){
- ptheta2 <- as.vector(prevparam[[2]])
- ptheta1 <- as.vector(prevparam[[1]])
- }
- else{
- yuima.warn("names of prevparam are strange.")
- return(NULL)
- }
-
- }else{
- yuima.warn("Conflict in parameter specification method.")
- return(NULL)
- }
- }
- ## END prevparam handling
-
- if(missing(h)){
- yuima.warn("length of each time is missing.")
- return(NULL)
- }
- if(length(yuima at model@parameter at drift)!=length(theta2)){
- yuima.warn("length of drift parameter is strange.")
- return(NULL)
- }
- if(length(yuima at model@parameter at diffusion)!=length(theta1)){
- yuima.warn("length of diffusion parameter is strange.")
- return(NULL)
- }
-
- d.size <- yuima at model@equation.number
- n <- length(yuima)[1]
- X <- as.matrix(onezoo(yuima))
- deltaX <- matrix(0, n-1, d.size)
- for(t in 1:(n-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(n-1)
- for(t in 1:(n-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)
- })
-
-
-## ml.ql by newton method.
-newton.ml.ql <- function(yuima, theta2, theta1, h, iteration=1, param.only=FALSE, verbose=FALSE, ...){
-
- ## get param
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- modelstate <- yuima at model@state.variable
- modelpara.drift <- yuima at model@parameter at drift
- modelpara.diff <- yuima at model@parameter at diffusion
-
- ## get expression of functions ##
- ## a
- a <- yuima at model@drift
-
- ## b
- b <- yuima at model@diffusion
-
- ## B = b %*% t(b)
- if(length(b)>=2){
- ## expression matrix calculation
- B <- list()
- for( i in 1:d.size){
- for( j in i:d.size){
- tmp <- NULL
- for(l in 1:r.size){
- B.il <- as.character(b[[i]][l])
- B.jl <- as.character(b[[j]][l])
- if(l==1){
- tmp <- paste(B.il, "*", B.jl)
- }else{
- tmp <- paste(tmp, "+", B.il, "*", B.jl)
- }
- }
- ## update B list
- B[[ (i-1)*d.size + j ]] <- parse(text=tmp)
- if(i!=j) B[[ (j-1)*d.size + i ]] <- parse(text=tmp)
- }
- }
- dim(B) <- c(d.size, d.size)
- }else{
- b <- yuima at model@diffusion[[1]]
- B <- parse(text=paste(as.character(b),
- " * ", as.character(b), sep=""))
- B <- list(B)
- dim(B) <- c(1,1)
- }
-
- ## some func def about B
- deriv.B <- function(B, var=character()){
- B.deriv <- B
- for(i in 1:nrow(B)){
- for( j in 1:ncol(B) ){
- B.deriv[i,j][[1]] <- D(B[i,j][[1]], var)
- }
- }
- return(B.deriv)
- }
- eval.B <- function(B, theta1, theta2, Xt.iMinus1, ...){
- ##assign variable
- for(j in 1:length(Xt.iMinus1)){
- assign(modelstate[j], Xt.iMinus1[j])
- }
- for(j in 1:length(theta1)){
- assign(modelpara.diff[j], theta1[j])
- }
- for(j in 1:length(theta2)){
- assign(modelpara.drift[j], theta2[j])
- }
- B.subs <- matrix(0, nrow(B), ncol(B))
- for( i in 1:nrow(B) ){
- for( j in 1:ncol(B) ){
- B.subs[i,j] <- eval(B[i,j][[1]])
- }
- }
- return(B.subs)
- }
- eval.inv.B <- function(B, theta1, theta2, Xt.iMinus1, ...)
- return(solve(eval.B(B, theta1, theta2, Xt.iMinus1, ...)))
- ## END
-
- ## some func about a
- deriv.a <- function(a, var=character()){
- a.deriv <- NULL
- for(i in 1:length(a)){
- a.deriv <- c(a.deriv, as.expression(D(a[i], var)))
- }
- return(a.deriv)
- }
-
- eval.a <- function(a, theta1, theta2, Xt.iMinus1, ...){
- ##assign variable
- for(j in 1:length(Xt.iMinus1)){
- assign(modelstate[j], Xt.iMinus1[j])
- }
- for(j in 1:length(theta1)){
- assign(modelpara.diff[j], theta1[j])
- }
- for(j in 1:length(theta2)){
- assign(modelpara.drift[j], theta2[j])
- }
- a.subs <- matrix(0, length(a), 1)
- for(i in 1:length(a)){
- a.subs[i,] <- eval(a[i])
- }
- return(a.subs)
- }
- ##END
-
- ## dGi_theta1
- dGi.theta1 <- function(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...){
- ##d_theta1 log(detB)+d_theta1 1/2h*(-)T%*%B^(-1)%*%(-)
- ## 2nd term d_theta1 log(det(B))
- term.2 <- matrix(0, length(theta1), 1)
- for( j in 1:length(theta1)){
- term.2[j,1] <- sum(diag(eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.B(deriv.B(B, modelpara.diff[j]), theta1, theta2, Xt.iMinus1)))
- }
-
- ## 3rd term d_theta1 1/2h *(-)B^(-1)(-)
- term.3 <- matrix(0, length(theta1),1)
- for( j in 1:length(theta1)){
- tmp <- delta.i.X - h * eval.a(a, theta1, theta2, Xt.iMinus1)
- term.3[j,1] <- -1/(2*h) * t(tmp) %*% eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.B(deriv.B(B, modelpara.diff[j]), theta1, theta2, Xt.iMinus1) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*% tmp
- }
-
- ret <- 1/2*term.2+term.3
- return(ret)
- }
-
- ## dH_theta1
- dH.theta1 <- function(theta1, theta2, h, yuima, B, a, ...){
- ret <- matrix(0, length(theta1), 1)
- data <- as.matrix(onezoo(yuima))
- for( i in 1:(nrow(data)-1)){
- ##Xt.iMinus1
- Xt.iMinus1 <- data[i,]
- ##delta.i.X
- delta.i.X <- data[i+1,] - data[i,]
- ##calc
- ret <- ret + dGi.theta1(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...)
- }
- ret <- -ret
- return(ret)
- }
-
-
- ## d2Gi_theta1
- d2Gi.theta1 <- function(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a,...){
- ##2nd term
- term.2 <- matrix(0, length(theta1), length(theta1))
- for(j in 1:length(theta1)){
- for(k in 1:length(theta1)){
- tmp <- -eval.inv.B(B,theta1, theta2, Xt.iMinus1) %*%
- eval.B(deriv.B(B,modelpara.diff[k]), theta1, theta2, Xt.iMinus1) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.B(deriv.B(B,modelpara.diff[j]), theta1, theta2, Xt.iMinus1) +
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.B( deriv.B(deriv.B(B,modelpara.diff[j]), modelpara.diff[k]),
- theta1, theta2, Xt.iMinus1)
- term.2[j,k] <- sum(diag(tmp)) / 2
- }
- }
- ##3rd term
- term.3 <- matrix(0, length(theta1), length(theta1))
- for(j in 1:length(theta1)){
- for(k in 1:length(theta1)){
- tmp <- -2 * eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.B( deriv.B(B, modelpara.diff[k]), theta1, theta2, Xt.iMinus1 ) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.B( deriv.B(B, modelpara.diff[j]), theta1, theta2, Xt.iMinus1) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) +
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.B( deriv.B(deriv.B(B,modelpara.diff[j]), modelpara.diff[k]),
- theta1, theta2, Xt.iMinus1 ) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1)
- tmp2 <- delta.i.X -h * eval.a(a, theta1, theta2, Xt.iMinus1)
- term.3[j,k] <- - 1 / (2*h) * t(tmp2) %*% tmp %*% tmp2
- }
- }
-
- ##ret
- ret <- term.2+term.3
- return(ret)
- }
-
- ## d2H_theta1
- d2H.theta1 <-function(theta1, theta2, h, yuima, B, a, ...){
- ret <- matrix(0, length(theta1), length(theta1))
- data <- as.matrix(onezoo(yuima))
- for(i in 1:(nrow(data)-1)){
- ##Xt.iMinus1
- Xt.iMinus1 <- data[i,]
- ##delta.i.X
- delta.i.X <- data[i+1,] - data[i,]
- ##calc
- ret <- ret + d2Gi.theta1(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a,...)
- }
- ret <- -ret
- return(ret)
- }
-
- ## dGi_theta2
- dGi.theta2 <- function(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...){
- ##calc
- ret <- matrix(0, length(theta2), 1)
- for( j in 1:length(theta2) ){
- ret[j,1] <- -t(delta.i.X - h * eval.a(a,theta1, theta2, Xt.iMinus1)) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.a( deriv.a(a, modelpara.drift[j]), theta1, theta2, Xt.iMinus1)
- }
-
- return(ret)
- }
-
- ## dH_theta2
- dH.theta2 <-function(theta1, theta2, h, yuima, B, a, ...){
- ret <- matrix(0, length(theta2), 1)
- data <- as.matrix(onezoo(yuima))
- for( i in 1:(nrow(data)-1)){
- ##Xt.iMinus1
- Xt.iMinus1 <- data[i,]
- ##delta.i.X
- delta.i.X <- data[i+1,] - data[i,]
- ##calc
- ret <- ret + dGi.theta2(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...)
- #cat("dH.theta2-tmp:", ret, "\n")
- }
- ret <- -ret
- return(ret)
- }
-
- ## d2Gi_theta2
- d2Gi.theta2 <- function(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...){
- ##calc
- ret <- matrix(0, length(theta2), length(theta2))
- for(j in 1:length(theta2)){
- for(k in 1:length(theta2)){
- ret[j,k] <- - t(delta.i.X) %*% eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.a( deriv.a( deriv.a(a,modelpara.drift[j]), modelpara.drift[k]),
- theta1, theta2, Xt.iMinus1) +
- h *
- t(eval.a(deriv.a(a, modelpara.drift[j]), theta1, theta2, Xt.iMinus1)) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.a(deriv.a(a, modelpara.drift[k]), theta1, theta2, Xt.iMinus1) +
- h *
- t(eval.a(a, theta1, theta2, Xt.iMinus1)) %*%
- eval.inv.B(B, theta1, theta2, Xt.iMinus1) %*%
- eval.a( deriv.a(deriv.a(a,modelpara.drift[j]), modelpara.drift[k]),
- theta1, theta2, Xt.iMinus1)
- }
- }
- return(ret)
- }
-
- ## d2H_theta2
- d2H.theta2 <- function(theta1, theta2, h, yuima, B, a, ...){
- ret <- matrix(0, length(theta2), length(theta2))
- data <- as.matrix(onezoo(yuima))
- for(i in 1:(nrow(data)-1)){
- ##Xt.iMinus1
- Xt.iMinus1 <- data[i,]
- ##delta.i.X
- delta.i.X <- data[i+1,] - data[i,]
- ##calc
- ret <- ret + d2Gi.theta2(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...)
- }
- ret <- -ret
- return(ret)
- }
- ## END function define ##
-
- ## newton algorithm main ##
- if(!param.only){
- if(verbose){
- #cat("theta2 init:", theta2, "\n")
- cat("theta1 init:", theta1, "\n")
- cat("theta2 init:", theta2, "\n")
- }
-
- for(ite in 1:iteration){
-
- #dHtheta2 <- dH.theta2(theta1, theta2, h, yuima, B, a, ...)
- #d2Htheta2 <- d2H.theta2(theta1, theta2, h, yuima, B, a, ...)
- #theta2 <- as.vector(theta2 - solve(d2Htheta2) %*% dHtheta2)
-
- dHtheta1 <- dH.theta1(theta1, theta2, h, yuima, B, a, ...)
- d2Htheta1 <- d2H.theta1(theta1, theta2, h, yuima, B, a, ...)
- theta1 <- as.vector(theta1 - solve(d2Htheta1) %*% dHtheta1)
-
- dHtheta2 <- dH.theta2(theta1, theta2, h, yuima, B, a, ...)
- d2Htheta2 <- d2H.theta2(theta1, theta2, h, yuima, B, a, ...)
- theta2 <- as.vector(theta2 - solve(d2Htheta2) %*% dHtheta2)
-
- if(verbose){
- cat("\n## Iteration", ite, "##\n")
- #cat("theta2 new:", theta2, "\n")
- cat("theta1 new:", theta1, "\n")
- cat("theta2 new:", theta2, "\n")
- }
- }
- }
- ## END newtom algorithm ##
-
- ##calc Sigma1, Sigma2, Hessian?##
- Sigma1 <- - solve(d2H.theta1(theta1, theta2, h, yuima, B, a, ...))
- Sigma2 <- - solve(d2H.theta2(theta1, theta2, h, yuima, B, a, ...))
- hessian <- 1 ## Not Implemented yet!
- ##END
-
- ##temp
- return(list(theta1.new=theta1, theta2.new=theta2, Sigma1=Sigma1, Sigma2=Sigma2, hessian=hessian))
-
-}
-## END ml.ql by newton method
-
-##::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,
- method="optim",
- param, interval)
- standardGeneric("ml.ql")
- )
-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,
- method="optim", #(BFGS, Newton)
- param, interval){
- if( missing(yuima)){
- yuima.warn("yuima object is missing.")
- return(NULL)
- }
-
- ## param handling
- if( missing(param) ){
- if( missing(theta2) || missing(theta1) ){
- yuima.warn("parameters of yuima.model is missing.")
- return(NULL)
- }
- }else{
- if( missing(theta2) && missing(theta1) ){
- if( !is.list(param) ){
- yuima.warn("param must be list.")
- return(NULL)
- }
-
- if( length(param)!=2){
- yuima.warn("length of param is strange.")
- return(NULL)
- }
-
- ## get theta2 and theta1 from param
- if( is.null(names(param)) ){
- theta2 <- as.vector(param[[1]])
- theta1 <- as.vector(param[[2]])
- }
- else if( sum( names(param)==c("theta2", "theta1") ) == 2 ){
- theta2 <- as.vector(param[[1]])
- theta1 <- as.vector(param[[2]])
- }
- else if( sum( names(param)==c("theta1", "theta2") ) == 2 ){
- theta2 <- as.vector(param[[2]])
- theta1 <- as.vector(param[[1]])
- }
- else{
- yuima.warn("names of param are strange.")
- return(NULL)
- }
- }else{
- yuima.warn("Conflict in parameter specification method.")
- return(NULL)
- }
- }
- ## END param handling
-
-
- if(length(yuima at model@parameter at drift)!=length(theta2)){
- yuima.warn("length of drift parameter is strange.")
- return(NULL)
- }
- if(length(yuima at model@parameter at diffusion)!=length(theta1)){
- yuima.warn("length of diffusion parameter is strange.")
- return(NULL)
- }
- if( missing(h)){
- yuima.warn("length of each time is missing.")
- return(NULL)
- }
-
- ## interval handling
- if( !missing(interval) ){
- if( missing(theta2.lim) && missing(theta2.lim) ){
- if( !is.list(interval) ){
- yuima.warn("interval must be list.")
- return(NULL)
- }
-
- theta2.len <- length(yuima at model@parameter at drift)
- theta1.len <- length(yuima at model@parameter at diffusion)
-
- if( length(interval) != (theta2.len+theta1.len) ){
- yuima.warn("length of interval is strange.")
- return(NULL)
- }
-
- ## get theta2.lim and theta1.lim from interval
- theta2.lim <- NULL
- theta1.lim <- NULL
- for(i in 1:theta2.len){
- theta2.lim <- rbind(theta2.lim, interval[[i]])
- }
- for(i in 1:theta1.len){
- theta1.lim <- rbind(theta1.lim, interval[[theta2.len+i]])
- }
- }else{
- yuima.warn("Conflict in parameter specification method.")
- return(NULL)
- }
- }
- ## END interval handling
-
- 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(method=="Newton"){
- opt <- newton.ml.ql(yuima, theta2, theta1, h,
- iteration=10, param.only=FALSE, verbose=print)
-
- coef <- c(opt$theta2.new, opt$theta1.new)
- min <- ql(yuima, opt$theta2.new, opt$theta1.new, h, print, param)
-
- }else{ ## optim
- 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")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 141
More information about the Yuima-commits
mailing list