[Yuima-commits] r57 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 19 23:54:01 CET 2010
Author: hinohide
Date: 2010-01-19 23:54:00 +0100 (Tue, 19 Jan 2010)
New Revision: 57
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/R/yuima.data.R
pkg/yuima/man/adaBayes.Rd
pkg/yuima/man/quasi-likelihood.Rd
Log:
ml.ql output form modifications
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-01-08 03:27:38 UTC (rev 56)
+++ pkg/yuima/DESCRIPTION 2010-01-19 22:54:00 UTC (rev 57)
@@ -1,12 +1,12 @@
-Package: yuima
-Type: Package
-Title: The YUIMA Project package
-Version: 0.0.85
-Date: 2010-1-8
-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
- for Stochastic Differential Equations.
-License: GPL-2
-URL: http://R-Forge.R-project.org/projects/yuima/
+Package: yuima
+Type: Package
+Title: The YUIMA Project package
+Version: 0.0.82
+Date: 2010-01-01
+Depends: methods, zoo, adapt, stats4
+Author: YUIMA Project Team.
+Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
+Description: The YUIMA Project for Simulation and Inference
+ for Stochastic Differential Equations.
+License: GPL-2
+URL: http://R-Forge.R-project.org/projects/yuima/
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2010-01-08 03:27:38 UTC (rev 56)
+++ pkg/yuima/R/quasi-likelihood.R 2010-01-19 22:54:00 UTC (rev 57)
@@ -582,21 +582,164 @@
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))
+ 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))
+ return(ql.grad(yuima, theta2=theta[1:length(theta2)],
+ theta1=theta[(length(theta2)+1):length(theta)], h=h, print=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
+ }
+ 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 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
+ }
+ ## END constrOptim2 ( original code is constrOptim() in Package stat )
+
+ ## optimize
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)
+ 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 <- constrOptim(c(theta2, theta1), ql.opt, NULL, ui=ui, ci=ci, control=list(fnscale=-1), outer.iterations=500)
+ 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()
+ call <- match.call()
+
+ ## set par
+ coef <- opt$par
+
+ ## 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
+
+ ## set min
+ min <- opt$value
+
+ ## 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")
+ }
+ names(fullcoef) <- coef.name
+ ## END set param name
+
+ ## 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-01-08 03:27:38 UTC (rev 56)
+++ pkg/yuima/R/yuima.data.R 2010-01-19 22:54:00 UTC (rev 57)
@@ -128,8 +128,8 @@
setMethod("plot","yuima",
function(x,y,xlab=x at model@time.variable,ylab=x at model@solve.variable,...){
if(length(x at model@time.variable)==0) {
- plot(x at data,y,...)
+ plot(x at data,...)
} else {
- plot(x at data,y,xlab=xlab,ylab=ylab,...)
+ plot(x at data,xlab=xlab,ylab=ylab,...)
}
})
Modified: pkg/yuima/man/adaBayes.Rd
===================================================================
--- pkg/yuima/man/adaBayes.Rd 2010-01-08 03:27:38 UTC (rev 56)
+++ pkg/yuima/man/adaBayes.Rd 2010-01-19 22:54:00 UTC (rev 57)
@@ -56,7 +56,8 @@
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 <- 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
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-01-08 03:27:38 UTC (rev 56)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-01-19 22:54:00 UTC (rev 57)
@@ -65,7 +65,7 @@
print("True param")
print("theta2 = .3, theta1 = .1")
print("ML estimator")
-opt$par
+print(opt)
## another way of parameter specification
##interval <- list(theta2.lim=c(0,1), theta1.lim=c(0,1))
@@ -75,7 +75,7 @@
##print("True param")
##print("theta2 = .3, theta1 = .1")
##print("ML estimator")
-##opt$par
+##print(opt)
system.time(
opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), BFGS=TRUE)
@@ -83,7 +83,7 @@
print("True param")
print("theta2 = .3, theta1 = .1")
print("ML estimator")
-opt$par
+print(opt)
##multidimension case
##dXt^e = - drift.matrix * Xt^e * dt + diff.matrix * dWt
@@ -122,19 +122,19 @@
system.time(
opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim)
)
-opt$par
+print(opt)
## 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)
#)
-#opt$par
+#print(opt)
system.time(
opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim, BFGS=TRUE)
)
-opt$par
+print(opt)
}
More information about the Yuima-commits
mailing list