[Yuima-commits] r67 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 20 03:17:52 CET 2010
Author: hinohide
Date: 2010-03-20 03:17:52 +0100 (Sat, 20 Mar 2010)
New Revision: 67
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/AllClasses.R
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/R/yuima.data.R
pkg/yuima/man/adaBayes.Rd
pkg/yuima/man/quasi-likelihood.Rd
pkg/yuima/man/setData.Rd
pkg/yuima/man/yuima-class.Rd
pkg/yuima/man/yuima.data-class.Rd
Log:
modified ml.ql
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/DESCRIPTION 2010-03-20 02:17:52 UTC (rev 67)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 0.0.85
-Date: 2010-01-25
+Version: 0.0.83
+Date: 2010-03-20
Depends: methods, zoo, adapt, stats4
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/NAMESPACE 2010-03-20 02:17:52 UTC (rev 67)
@@ -3,6 +3,7 @@
##importFrom("stats", "end", "time", "start")
importFrom("graphics", "plot")
importFrom("zoo")
+importFrom(stats, confint)
exportClasses("yuima",
"yuima.data",
@@ -23,6 +24,7 @@
"cce",
"poisson.random.sampling",
"get.zoo.data",
+ "get.mle",
"initialize",
"ql",
"rql",
@@ -38,7 +40,10 @@
"simFunctional",
"F0",
"Fnorm",
- "asymptotic_term"
+ "asymptotic_term",
+
+ "get.mle"
+
)
## function which we want to expose to the user
@@ -63,6 +68,7 @@
export(poisson.random.sampling)
export(get.zoo.data)
+export(get.mle)
export(ql,rql,ml.ql)
export(adaBayes)
@@ -80,3 +86,4 @@
export(Fnorm)
export(asymptotic_term)
+
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/R/AllClasses.R 2010-03-20 02:17:52 UTC (rev 67)
@@ -42,7 +42,8 @@
# classes
setClass("yuima.data", representation(original.data = "ANY",
- zoo.data = "ANY"
+ zoo.data = "ANY",
+ mle="ANY"
)
)
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/R/quasi-likelihood.R 2010-03-20 02:17:52 UTC (rev 67)
@@ -455,6 +455,319 @@
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")
+ }
+
+ 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)
+
+ if(verbose){
+ cat("\n## Iteration", ite, "##\n")
+ cat("theta2 new:", theta2, "\n")
+ cat("theta1 new:", theta1, "\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.
@@ -464,11 +777,17 @@
##::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, BFGS=FALSE, param, interval)
+ 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, BFGS=FALSE, param, interval){
+ 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)){
cat("\nyuima object is missing.\n")
return(NULL)
@@ -562,25 +881,6 @@
}
## END interval handling
- 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))
@@ -589,157 +889,106 @@
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)
- ## constrOptim2 ( original code is constrOptim() in Package stat )
- constrOptim2 <- function (theta, f, grad, ui, ci, mu = 1e-04, control = list(),
- method = if (is.null(grad)) "Nelder-Mead" else "BFGS",
- outer.iterations = 100, outer.eps = 1e-05, ...)
- {
- if (!is.null(control$fnscale) && control$fnscale < 0)
- mu <- -mu
- R <- function(theta, theta.old, ...) {
- ui.theta <- ui %*% theta
- gi <- ui.theta - ci
- if (any(gi < 0))
- return(NaN)
- gi.old <- ui %*% theta.old - ci
- bar <- sum(gi.old * log(gi) - ui.theta)
- if (!is.finite(bar))
- bar <- -Inf
- f(theta, ...) - mu * bar
+ coef <- c(opt$theta2.new, opt$theta1.new)
+ min <- ql(yuima, opt$theta2.new, opt$theta1.new, h, print, param)
+
+ #}else if(method=="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{ ## 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")
+ return(NULL)
}
- dR <- function(theta, theta.old, ...) {
- ui.theta <- ui %*% theta
- gi <- drop(ui.theta - ci)
- gi.old <- drop(ui %*% theta.old - ci)
- dbar <- colSums(ui * gi.old/gi - ui)
- grad(theta, ...) - mu * dbar
- }
- if (any(ui %*% theta - ci <= 0))
- stop("initial value not feasible")
- obj <- f(theta, ...)
- r <- R(theta, theta, ...)
- for (i in 1L:outer.iterations) {
- obj.old <- obj
- r.old <- r
- theta.old <- theta
- fun <- function(theta, ...) {
- R(theta, theta.old, ...)
- }
- if(!is.null(grad)){
- gradient <- function(theta, ...) {
- dR(theta, theta.old, ...)
- }
- }else{
- gradient <- 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)]
+
+
+ opt <- constrOptim(c(theta2, theta1), ql.opt, NULL, ui=ui,
+ ci=ci, control=list(fnscale=-1), outer.iterations=500)
+ coef <- opt$par
+ min <- opt$value
- ## if hessian calculation is failure, hessian TRUE -> FALSE, and redo.
- a <- try( optim(theta.old, fun, gradient, control = control,
- method = method, hessian=TRUE, ...), silent=TRUE )
- if(class(a) == "try-error"){
- a <- optim(theta.old, fun, gradient, control=control,
- method=method, hessian=FALSE, ...)
- }
- ## END
-
- r <- a$value
- if (is.finite(r) && is.finite(r.old) && abs(r - r.old)/(outer.eps +
- abs(r - r.old)) < outer.eps)
- break
- theta <- a$par
- obj <- f(theta, ...)
- if (obj > obj.old)
- break
- }
- if (i == outer.iterations) {
- a$convergence <- 7
- a$message <- "Barrier algorithm ran out of iterations and did not converge"
- }
- if (mu > 0 && obj > obj.old) {
- a$convergence <- 11
- a$message <- paste("Objective function increased at outer iteration", i)
- }
- if (mu < 0 && obj < obj.old) {
- a$convergence <- 11
- a$message <- paste("Objective function decreased at outer iteration", i)
- }
- a$outer.iterations <- i
- a$barrier.value <- a$value
- a$value <- f(a$par, ...)
- a$barrier.value <- a$barrier.value - a$value
- a
+ theta2 <- coef[1:length(yuima at model@parameter at drift)]
+ theta1 <- coef[(length(yuima at model@parameter at drift)+1):length(coef)]
+ opt2 <- newton.ml.ql(yuima, theta2, theta1, h,
+ iteration=0, param.only=TRUE, verbose=print)
+ opt$Sigma1 <- opt2$Sigma1
+ opt$Sigma2 <- opt2$Sigma2
+ opt$hessian <- opt2$hessian
+
+ if(opt$convergence != 0){
+ print("WARNING:optimization did not converge.")
}
- ## END constrOptim2 ( original code is constrOptim() in Package stat )
-
- ## optimize
- if(BFGS){
- opt <- constrOptim2(c(theta2, theta1), ql.opt, ql.grad.opt, ui=ui, ci=ci,
- control=list(fnscale=-1), method="BFGS", outer.iterations=500)
- }else{
- opt <- constrOptim2(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.")
- }
- ## END optimize
-
- ## opt is converted into the mle form
- ## get call()
+ ##convert to mle object
call <- match.call()
+ names(coef) <- yuima at model@parameter at all
+ fullcoef <- coef
+ method <- method
+
+ hessian <- opt$hessian
+ vcov <- if(length(coef)) solve(hessian)
+ else matrix(numeric(0L), 0L, 0L)
- ## set par
- coef <- opt$par
+ opt <- new("mle", call=call, coef=coef, fullcoef=fullcoef,
+ vcov=vcov, min=min, details=opt, minuslogl=ql.opt,
+ method=method)
+ ##END convert to mle
+
+ yuima at data@mle <- opt
- ## get vcov
- if( !is.null(opt$hessian) ){
- if(length(coef))
- vcov <- solve(opt$hessian)
- else
- vcov <- matrix(numeric(0L), 0L, 0L)
- }else{
- opt$hessian <- "calculation error"
- vcov <- matrix()
- }
- ## END get vcov
+ return(yuima)
+ })
+
+
+setMethod("confint", "yuima",
+ function(object, parm, level=0.95, ...){
+ yuima <- object
+ #print("This is yuima confint !")
+ Sigma1 <- yuima at data@mle at details$Sigma1
+ Sigma2 <- yuima at data@mle at details$Sigma2
+ #theta2 <- yuima at data@mle at coef[yuima at model@parameter at drift]
+ #theta1 <- yuima at data@mle at coef[yuima at model@parameter at diffusion]
+ coef <- yuima at data@mle at coef
- ## set min
- min <- opt$value
+ Calpha2 <- qnorm(level)
+
+ ## make matrix
+ ret <- matrix(0, length(yuima at model@parameter at all), 2)
+ rownames(ret) <- yuima at model@parameter at all
+ a <- (1 - level)/2
+ a <- c(a, 1 - a)
+ colnames(ret) <- paste(round(100 * a, 1), "%")
- ## set param name
- fullcoef <- coef
- coef.name <- NULL
- theta2.len <- length(yuima at model@parameter at drift)
- if( theta2.len != 1){
- for( i in 1:theta2.len){
- coef.name <- c(coef.name, paste("theta2.", i, sep=""))
- }
- }else{
- coef.name <- "theta2"
- }
- theta1.len <- length(yuima at model@parameter at diffusion)
- if( theta1.len != 1){
- for( i in 1:theta1.len){
- coef.name <- c(coef.name, paste("theta1.", i, sep=""))
- }
- }else{
- coef.name <- c(coef.name, "theta1")
+ ## get confint
+ Sigma.diag <- c( diag(Sigma2), diag(Sigma1) )
+ names(Sigma.diag) <- yuima at model@parameter at all
+ for(param in yuima at model@parameter at all){
+ ret[param,] <- c( (coef[param]-sqrt(Sigma.diag[param])*Calpha2),
+ (coef[param]+sqrt(Sigma.diag[param])*Calpha2))
}
- names(fullcoef) <- coef.name
- ## END set param name
+
+ return(ret)
- ## set method name
- method <- if(BFGS) "BFGS" else "Nelder-Mead"
-
- ## make mle object
- opt <- new("mle", call=call, coef=coef, fullcoef=unlist(fullcoef),
- vcov=vcov, min=min, details=opt, minuslogl=ql.opt,
- method=method)
-
- ## END opt convert
-
- return(opt)
})
-
Modified: pkg/yuima/R/yuima.data.R
===================================================================
--- pkg/yuima/R/yuima.data.R 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/R/yuima.data.R 2010-03-20 02:17:52 UTC (rev 67)
@@ -41,8 +41,18 @@
setMethod("get.zoo.data", signature(x="yuima.data"),
function(x){
return(x at zoo.data)
- })
+ })
+setGeneric("get.mle",
+ function(x)
+ standardGeneric("get.mle")
+ )
+
+setMethod("get.mle", signature(x="yuima.data"),
+ function(x){
+ return(x at mle)
+ })
+
# following funcs are basic generic funcs
setGeneric("plot",
@@ -111,6 +121,11 @@
# same methods for 'yuima'. Depend on same methods for 'data'
+setMethod("get.mle", "yuima",
+ function(x){
+ return(get.mle(x at data))
+ })
+
setMethod("get.zoo.data", "yuima",
function(x){
return(get.zoo.data(x at data))
Modified: pkg/yuima/man/adaBayes.Rd
===================================================================
--- pkg/yuima/man/adaBayes.Rd 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/man/adaBayes.Rd 2010-03-20 02:17:52 UTC (rev 67)
@@ -56,18 +56,18 @@
QL <- 0
param <- numeric(2)
for( i in 1:5){
- PAR <- ml.ql(yuima.mod,theta2=rand[i,1],theta1=rand[i,2],h=(n^(-2/3)))
- PAR <- PAR at details
- if(PAR$value > QL){
- QL <- PAR$value
- param <- PAR$par
+ PAR <- ml.ql(yuima.mod,theta2=rand[i,1],theta1=rand[i,2],h=(n^(-2/3)))
+ PAR <- get.mle(PAR)
+ if(PAR at min > QL){
+ QL <- PAR at min
+ param <- PAR at coef
}
}
BE <- adaBayes(yuima.mod,h=(n-1)^(-2/3),prior2=prior2,prior1=prior1,domain2=domain,domain1=domain,theta2.offset=param[1],theta1.offset=param[2])
print("True param")
print("theta2 = .6, theta1 = .3")
print("ML estimate")
- PAR$par
+ PAR at coef
print("Bayes estimate")
BE
}
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-03-20 02:17:52 UTC (rev 67)
@@ -11,7 +11,7 @@
\description{Calculate the quasi-likelihood and estimate of the parameters of the
stochastic differential equation by the maximum likelihood method.}
\usage{
-ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,BFGS=FALSE,param,interval)
+ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,method,param,interval)
ql(yuima,theta2,theta1,h,print=FALSE,param)
rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
}
@@ -23,7 +23,7 @@
parameters. Vector can be available only if theta is a scalar.}
\item{ptheta2,ptheta1}{}
\item{print}{you can see a progress of the estimation when print is TRUE.}
- \item{BFGS}{}
+ \item{method}{}
\item{param}{}
\item{interval}{}
\item{prevparam}{}
@@ -60,30 +60,30 @@
##QL
system.time(
-opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
+yuima <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
)
print("True param")
print("theta2 = .3, theta1 = .1")
print("ML estimator")
-print(opt)
+get.mle(yuima)@coef
## another way of parameter specification
##interval <- list(theta2.lim=c(0,1), theta1.lim=c(0,1))
##system.time(
-##opt <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
+##yuima <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
##)
##print("True param")
##print("theta2 = .3, theta1 = .1")
##print("ML estimator")
-##print(opt)
+#get.mle(yuima)@coef
system.time(
-opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), BFGS=TRUE)
+yuima <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newtoon")
)
print("True param")
print("theta2 = .3, theta1 = .1")
print("ML estimator")
-print(opt)
+get.mle(yuima)@coef
##multidimension case
##dXt^e = - drift.matrix * Xt^e * dt + diff.matrix * dWt
@@ -120,24 +120,57 @@
theta1.lim <- t( matrix( c(theta1.1.lim, theta1.2.lim), 2, 2) )
system.time(
-opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim)
+yuima <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim)
)
-print(opt)
+get.mle(yuima)@coef
## another way of parameter specification
#interval <- list(theta2.1.lim, theta2.2.lim, theta1.1.lim, theta1.2.lim)
#system.time(
-#opt <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
+#yuima <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
#)
-#print(opt)
+#get.mle(yuima)@coef
-system.time(
-opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim, BFGS=TRUE)
-)
-print(opt)
+
+##dXt^e = -theta2 * Xt^e * dt + theta1 * dWt
+diff.matrix <- matrix(c("theta1"), 1, 1)
+ymodel <- setModel(drift=c("(-1)*theta2*x"), diffusion=diff.matrix,
+ time.variable="t", state.variable="x", solve.variable="x")
+n <- 100
+
+ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
+yuima <- setYuima(model=ymodel, sampling=ysamp)
+set.seed(123)
+yuima <- simulate(yuima, xinit=1, true.parameter=c(0.5, 0.8))
+#0.3,0.1
+h <- 1/((n)^(2/3))
+
+onezoo <- function(ydata) {
+ dat <- get.zoo.data(ydata)
+ dats <- dat[[1]]
+ if(length(dat)>1) {
+ for(i in 2:(length(dat))) {
+ dats <- merge(dats,dat[[i]])
+ }
+ }
+ return(dats)
}
+theta2 <- 0.1#0.65
+theta1 <- 0.2#0.7
+print("use Newton method")
+res1 <- ml.ql(yuima, theta2, theta1, h, method="Newton", print=TRUE)
+print(confint(res1))
+
+
+print("use optim")
+res2 <- ml.ql(yuima, theta2, theta1, h, print=FALSE)
+print(confint(res2))
+
+
+}
+
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ts}
Modified: pkg/yuima/man/setData.Rd
===================================================================
--- pkg/yuima/man/setData.Rd 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/man/setData.Rd 2010-03-20 02:17:52 UTC (rev 67)
@@ -3,6 +3,7 @@
\alias{get.zoo.data}
\alias{dim}
\alias{length}
+\alias{get.mle}
\title{
Set and access data of an object of type "yuima.data" or "yuima".
@@ -14,6 +15,8 @@
\code{\link{yuima.data-class}} object. (Note: value is a \code{list} of
\code{\link{zoo}} objects).
+\code{get.mle} return the content of mle class of ml.ql result.
+
\code{plot} plot method for object of \code{\link{yuima.data-class}} or
\code{\link{yuima-class}}.
@@ -52,6 +55,11 @@
returns the content of the slot \code{zoo.data} of \code{x} if \code{x}
is of \code{\link{yuima.data-class}} or the content of
\code{x at data@zoo.data} if \code{x} is of \code{\link{yuima-class}}.
+
+The function \code{get.mle}
+returns the content of the slot \code{mle} of \code{x} if \code{x}
+is of \code{\link{yuima.data-class}} or the content of
+\code{x at data@mle} if \code{x} is of \code{\link{yuima-class}}.
}
\value{
\item{value}{a list of object(s) of \code{\link{yuima.data-class}} for
Modified: pkg/yuima/man/yuima-class.Rd
===================================================================
--- pkg/yuima/man/yuima-class.Rd 2010-03-09 04:57:55 UTC (rev 66)
+++ pkg/yuima/man/yuima-class.Rd 2010-03-20 02:17:52 UTC (rev 67)
@@ -2,6 +2,7 @@
\docType{class}
\alias{yuima-class}
+\alias{get.mle,yuima-method}
\alias{get.zoo.data,yuima-method}
\alias{plot,yuima-method}
\alias{dim,yuima-method}
@@ -68,7 +69,7 @@
\item{plot}{\code{signature(x = "yuima", \dots)}: calls
\code{\link{plot}} from the \code{\link{zoo}} package with argument
\code{x at data@zoo.data}. Additional arguments \code{\dots} are passed
- as is to the \code{\link{plot}} function.}
+ as is to the \code{\link{plot}} function.}
\item{dim}{\code{signature(x = "yuima")}: the number of SDEs in the
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 67
More information about the Yuima-commits
mailing list