[Yuima-commits] r739 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 5 05:52:55 CET 2021
Author: eguchi
Date: 2021-02-05 05:52:54 +0100 (Fri, 05 Feb 2021)
New Revision: 739
Modified:
pkg/yuima/R/adaBayes.R
Log:
modified
Modified: pkg/yuima/R/adaBayes.R
===================================================================
--- pkg/yuima/R/adaBayes.R 2020-12-20 19:24:37 UTC (rev 738)
+++ pkg/yuima/R/adaBayes.R 2021-02-05 04:52:54 UTC (rev 739)
@@ -1,560 +1,859 @@
-##::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)){
- 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)){
- lower = append(lower,-Inf)
- upper = append(upper,Inf)
- }
- }
- # 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(missing(lower) || missing(upper)){
- lower <- priorLower
- upper <- priorUpper
- }
- else{
- 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(lower)){
- lower <- as.list(lower)
- }
- if(!is.list(upper)){
- upper <- as.list(upper)
- }
-
- 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)
-
- pp<-0
- if(env$h >= 1){
- qq <- 2/G
- }else{
- while(1){
- if(n*env$h^pp < 0.1) break
- pp <- pp + 1
- }
- qq <- max(pp,2/G)
- }
-
- C.temper.diff <- n_0^(2/(qq*G)-1) #this is used in pg.
- C.temper.drift <- (n_0*env$h)^(2/(qq*G)-1) #this is used in pg.
-
- mle <- qmle(yuima, "start"=start, "lower"=lower,"upper"=upper, "method"="L-BFGS-B",rcpp=rcpp)
- start <- as.list(mle at coef)
-
- integ <- function(idx.fixed=NULL,f=f,start=start,par=NULL,hessian=FALSE,upper,lower){
- if(length(idx.fixed)==0){
- intf <- hcubature(f,lowerLimit=unlist(lower),upperLimit=unlist(upper),fDim=(length(upper)+1))$integral
- }else{
- intf <- hcubature(f,lowerLimit=unlist(lower[-idx.fixed]),upperLimit=unlist(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){
-
+##::quasi-bayes function
+
+
+setGeneric("adaBayes",
+ function(yuima, start,prior,lower,upper, method="mcmc",iteration=NULL,mcmc,rate=1.0,
+ rcpp=TRUE,algorithm="randomwalk",center=NULL,sd=NULL,rho=NULL,path=FALSE
+ )
+ standardGeneric("adaBayes")
+)
+#
+# adabayes<- setClass("adabayes",contains = "mle",
+# slots = c(mcmc="list", accept_rate="list",coef = "numeirc",
+# call="call",vcov="matrix",fullcoef="numeric",model="yuima.model"))
+adabayes<- setClass("adabayes",
+ #contains = "mle",
+ slots = c(mcmc="list", accept_rate="list",coef = "numeric",call="call",vcov="matrix",fullcoef="numeric"))
+
+setMethod("adaBayes",
+ "yuima",
+ function(yuima, start,prior,lower,upper, method="mcmc",iteration=NULL,mcmc,rate=1.0,rcpp=TRUE,
+ algorithm="randomwalk",center=NULL,sd=NULL,rho=NULL,path=FALSE
+ )
+ {
+
+
+
+ if(length(iteration)>0){mcmc=iteration}
+ mean=unlist(center)
+ 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)){
+ 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)){
+ lower = append(lower,-Inf)
+ upper = append(upper,Inf)
+ }
+ }
+ # 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))
+
+
+ }
+
+ }
+
+ #20200518kaino
+ names(priorLower)<-names(pdlist)
+ names(priorUpper)<-names(pdlist)
+
+ if(missing(lower) || missing(upper)){
+ lower <- priorLower
+ upper <- priorUpper
+ }
+ else{
+ #20200518kaino
+ #if(sum(unlist(lower)<priorLower) + sum(unlist(upper)>priorUpper) > 0){
+ if(sum(unlist(lower)<priorLower[names(lower)]) + sum(unlist(upper)>priorUpper[names(upper)]) > 0){
+ yuima.stop("lower&upper of prior are out of parameter space.")
+ }
+ }
+
+ #20200518
+ #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(lower)){
+ lower <- as.list(lower)
+ }
+ if(!is.list(upper)){
+ upper <- as.list(upper)
+ }
+
+ 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")
+
+
+ fulcoef <- NULL
+
+ if(length(diff.par)>0)
+ fulcoef <- diff.par
+
+ if(length(drift.par)>0)
+ fulcoef <- c(fulcoef, drift.par)
+
+ if(length(yuima at model@parameter at common)>0){
+ fulcoef<-unique(fulcoef)
+ }
+
+
+ fixed.par <- names(fixed)
+
+ if (any(!(fixed.par %in% fulcoef)))
+ yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
+
+ nm <- names(start)
+ npar <- length(nm)
+
+ oo <- match(nm, fulcoef)
+ 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)){
+ #20200525kaino
+ #pdlist <- pdlist[order(oo)]
+ pdlist <- pdlist[names(start)]
+ }
+ nm <- names(start)
+
+ idx.diff <- match(diff.par, nm)
+ idx.drift <- match(drift.par,nm)
+ if(length(common.par)>0){
+ idx.common=match(common.par,nm)
+ idx.drift=setdiff(idx.drift,idx.common)
+ }
+ 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
+ #20200601kaino
+ if (is.CARMA(yuima)){
+ # 24/12
+ d.size <-1
+ }
+ 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)
+ #20200601kaino
+ if (is.CARMA(yuima)){
+ #24/12 If we consider a carma model,
+ # the observations are only the first column of env$X
+ # assign("X", as.matrix(onezoo(yuima)), envir=env)
+ # env$X<-as.matrix(env$X[,1])
+ # assign("deltaX", matrix(0, n-1, d.size)[,1], envir=env)
+ env$X<-as.matrix(env$X[,1])
+ env$deltaX<-as.matrix(env$deltaX[,1])
+ assign("time.obs",length(env$X),envir=env)
+ assign("p", yuima at model@info at p, envir=env)
+ assign("q", yuima at model@info at q, envir=env)
+ assign("V_inf0", matrix(diag(rep(1,env$p)),env$p,env$p), 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)
+
+ pp<-0
+ if(env$h >= 1){
+ qq <- 2/G
+ }else{
+ while(1){
+ if(n*env$h^pp < 0.1) break
+ pp <- pp + 1
+ }
+ qq <- max(pp,2/G)
+ }
+
+ C.temper.diff <- n_0^(2/(qq*G)-1) #this is used in pg.
+ C.temper.drift <- (n_0*env$h)^(2/(qq*G)-1) #this is used in pg.
+
+ mle <- qmle(yuima, "start"=start, "lower"=lower,"upper"=upper, "method"="L-BFGS-B",rcpp=rcpp)
+ start <- as.list(mle at coef)
+
+ sigma.diff=NULL
+ sigma.drift=NULL
+
+ if(length(sd)<npar){
+ sigma=mle at vcov;
+ }
+ # if(length(diff.par)>0){
+ # sigma.diff=diag(1,length(diff.par))
+ # for(i in 1:length(diff.par)){
+ # sigma.diff[i,i]=sd[[diff.par[i]]]
+ # }
+ # }
+ #
+ # if(length(drift.par)>0){
+ # sigma.drift=diag(1,length(drift.par))
+ # for(i in 1:length(drift.par)){
+ # sigma.drift[i,i]=sd[[drift.par[i]]]
+ # }
+ # }
+ #}
+
+
+ # mu.diff=NULL
+ # mu.drift=NULL
+ #
+ #
+ # if(length(mean)<npar){
+ # mu.diff=NULL
+ # mu.drift=NULL}
+ #else{
+ # if(length(diff.par)>0){
+ # for(i in 1:length(diff.par)){
+ # mu.diff[i]=mean[[diff.par[i]]]
+ # }
+ # }
+ #
+ # if(length(drift.par)>0){
+ # for(i in 1:length(drift.par)){
+ # mu.drift[i]=mean[[drift.par[i]]]
+ # }
+ # }
+ # }
+
+ ####mpcn make proposal
+ sqn<-function(x,LL){
+ vv=x%*%LL
+ zz=sum(vv*vv)
+ if(zz<0.0000001){
+ zz=0.0000001
+ }
+ return(zz)
+ }
+
+ # mpro<-function(mu,sample,low,up,rho,LL,L){
+ # d=length(mu);
+ # tmp=mu+sqrt(rho)*(sample-mu);
+ # tmp = mu+sqrt(rho)*(sample-mu);
+ # dt=(sample-mu)%*%LL%*%(sample-mu)
+ # tmp2 = 2.0/dt;
+ # while(1){
+ # prop = tmp+rnorm(d)*L*sqrt((1.0-rho)/rgamma(1,0.5*d,tmp2));
+ # if((sum(low>prop)+sum(up<prop))==0) break;
+ # }
+ # return(prop)
+ # }
+
+ make<-function(mu,sample,low,up,rho,LL,L){
+
+ d=length(mu);
+ tmp=mu+sqrt(rho)*(sample-mu);
+ dt=sqn(sample-mu,LL);
+ tmp2 = 2.0/dt;
+
+ prop = tmp+rnorm(d)%*%L*sqrt((1.0-rho)/rgamma(1,0.5*d,scale=tmp2));
+
+ # for(i in 1:100){
+ # prop = tmp+rnorm(d)*sqrt((1.0-rho)/rgamma(1,0.5*d,tmp2));
+ # if((sum(low>prop)+sum(up<prop))==0){break;}
+ # flg=flg+1;
+ # }
+ # if(flg>100){print("error")}else{
+ # return(prop)
+ # }
+ return(prop)
+ }
+
+
+ integ <- function(idx.fixed=NULL,f=f,start=start,par=NULL,hessian=FALSE,upper,lower){
+ if(length(idx.fixed)==0){
+ intf <- hcubature(f,lowerLimit=unlist(lower),upperLimit=unlist(upper),fDim=(length(upper)+1))$integral
+ }else{
+ intf <- hcubature(f,lowerLimit=unlist(lower[-idx.fixed]),upperLimit=unlist(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){
+ #vcov=vcov[nm,nm];#Song add
+ if(setequal(unlist(drift.par),unlist(diff.par))&(length(mean)>length(diff.par))){mean=mean[1:length(diff.par)]}
+ if(length(idx.fixed)==0){#only have drift or diffusion part
+ 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){
+
+ acc=0
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))
- # }
+ x_c <- mean
+ if(path){
+ X <- matrix(0,mcmc,length(mean))
+ acc=0
+ }
+
+ nm=length(mean)
+ if(length(vcov)!=(nm^2)){
+ vcov=diag(1,nm,nm)*vcov[1:nm,1:nm]
+ }
+
+ sigma=NULL
+ if(length(sd)>0&length(sd)==npar){
+ sigma=diag(sd,npar,npar);
+ if(length(idx.fixed)>0){
+ vcov=sigma[-idx.fixed,-idx.fixed]
+ }else{
+ vcov=sigma;
+ }
+ }
+
+ p_c <- p(x_c)
+ val <- f(x_c)
+ if(path) X[1,] <- 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){
- while(1){
- x_n <- rmvnorm(1,x_c,sqrt(vcov))
- if(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0) break
- }
+ if(length(mean)>1){
+
+ x_n <- rmvnorm(1,x_c,vcov)
+ #if(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0) break
}else{
- while(1){
- x_n <- rnorm(1,x_c,sqrt(vcov))
- if(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0) break
- }
+ x_n <- rnorm(1,x_c,sqrt(vcov))
+ #if(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0) break
}
-
- #if(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0){
- #q_n <- dnorm(x_n,x_c,sqrt(vcov))
- p_n <- p(x_n)
- #q_c <- dnorm(x_c,x_n,sqrt(vcov))
- #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(sum(x_n<lowerLimit)==0 & sum(x_n>upperLimit)==0){
+ p_n <- p(x_n)
+ u <- log(runif(1))
a <- p_n-p_c
- #asd[(i+1),2:3] <- c(p_c, p_n)
- if(u<a){
- p_c <- p_n
- #q_c <- q_n
- x_c <- x_n
- }
+ if(u<a){
+ p_c <- p_n
+ # q_c <- q_n
+ x_c <- x_n
+ acc=acc+1
+ }
+ }
+
+ if(path) X[i+1,] <- x_c
#}
- val <- val+f(x_c)
- }
- return(unlist(val/mcmc))
- }
- else if(tolower(algorithm)=="mpcn"){ #MpCN
- x_n <- mean
- val <- mean
- logLik_old <- p(x_n)+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(prop)+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_0,length(diff.par)),rep(1 / (n_0 * 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){
+ val <- val+f(x_c)
+ }
+ if(path){
+ return(list(val=unlist(val/mcmc),X=X,acc=acc/mcmc))
+ }else{
+ return(list(val=unlist(val/mcmc)))
+ }
+ }
+ else if(tolower(algorithm)=="mpcn"){ #MpCN
+ if(path) X <- matrix(0,mcmc,length(mean))
+ if((sum(idx.diff==idx.fixed)==0)){
+ if(length(center)>0){
+ if(length(idx.fixed)>0){
+ mean =unlist(center[-idx.fixed])
+ }else{
+ mean =unlist(center) ##drift or diffusion parameter only
+ }
+ }
+ }
+
+ # if((sum(idx.diff==idx.fixed)>0)){
+ # if(length(center)>0){
+ # mean =unlist(center[-idx.fixed])
+ # }
+ # }
+ x_n <- mean
+ val <- mean
+
+ L=diag(1,length(mean));LL=L;
+ nm=length(mean)
+ LL=vcov;
+ if(length(vcov)!=(nm^2)){
+ LL=diag(1,nm,nm)
+ }
+ if(length(sd)>0&length(sd)==npar){
+ sigma=diag(sd,npar,npar);
+ if(length(idx.fixed)>0){
+ LL=sigma[-idx.fixed,-idx.fixed]
+ }else{
+ LL=sigma;
+ }
+ }
+
+
+ # if(length(pre_conditionnal)==1){
+ # LL=t(chol(solve(vcov)));L=chol(vcov);
+ # }
+ if(length(rho)==0){rho=0.8}
+
+ logLik_old <- p(x_n)+0.5*length(mean)*log(sqn(x_n-mean,LL))
+ if(path) X[1,] <- x_n
+
+
+ for(i in 1:(mcmc-1)){
+ #browser()
+ prop <- make(mean,x_n,unlist(lowerLimit),unlist(upperLimit),rho,LL,L)
+ if((sum(prop<lowerLimit)+sum(prop>upperLimit))==0){
+ logLik_new <- p(prop)+0.5*length(mean)*log(sqn(prop-mean,LL))
+ u <- log(runif(1))
+ if( logLik_new-logLik_old > u){
+ x_n <- prop
+ logLik_old <- logLik_new
+ acc=acc+1
+ }
+ }
+ if(path) X[i+1,] <- x_n
+ val <- val+f(x_n)
+ }
+ #return(unlist(val/mcmc))
+ if(path){
+ return(list(val=unlist(val/mcmc),X=X,acc=acc/mcmc))
+ }else{
+ return(list(val=unlist(val/mcmc)))
+ }
+ }
+ }
+
+ #print(mle at coef)
+ flagNotPosDif <- 0
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 739
More information about the Yuima-commits
mailing list