[Yuima-commits] r534 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Dec 6 03:55:12 CET 2016
Author: kamatani
Date: 2016-12-06 03:55:11 +0100 (Tue, 06 Dec 2016)
New Revision: 534
Modified:
pkg/yuima/R/adaBayes.R
Log:
update adaBayes mcmc components
Modified: pkg/yuima/R/adaBayes.R
===================================================================
--- pkg/yuima/R/adaBayes.R 2016-12-02 14:44:51 UTC (rev 533)
+++ pkg/yuima/R/adaBayes.R 2016-12-06 02:55:11 UTC (rev 534)
@@ -1,460 +1,504 @@
-##::quasi-bayes function
-
-
-setGeneric("adaBayes",
- function(yuima, start,prior,lower,upper, method="nomcmc")
-standardGeneric("adaBayes")
-)
-setMethod("adaBayes", "yuima",
- function(yuima, start,prior,lower,upper, method="nomcmc")
-{
-
- joint <- FALSE
- fixed <- numeric(0)
- print <- FALSE
-
- call <- match.call()
-
- if( missing(yuima))
- yuima.stop("yuima object is missing.")
-
-## param handling
-
-## FIXME: maybe we should choose initial values at random within lower/upper
-## at present, qmle stops
- if( missing(start) )
- yuima.stop("Starting values for the parameters are missing.")
-
- diff.par <- yuima at model@parameter at diffusion
- drift.par <- yuima at model@parameter at drift
- jump.par <- yuima at model@parameter at jump
- measure.par <- yuima at model@parameter at measure
- common.par <- yuima at model@parameter at common
-
-## BEGIN Prior construction
- if(!missing(prior)){
- lower = numeric(0)
- upper = numeric(0)
- pdlist <- numeric(length(yuima at model@parameter at all))
- names(pdlist) <- yuima at model@parameter at all
- for(i in 1: length(pdlist)){
- if(prior[[names(pdlist)[i]]]$measure.type=="code"){
- expr <- prior[[names(pdlist)[i]]]$df
- code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", expr, perl=TRUE))
- args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", expr, perl=TRUE)), ","))
- pdlist[i] <- switch(code,
- dunif=paste("function(z){return(dunif(z, ", args[2], ", ", args[3],"))}"),
- dnorm=paste("function(z){return(dnorm(z,", args[2], ", ", args[3], "))}"),
- dbeta=paste("function(z){return(dbeta(z, ", args[2], ", ", args[3], "))}"),
- dgamma=paste("function(z){return(dgamma(z, ", args[2], ", ", args[3], "))}")
- )
- qf <- switch(code,
- dunif=paste("function(z){return(qunif(z, ", args[2], ", ", args[3],"))}"),
- dnorm=paste("function(z){return(qnorm(z,", args[2], ", ", args[3], "))}"),
- dbeta=paste("function(z){return(qbeta(z, ", args[2], ", ", args[3], "))}"),
- dgamma=paste("function(z){return(qgamma(z, ", args[2], ", ", args[3], "))}")
- )
- lower = append(lower,eval(parse("text"=qf))(0.05))
- upper = append(upper,eval(parse("text"=qf))(0.95))
- }
-
- }
- names(lower) <- names(pdlist)
- names(upper) <- names(pdlist)
-
- pd <- function(param){
- value <- 1
- for(i in 1:length(pdlist)){
- value <- value*eval(parse(text=pdlist[[i]]))(param[[i]])
- }
- return(value)
- }
- }else{
- pd <- function(param) return(1)
- }
-## END Prior construction
- JointOptim <- joint
- if(length(common.par)>0){
- JointOptim <- TRUE
- yuima.warn("Drift and diffusion parameters must be different. Doing
- joint estimation, asymptotic theory may not hold true.")
- }
-
-
- if(length(jump.par)+length(measure.par)>0)
- yuima.stop("Cannot estimate the jump models, yet")
-
- if(!is.list(start))
- yuima.stop("Argument 'start' must be of list type.")
-
- fullcoef <- NULL
-
- if(length(diff.par)>0)
- fullcoef <- diff.par
-
- if(length(drift.par)>0)
- fullcoef <- c(fullcoef, drift.par)
-
- npar <- length(fullcoef)
-
- fixed.par <- names(fixed)
-
- if (any(!(fixed.par %in% fullcoef)))
- yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
-
- nm <- names(start)
- oo <- match(nm, fullcoef)
- if(any(is.na(oo)))
- yuima.stop("some named arguments in 'start' are not arguments to the supplied yuima model")
- start <- start[order(oo)]
- nm <- names(start)
-
- idx.diff <- match(diff.par, nm)
- idx.drift <- match(drift.par, nm)
- idx.fixed <- match(fixed.par, nm)
-
- tmplower <- as.list( rep( -Inf, length(nm)))
- names(tmplower) <- nm
- if(!missing(lower)){
- idx <- match(names(lower), names(tmplower))
- if(any(is.na(idx)))
- yuima.stop("names in 'lower' do not match names fo parameters")
- tmplower[ idx ] <- lower
- }
- lower <- tmplower
-
- tmpupper <- as.list( rep( Inf, length(nm)))
- names(tmpupper) <- nm
- if(!missing(upper)){
- idx <- match(names(upper), names(tmpupper))
- if(any(is.na(idx)))
- yuima.stop("names in 'lower' do not match names fo parameters")
- tmpupper[ idx ] <- upper
- }
- upper <- tmpupper
-
-
-
-
- d.size <- yuima at model@equation.number
- n <- length(yuima)[1]
-
- env <- new.env()
- assign("X", as.matrix(onezoo(yuima)), envir=env)
- assign("deltaX", matrix(0, n-1, d.size), envir=env)
- assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=env)
-
- assign("Cn.r", rep(1,n-1), envir=env)
-
- for(t in 1:(n-1))
- env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
-
- assign("h", deltat(yuima at data@zoo.data[[1]]), envir=env)
-
- mle <- qmle(yuima, "start"=start, "lower"=lower,"upper"=upper, "method"="L-BFGS-B")
- integ <- function(idx.fixed=NULL,f=f,start=start,par=NULL,hessian=FALSE,upper,lower){
- intf <- adaptIntegrate(f,lowerLimit=lower[-idx.fixed],upperLimit=upper[-idx.fixed],fDim=(length(upper[-idx.fixed])+1))$integral
- return(intf[-1]/intf[1])
- }
- mcinteg <- function(idx.fixed=NULL,f=f,p,start=start,par=NULL,hessian=FALSE,upper,lower,mean,vcov){
- intf <- mcIntegrate(f,p,lowerLimit=lower[-idx.fixed],upperLimit=upper[-idx.fixed],mean[-idx.fixed],vcov[-idx.fixed,-idx.fixed],mcmc=1000)
- return(intf)
- }
-
- mcIntegrate <- function(f,p, lowerLimit, upperLimit,mean,vcov,mcmc=1000){
- x_c <- mean
- p_c <- p(mean)
- val <- f(x_c)
-
- if(length(mean)>1){
- x <- rmvnorm(mcmc-1,mean,vcov)
- q <- dmvnorm(x,mean,vcov)
- q_c <- dmvnorm(mean,mean,vcov)
- }else{
- x <- rnorm(mcmc-1,mean,sqrt(vcov))
- q <- dnorm(x,mean,sqrt(vcov))
- q_c <- dnorm(mean,mean,sqrt(vcov))
- }
-
- for(i in 1:(mcmc-1)){
- if(length(mean)>1){x_n <- x[i,]}else{x_n <- x[i]}
- if(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0){
- q_n <- q[i]
- p_n <- p(x_n)
- u <- runif(1)
- a <- (p_n*q_c)/(p_c*q_n)
-
- if(u<a){
- p_c <- p_n
- q_c <- q_n
- x_c <- x_n
- }
- }
- val <- val+f(x_c)
- }
- return(unlist(val/mcmc))
- }
-
- print(mle at coef)
- tmp <- minusquasilogl(yuima=yuima, param=mle at coef, print=print, env)
-
- g <- function(p,fixed,idx.fixed){
- mycoef <- mle at coef
- if(length(idx.fixed)>0){
- mycoef[-idx.fixed] <- p
- mycoef[idx.fixed] <- fixed
- }else{
- names(mycoef) <- nm
- }
-
- return(c(1,p)*exp(-minusquasilogl(yuima=yuima, param=mycoef, print=print, env)+tmp)*pd(param=mycoef))
- }
-
- pg <- function(p,fixed,idx.fixed){
- mycoef <- start
- if(length(idx.fixed)>0){
- mycoef[-idx.fixed] <- p
- mycoef[idx.fixed] <- fixed
- }else{
- names(mycoef) <- nm
- }
-
- return(exp(-minusquasilogl(yuima=yuima, param=mycoef, print=print, env)+tmp)*pd(param=mycoef))
- }
-
- idf <- function(p){return(p)}
-
-# fj <- function(p) {
-# mycoef <- as.list(p)
-# names(mycoef) <- nm
-# mycoef[fixed.par] <- fixed
-# minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
-# }
-
- oout <- NULL
- HESS <- matrix(0, length(nm), length(nm))
- colnames(HESS) <- nm
- rownames(HESS) <- nm
- HaveDriftHess <- FALSE
- HaveDiffHess <- FALSE
- if(length(start)){
-# if(JointOptim){ ### joint optimization
-# if(length(start)>1){ #multidimensional optim
-# oout <- optim(start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
-# HESS <- oout$hessian
-# HaveDriftHess <- TRUE
-# HaveDiffHess <- TRUE
-# } else { ### one dimensional optim
-# opt1 <- optimize(f, ...) ## an interval should be provided
-# opt1 <- list(par=integ(f=f,upper=upper,lower=lower,fDim=length(lower)+1),objective=0)
-# oout <- list(par = opt1$minimum, value = opt1$objective)
-# } ### endif( length(start)>1 )
-# } else { ### first diffusion, then drift
- theta1 <- NULL
-
- old.fixed <- fixed
- old.start <- start
-
- if(length(idx.diff)>0){
-## DIFFUSION ESTIMATIOn first
-
- old.fixed <- fixed
- old.start <- start
- new.start <- start[idx.diff] # considering only initial guess for diffusion
- new.fixed <- fixed
- if(length(idx.drift)>0)
- new.fixed[nm[idx.drift]] <- start[idx.drift]
- fixed <- new.fixed
- fixed.par <- names(fixed)
- idx.fixed <- match(fixed.par, nm)
- names(new.start) <- nm[idx.diff]
-
- mydots <- as.list(call)[-(1:2)]
- mydots$fixed <- NULL
- mydots$fn <- as.name("f")
- mydots$start <- NULL
- mydots$par <- unlist(new.start)
- mydots$hessian <- FALSE
- mydots$upper <- unlist( upper[ nm[idx.diff] ])
- mydots$lower <- unlist( lower[ nm[idx.diff] ])
- f <- function(p){return(g(p,fixed,idx.fixed))}
- pf <- function(p){return(pg(p,fixed,idx.fixed))}
- if(length(mydots$par)>1){
-# oout <- do.call(optim, args=mydots)
- if(method=="mcmc"){
- oout <- list(par=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=diag(diag(mle at vcov))))
- }else{
- oout <- list(par=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower,start=start))
- }
- } else {
- mydots$f <- mydots$fn
- mydots$fn <- NULL
- mydots$par <- NULL
- mydots$hessian <- NULL
- mydots$method <- NULL
- mydots$interval <- as.numeric(c(unlist(lower[diff.par]),unlist(upper[diff.par])))
- mydots$lower <- NULL
- mydots$upper <- NULL
-# opt1 <- do.call(optimize, args=mydots)
- if(method=="mcmc"){
- opt1 <- list(minimum=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=diag(diag(mle at vcov))))
- }else{
- opt1 <- list(minimum=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower))
- }
-
- theta1 <- opt1$minimum
- names(theta1) <- diff.par
-# oout <- list(par = theta1, value = opt1$objective)
- oout <- list(par=theta1,value=0)
- }
- theta1 <- oout$par
- } ## endif(length(idx.diff)>0)
-
- theta2 <- NULL
-
- if(length(idx.drift)>0){
-## DRIFT estimation with first state diffusion estimates
- fixed <- old.fixed
- start <- old.start
- new.start <- start[idx.drift] # considering only initial guess for drift
- new.fixed <- fixed
- new.fixed[names(theta1)] <- theta1
- fixed <- new.fixed
- fixed.par <- names(fixed)
- idx.fixed <- match(fixed.par, nm)
- names(new.start) <- nm[idx.drift]
-
- mydots <- as.list(call)[-(1:2)]
- mydots$fixed <- NULL
- mydots$fn <- as.name("f")
- mydots$start <- NULL
- mydots$par <- unlist(new.start)
- mydots$hessian <- FALSE
- mydots$upper <- unlist( upper[ nm[idx.drift] ])
- mydots$lower <- unlist( lower[ nm[idx.drift] ])
- f <- function(p){return(g(p,fixed,idx.fixed))}
- pf <- function(p){return(pg(p,fixed,idx.fixed))}
-
- if(length(mydots$par)>1){
-# oout1 <- do.call(optim, args=mydots)
- if(method=="mcmc"){
- oout1 <- list(par=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=diag(diag(mle at vcov))))
- }else{
- oout1 <- list(par=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower))
- }
- } else {
- mydots$f <- mydots$fn
- mydots$fn <- NULL
- mydots$par <- NULL
- mydots$hessian <- NULL
- mydots$method <- NULL
- mydots$interval <- as.numeric(c(lower[drift.par],upper[drift.par]))
-# opt1 <- do.call(optimize, args=mydots)
- if(method=="mcmc"){
- opt1 <- list(minimum=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=diag(diag(mle at vcov))))
- }else{
- opt1 <- list(minimum=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower))
- }
- theta2 <- opt1$minimum
- names(theta2) <- drift.par
- oout1 <- list(par = theta2, value = as.numeric(opt1$objective))
- }
- theta2 <- oout1$par
- } ## endif(length(idx.drift)>0)
- oout1 <- list(par= c(theta1, theta2))
- names(oout1$par) <- c(diff.par,drift.par)
- oout <- oout1
-
-# } ### endif JointOptim
- } else {
- list(par = numeric(0L), value = f(start))
- }
-
-
- fDrift <- function(p) {
- mycoef <- as.list(p)
- names(mycoef) <- drift.par
- mycoef[diff.par] <- coef[diff.par]
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
- }
-
- fDiff <- function(p) {
- mycoef <- as.list(p)
- names(mycoef) <- diff.par
- mycoef[drift.par] <- coef[drift.par]
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
- }
-
- coef <- oout$par
- control=list()
- par <- coef
- names(par) <- c(diff.par, drift.par)
- nm <- c(diff.par, drift.par)
-
-# print(par)
-# print(coef)
- conDrift <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(drift.par)),
- ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
- factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
- conDiff <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(diff.par)),
- ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
- factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
-
-# nmsC <- names(con)
-# if (method == "Nelder-Mead")
-# con$maxit <- 500
-# if (method == "SANN") {
-# con$maxit <- 10000
-# con$REPORT <- 100
-# }
-# con[(namc <- names(control))] <- control
-# if (length(noNms <- namc[!namc %in% nmsC]))
-# warning("unknown names in control: ", paste(noNms, collapse = ", "))
-# if (con$trace < 0)
-# warning("read the documentation for 'trace' more carefully")
-# else if (method == "SANN" && con$trace && as.integer(con$REPORT) ==
-# 0)
-# stop("'trace != 0' needs 'REPORT >= 1'")
-# if (method == "L-BFGS-B" && any(!is.na(match(c("reltol",
-# "abstol"), namc))))
-# warning("method L-BFGS-B uses 'factr' (and 'pgtol') instead of 'reltol' and 'abstol'")
-# npar <- length(par)
-# if (npar == 1 && method == "Nelder-Mead")
-# warning("one-diml optimization by Nelder-Mead is unreliable: use optimize")
-#
- if(!HaveDriftHess & (length(drift.par)>0)){
- #hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
- hess2 <- optimHess(coef[drift.par], fDrift, NULL, control=conDrift)
- HESS[drift.par,drift.par] <- hess2
- }
-
- if(!HaveDiffHess & (length(diff.par)>0)){
- #hess1 <- .Internal(optimhess(coef[diff.par], fDiff, NULL, conDiff))
- hess1 <- optimHess(coef[diff.par], fDiff, NULL, control=conDiff)
- HESS[diff.par,diff.par] <- hess1
- }
-
- oout$hessian <- HESS
-
- vcov <- if (length(coef))
- solve(oout$hessian)
- else matrix(numeric(0L), 0L, 0L)
-
- mycoef <- as.list(coef)
- names(mycoef) <- nm
- mycoef[fixed.par] <- fixed
-
- min <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
-
- new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
-# vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
- vcov = vcov, details = oout,
- method = method)
-}
-)
+##::quasi-bayes function
+
+
+setGeneric("adaBayes",
+function(yuima, start,prior,lower,upper, method="mcmc",mcmc=1000,rate=1.0,rcpp=TRUE,algorithm="randomwalk")
+standardGeneric("adaBayes")
+)
+setMethod("adaBayes", "yuima",
+function(yuima, start,prior,lower,upper, method="mcmc",mcmc=1000,rate=1.0,rcpp=TRUE,algorithm="randomwalk")
+{
+
+
+
+
+
+ joint <- FALSE
+ fixed <- numeric(0)
+ print <- FALSE
+
+ call <- match.call()
+
+ if( missing(yuima))
+ yuima.stop("yuima object is missing.")
+
+ ## param handling
+
+ ## FIXME: maybe we should choose initial values at random within lower/upper
+ ## at present, qmle stops
+
+ if(missing(lower) || missing(upper)){
+ yuima.stop("lower or upper is missing.")
+ }
+
+ diff.par <- yuima at model@parameter at diffusion
+ drift.par <- yuima at model@parameter at drift
+ jump.par <- yuima at model@parameter at jump
+ measure.par <- yuima at model@parameter at measure
+ common.par <- yuima at model@parameter at common
+
+ ## BEGIN Prior construction
+ if(!missing(prior)){
+ priorLower = numeric(0)
+ priorUpper = numeric(0)
+ pdlist <- numeric(length(yuima at model@parameter at all))
+ names(pdlist) <- yuima at model@parameter at all
+ for(i in 1: length(pdlist)){
+ if(prior[[names(pdlist)[i]]]$measure.type=="code"){
+ expr <- prior[[names(pdlist)[i]]]$df
+ code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", expr, perl=TRUE))
+ args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", expr, perl=TRUE)), ","))
+ pdlist[i] <- switch(code,
+ dunif=paste("function(z){return(dunif(z, ", args[2], ", ", args[3],"))}"),
+ dnorm=paste("function(z){return(dnorm(z,", args[2], ", ", args[3], "))}"),
+ dbeta=paste("function(z){return(dbeta(z, ", args[2], ", ", args[3], "))}"),
+ dgamma=paste("function(z){return(dgamma(z, ", args[2], ", ", args[3], "))}"),
+ dexp=paste("function(z){return(dexp(z, ", args[2], "))}")
+ )
+ qf <- switch(code,
+ dunif=paste("function(z){return(qunif(z, ", args[2], ", ", args[3],"))}"),
+ dnorm=paste("function(z){return(qnorm(z,", args[2], ", ", args[3], "))}"),
+ dbeta=paste("function(z){return(qbeta(z, ", args[2], ", ", args[3], "))}"),
+ dgamma=paste("function(z){return(qgamma(z, ", args[2], ", ", args[3], "))}"),
+ dexp=paste("function(z){return(qexp(z, ", args[2], "))}")
+ )
+ priorLower = append(priorLower,eval(parse("text"=qf))(0.00))
+ priorUpper = append(priorUpper,eval(parse("text"=qf))(1.00))
+
+
+ }
+
+ }
+ if(sum(unlist(lower)<priorLower) + sum(unlist(upper)>priorUpper) > 0){
+ yuima.stop("lower&upper of prior are out of parameter space.")
+ }
+
+ names(lower) <- names(pdlist)
+ names(upper) <- names(pdlist)
+
+
+
+ pd <- function(param){
+ value <- 1
+ for(i in 1:length(pdlist)){
+ value <- value*eval(parse(text=pdlist[[i]]))(param[[i]])
+ }
+ return(value)
+ }
+ }else{
+ pd <- function(param) return(1)
+ }
+ ## END Prior construction
+
+ if(!is.list(start) || (sum(unlist(start)<unlist(lower))+sum(unlist(start)>unlist(upper))>0)){
+ #cannot use "missing(start)"
+ start <- lower
+ start[1:length(start)] <- runif(length(start),unlist(lower),unlist(upper))
+ #yuima.warn("param.init is out of parameter space.redefigned init by runif.")
+ }
+
+ JointOptim <- joint
+ if(length(common.par)>0){
+ JointOptim <- TRUE
+ yuima.warn("Drift and diffusion parameters must be different. Doing
+ joint estimation, asymptotic theory may not hold true.")
+ }
+
+
+ if(length(jump.par)+length(measure.par)>0)
+ yuima.stop("Cannot estimate the jump models, yet")
+
+
+ fullcoef <- NULL
+
+ if(length(diff.par)>0)
+ fullcoef <- diff.par
+
+ if(length(drift.par)>0)
+ fullcoef <- c(fullcoef, drift.par)
+
+ npar <- length(fullcoef)
+
+ fixed.par <- names(fixed)
+
+ if (any(!(fixed.par %in% fullcoef)))
+ yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
+
+ nm <- names(start)
+ oo <- match(nm, fullcoef)
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'start' are not arguments to the supplied yuima model")
+ start <- start[order(oo)]
+ if(!missing(prior)){
+ pdlist <- pdlist[order(oo)]
+ }
+ nm <- names(start)
+
+ idx.diff <- match(diff.par, nm)
+ idx.drift <- match(drift.par, nm)
+ idx.fixed <- match(fixed.par, nm)
+ tmplower <- as.list( rep( -Inf, length(nm)))
+ names(tmplower) <- nm
+ if(!missing(lower)){
+ idx <- match(names(lower), names(tmplower))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmplower[ idx ] <- lower
+ }
+ lower <- tmplower
+
+ tmpupper <- as.list( rep( Inf, length(nm)))
+ names(tmpupper) <- nm
+ if(!missing(upper)){
+ idx <- match(names(upper), names(tmpupper))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmpupper[ idx ] <- upper
+ }
+ upper <- tmpupper
+
+
+
+
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ G <- rate
+ if(G<=0 || G>1){
+ yuima.stop("rate G should be 0 < G <= 1")
+ }
+ n_0 <- floor(n^G)
+ if(n_0 < 2) n_0 <- 2
+
+ #######data is reduced to n_0 before qmle(16/11/2016)
+ env <- new.env()
+ #assign("X", yuima at data@original.data[1:n_0,], envir=env)
+ assign("X", as.matrix(onezoo(yuima))[1:n_0,], envir=env)
+ assign("deltaX", matrix(0, n_0 - 1, d.size), envir=env)
+ assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=env)
+
+ assign("Cn.r", rep(1,n_0-1), envir=env)
+
+ for(t in 1:(n_0-1))
+ env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+
+ assign("h", deltat(yuima at data@zoo.data[[1]]), envir=env)
+
+ mle <- qmle(yuima, "start"=start, "lower"=lower,"upper"=upper, "method"="L-BFGS-B",rcpp=rcpp)
+ integ <- function(idx.fixed=NULL,f=f,start=start,par=NULL,hessian=FALSE,upper,lower){
+ if(length(idx.fixed)==0){
+ intf <- adaptIntegrate(f,lowerLimit=lower,upperLimit=upper,fDim=(length(upper)+1))$integral
+ }else{
+ intf <- adaptIntegrate(f,lowerLimit=lower[-idx.fixed],upperLimit=upper[-idx.fixed],fDim=(length(upper[-idx.fixed])+1))$integral
+ }
+ return(intf[-1]/intf[1])
+ }
+ mcinteg <- function(idx.fixed=NULL,f=f,p,start=start,par=NULL,hessian=FALSE,upper,lower,mean,vcov,mcmc){
+ if(length(idx.fixed)==0){
+ intf <- mcIntegrate(f,p,lowerLimit=lower,upperLimit=upper,mean,vcov,mcmc)
+ }else{
+ intf <- mcIntegrate(f,p,lowerLimit=lower[-idx.fixed],upperLimit=upper[-idx.fixed],mean[-idx.fixed],vcov[-idx.fixed,-idx.fixed],mcmc)
+ }
+ return(intf)
+ }
+
+ mcIntegrate <- function(f,p, lowerLimit, upperLimit,mean,vcov,mcmc){
+
+ if(algorithm=="randomwalk"){
+ x_c <- mean
+ p_c <- p(mean)
+ val <- f(x_c)
+
+ if(length(mean)>1){
+ x <- rmvnorm(mcmc-1,mean,vcov)
+ q <- dmvnorm(x,mean,vcov)
+ q_c <- dmvnorm(mean,mean,vcov)
+ }else{
+ x <- rnorm(mcmc-1,mean,sqrt(vcov))
+ q <- dnorm(x,mean,sqrt(vcov))
+ q_c <- dnorm(mean,mean,sqrt(vcov))
+ }
+
+ for(i in 1:(mcmc-1)){
+ if(length(mean)>1){x_n <- x[i,]}else{x_n <- x[i]}
+ if(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0){
+ q_n <- q[i]
+ p_n <- p(x_n)
+ #u <- runif(1)
+ #a <- (p_n*q_c)/(p_c*q_n)
+ u <- log(runif(1))
+ a <- p_n-p_c+log(q_c/q_n)
+ if(u<a){
+ p_c <- p_n
+ q_c <- q_n
+ x_c <- x_n
+ }
+ }
+ val <- val+f(x_c)
+ }
+ return(unlist(val/mcmc))
+ }
+ else if(algorithm=="MpCN"){
+ x_n <- mean
+ val <- mean
+ logLik_old <- p(mean)+0.5*length(mean)*log(sqnorm(x_n-mean))
+
+ for(i in 1:(mcmc-1)){
+ #browser()
+ prop <- makeprop(mean,x_n,unlist(lowerLimit),unlist(upperLimit))
+ logLik_new <- p(mean)+0.5*length(mean)*log(sqnorm(prop-mean))
+ u <- log(runif(1))
+ if( logLik_new-logLik_old > u){
+ x_n <- prop
+ logLik_old <- logLik_new
+ }
+ val <- val+f(x_n)
+ }
+ return(unlist(val/mcmc))
+ }
+ }
+
+ #print(mle at coef)
+
+
+ flagNotPosDif <- 0
+ for(i in 1:npar){
+ if(mle at vcov[i,i] <= 0) flagNotPosDif <- 1 #Check mle at vcov is positive difinite matrix
+ }
+ if(flagNotPosDif == 1){
+ mle at vcov <- diag(c(rep(1 / n,length(diff.par)),rep(1 / (n * env$h),length(drift.par)))) # Redifine mle at vcov
+ }
+
+
+ tmp <- minusquasilogl(yuima=yuima, param=mle at coef, print=print, env,rcpp=rcpp)
+
+ g <- function(p,fixed,idx.fixed){
+ mycoef <- mle at coef
+ if(length(idx.fixed)>0){
+ mycoef[-idx.fixed] <- p
+ mycoef[idx.fixed] <- fixed
+ }else{
+ names(mycoef) <- nm
+ }
+ return(c(1,p)*exp(-minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)+tmp)*pd(param=mycoef))
+ }
+
+ pg <- function(p,fixed,idx.fixed){
+ mycoef <- start
+ if(length(idx.fixed)>0){
+ mycoef[-idx.fixed] <- p
+ mycoef[idx.fixed] <- fixed
+ }else{
+ names(mycoef) <- nm
+ }
+ #return(exp(-minusquasilogl(yuima=yuima, param=mycoef, print=print, env)+tmp)*pd(param=mycoef))
+ return(-minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)+tmp+log(pd(param=mycoef)))#log
+ }
+
+ idf <- function(p){return(p)}
+
+ # fj <- function(p) {
+ # mycoef <- as.list(p)
+ # names(mycoef) <- nm
+ # mycoef[fixed.par] <- fixed
+ # minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ # }
+
+ oout <- NULL
+ HESS <- matrix(0, length(nm), length(nm))
+ colnames(HESS) <- nm
+ rownames(HESS) <- nm
+ HaveDriftHess <- FALSE
+ HaveDiffHess <- FALSE
+ if(length(start)){
+ # if(JointOptim){ ### joint optimization
+ # if(length(start)>1){ #multidimensional optim
+ # oout <- optim(start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
+ # HESS <- oout$hessian
+ # HaveDriftHess <- TRUE
+ # HaveDiffHess <- TRUE
+ # } else { ### one dimensional optim
+ # opt1 <- optimize(f, ...) ## an interval should be provided
+ # opt1 <- list(par=integ(f=f,upper=upper,lower=lower,fDim=length(lower)+1),objective=0)
+ # oout <- list(par = opt1$minimum, value = opt1$objective)
+ # } ### endif( length(start)>1 )
+ # } else { ### first diffusion, then drift
+ theta1 <- NULL
+
+ old.fixed <- fixed
+ old.start <- start
+
+ if(length(idx.diff)>0){
+ ## DIFFUSION ESTIMATIOn first
+ old.fixed <- fixed
+ old.start <- start
+ new.start <- start[idx.diff] # considering only initial guess for diffusion
+ new.fixed <- fixed
+ if(length(idx.drift)>0)
+ new.fixed[nm[idx.drift]] <- start[idx.drift]
+ fixed <- new.fixed
+ fixed.par <- names(fixed)
+ idx.fixed <- match(fixed.par, nm)
+ names(new.start) <- nm[idx.diff]
+
+ f <- function(p){return(g(p,fixed,idx.fixed))}
+ pf <- function(p){return(pg(p,fixed,idx.fixed))}
+ if(length(unlist(new.start))>1){
+ # oout <- do.call(optim, args=mydots)
+ if(method=="mcmc"){
+ oout <- list(par=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=mle at vcov,mcmc=mcmc))
+ }else{
+ oout <- list(par=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower,start=start))
+ }
+ } else {
+ # opt1 <- do.call(optimize, args=mydots)
+ if(method=="mcmc"){
+ opt1 <- list(minimum=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=mle at vcov,mcmc=mcmc))
+ }else{
+ opt1 <- list(minimum=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower))
+ }
+ theta1 <- opt1$minimum
+ names(theta1) <- diff.par
+ # oout <- list(par = theta1, value = opt1$objective)
+ oout <- list(par=theta1,value=0)
+ }
+ theta1 <- oout$par
+ #names(theta1) <- nm[idx.diff]
+ names(theta1) <- diff.par
+ } ## endif(length(idx.diff)>0)
+
+ theta2 <- NULL
+
+ if(length(idx.drift)>0){
+ ## DRIFT estimation with first state diffusion estimates
+ fixed <- old.fixed
+ start <- old.start
+ new.start <- start[idx.drift] # considering only initial guess for drift
+ new.fixed <- fixed
+ new.fixed[names(theta1)] <- theta1
+ fixed <- new.fixed
+ fixed.par <- names(fixed)
+ idx.fixed <- match(fixed.par, nm)
+ names(new.start) <- nm[idx.drift]
+
+ f <- function(p){return(g(p,fixed,idx.fixed))}
+ pf <- function(p){return(pg(p,fixed,idx.fixed))}
+
+ if(length(unlist(new.start))>1){
+ # oout1 <- do.call(optim, args=mydots)
+ if(method=="mcmc"){
+ oout1 <- list(par=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=mle at vcov,mcmc=mcmc))
+ }else{
+ oout1 <- list(par=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower))
+ }
+ } else {
+ # opt1 <- do.call(optimize, args=mydots)
+ if(method=="mcmc"){
+ opt1 <- list(minimum=mcinteg(idx.fixed=idx.fixed,f=idf,p=pf,upper=upper,lower=lower,mean=mle at coef,vcov=mle at vcov,mcmc=mcmc))
+ }else{
+ opt1 <- list(minimum=integ(idx.fixed=idx.fixed,f=f,upper=upper,lower=lower))
+ }
+ theta2 <- opt1$minimum
+ names(theta2) <- drift.par
+ oout1 <- list(par = theta2, value = as.numeric(opt1$objective))
+ }
+ theta2 <- oout1$par
+ } ## endif(length(idx.drift)>0)
+ oout1 <- list(par= c(theta1, theta2))
+ names(oout1$par) <- c(diff.par,drift.par)
+ oout <- oout1
+
+ # } ### endif JointOptim
+ } else {
+ list(par = numeric(0L), value = f(start))
+ }
+
+
+ fDrift <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- drift.par
+ mycoef[diff.par] <- coef[diff.par]
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
+ }
+
+ fDiff <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- diff.par
+ mycoef[drift.par] <- coef[drift.par]
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
+ }
+
+ coef <- oout$par
+ control=list()
+ par <- coef
+ names(par) <- c(diff.par, drift.par)
+ nm <- c(diff.par, drift.par)
+
+ # print(par)
+ # print(coef)
+ conDrift <- list(trace = 5, fnscale = 1,
+ parscale = rep.int(5, length(drift.par)),
+ ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
+ conDiff <- list(trace = 5, fnscale = 1,
+ parscale = rep.int(5, length(diff.par)),
+ ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
+
+ # nmsC <- names(con)
+ # if (method == "Nelder-Mead")
+ # con$maxit <- 500
+ # if (method == "SANN") {
+ # con$maxit <- 10000
+ # con$REPORT <- 100
+ # }
+ # con[(namc <- names(control))] <- control
+ # if (length(noNms <- namc[!namc %in% nmsC]))
+ # warning("unknown names in control: ", paste(noNms, collapse = ", "))
+ # if (con$trace < 0)
+ # warning("read the documentation for 'trace' more carefully")
+ # else if (method == "SANN" && con$trace && as.integer(con$REPORT) ==
+ # 0)
+ # stop("'trace != 0' needs 'REPORT >= 1'")
+ # if (method == "L-BFGS-B" && any(!is.na(match(c("reltol",
+ # "abstol"), namc))))
+ # warning("method L-BFGS-B uses 'factr' (and 'pgtol') instead of 'reltol' and 'abstol'")
+ # npar <- length(par)
+ # if (npar == 1 && method == "Nelder-Mead")
+ # warning("one-diml optimization by Nelder-Mead is unreliable: use optimize")
+ #
+ if(!HaveDriftHess & (length(drift.par)>0)){
+ #hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
+ hess2 <- optimHess(coef[drift.par], fDrift, NULL, control=conDrift)
+ HESS[drift.par,drift.par] <- hess2
+ }
+
+ if(!HaveDiffHess & (length(diff.par)>0)){
+ #hess1 <- .Internal(optimhess(coef[diff.par], fDiff, NULL, conDiff))
+ hess1 <- optimHess(coef[diff.par], fDiff, NULL, control=conDiff)
+ HESS[diff.par,diff.par] <- hess1
+ }
+
+ oout$hessian <- HESS
+
+ vcov <- if (length(coef))
+ solve(oout$hessian)
+ else matrix(numeric(0L), 0L, 0L)
+
+ mycoef <- as.list(coef)
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+
+ #min <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
+
+ new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ # vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
+ vcov = vcov, details = oout,
+ method = method)
+ }
+)
+
More information about the Yuima-commits
mailing list