[Yuima-commits] r82 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 25 10:48:36 CEST 2010
Author: kamatani
Date: 2010-06-25 10:48:35 +0200 (Fri, 25 Jun 2010)
New Revision: 82
Modified:
pkg/yuima/R/adaBayes.R
Log:
input method is updated as qmle
Modified: pkg/yuima/R/adaBayes.R
===================================================================
--- pkg/yuima/R/adaBayes.R 2010-06-23 00:29:17 UTC (rev 81)
+++ pkg/yuima/R/adaBayes.R 2010-06-25 08:48:35 UTC (rev 82)
@@ -1,164 +1,690 @@
-## Bayes estimator to calculate theta1 and theta2
-## arguments
-## sdemod : sdeModel object
-## h : sampling interval
-## prior1: prior distribution of theta1 (returned value is scalar)
-## prior2: prior distribution of theta2 (returned value is scalar)
-## domain1,domain2 : domain matrixes (or vectors) of theta1 and theta2 (length(theta1) x 2 matrix ,length(theta2) x 2 matrix)
-## ex.) theta1 <- c("theta11","theta12") , theta2 = c("theta21","theta22","theta23")
-## theta1.lower <- c(0,0) , theta1.upper <- c(1,2)
-## theta2.lower <- c(0.5,0.2,0.7) , theta2.upper <- c(10,5,15)
-## domain1 <- matrix(c(theta1.lower , theta1.upper ),2,2)
-## domain2 <- matrix(c(theta2.lower , theta2.upper ),3,2)
-## theta1.ans , theta2.ans : offset vector (or scalar) to calculate
+##::quasi-likelihood with prior
+##
+##
+
+## rmnorm
+rmnorm <- function (n, mean, cov=diag(length(mean)))
+{
+ if(length(mean)==1){
+ return(rnorm(n,mean,sqrt(cov)))
+ }else{
+ return(rmvnorm(n,mean,cov))
+ }
+}
+
+## dmnorm
+dmnorm <- function (x, mean, cov=diag(length(mean)))
+{
+ if(length(mean)==1){
+ return(dnorm(x,mean,sqrt(cov)))
+ }else{
+ return(dmvnorm(x,mean,cov))
+ }
+}
+
+## ml.ql by newton method.
+newton.ml.qlb <- 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)
+ }
+
+## Temporary dH2_theta1
+ d2Htemp.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
+ temp <- dGi.theta1(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...)
+ ret <- ret + temp%*%t(temp)
+ }
+ 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)
+ }
+
+## Temp dH_theta2
+ d2Htemp.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
+ temp <- dGi.theta2(theta1, theta2, h, Xt.iMinus1, delta.i.X, B, a, ...)
+ ret <- ret + temp%*%t(temp)
+#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 <- d2Htemp.theta2(theta1, theta2, h, yuima, B, a, ...)
+ theta2 <- as.vector(theta2 - solve(d2Htheta2) %*% dHtheta2)
+ dHtheta1 <- dH.theta1(theta1, theta2, h, yuima, B, a, ...)
+ d2Htheta1 <- d2Htemp.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, ...))
+ Sigma1temp <- - solve(d2Htemp.theta1(theta1, theta2, h, yuima, B, a, ...))
+ Sigma2 <- - solve(d2H.theta2(theta1, theta2, h, yuima, B, a, ...))
+ Sigma2temp <- - solve(d2Htemp.theta2(theta1, theta2, h, yuima, B, a, ...))
+ hessian <- 1 ## Not Implemented yet!
+##END
+##temp
+ return(list(theta1.new=theta1, theta2.new=theta2, Sigma1=Sigma1, Sigma1temp=Sigma1temp, Sigma2=Sigma2,Sigma2temp=Sigma2temp, hessian=hessian))
+
+}
+
+##::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.
+
+##::quasi-bayes function
+
+##::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("adaBayes",
- function(yuima,h,prior2,prior1,domain2,domain1,theta2.offset=0,theta1.offset=0)
+ function(yuima, print=FALSE, init, prior,propose,n.iter=100,lower,upper,n.burnin,method="nomcmc",mhtype="independent")
standardGeneric("adaBayes")
)
-setMethod("adaBayes",signature("yuima"),function(yuima,h,prior2,prior1,domain2 ,domain1 ,theta2.offset=0,theta1.offset=0){
- ## Hn function
- Hn <- function(yuima,h,theta1,theta2){
- d.size <- yuima at model@equation.number
- n.size <- dim(as.matrix(onezoo(yuima)))[1]
- tmp <- ((n.size*d.size/2) * log( (2*pi*h) ) + ql(yuima,theta2=theta2,theta1=theta1,h=h))
- return(tmp)
- }
-
- ## rHn function
- rHn <- function(yuima,h,theta1,theta2,theta2.offset,theta1.offset){
- d.size <- yuima at model@equation.number
- n.size <- dim(as.matrix(onezoo(yuima)))[1]
- tmp <- ((n.size*d.size/2) * log( (2*pi*h) ) + rql(yuima,theta2=theta2,theta1=theta1,theta2.offset,theta1.offset,h=h))
- return(tmp)
- }
-
- ## Bayes estimator for theta1
- bayes.theta1.tilda <- function(yuima,prior1=prior1){
- rHn.tmp <- function(yuima,h,theta1,theta2,theta2.offset,theta1.offset){
- tmp <- numeric( length(theta1) )
- for(i in 1:length(theta1)){
- rHn1 <- rHn( yuima , h , theta1[i] , theta2 ,theta2.offset,theta1.offset)
- tmp[i] <- exp( rHn1 )
- }
- if( is.infinite(sum(tmp)) != 0 ){
- stop("Hn is too big.\n Please check theta2.offset and theta1.offset")
- }
- return(tmp)
- }
- ## if theta1 is a vector , prior1 returns vector compliant with each theta1 value
- prior1.tmp <- function(theta1){
- tmp <- numeric(length(theta1))
- for( i in 1:length(theta1) ){
- tmp[i] <- prior1(theta1[i],domain1)
- }
- return(tmp)
- }
- ## numerator function
- Numerator.tmp <- function(theta1){
- tmp <- theta1 *rHn.tmp(yuima=yuima,h=h,theta1,theta2=1,theta2.offset=theta2.offset,theta1.offset=theta1.offset) * prior1.tmp(theta1)
- return(tmp)
- }
- ## denominator function
- Denominator.tmp <- function(theta1){
- tmp <- rHn.tmp(yuima=yuima,h=h,theta1,theta2=1,theta2.offset=theta2.offset,theta1.offset=theta1.offset) * prior1.tmp(theta1)
- return(tmp)
- }
- numerator <- integrate(Numerator.tmp,domain1[1,1],domain1[1,2])
- denominator <- integrate(Denominator.tmp,domain1[1,1],domain1[1,2])
- if(numerator$value == 0 || denominator$value == 0){
- stop("Quasi-likelihood is too small.\n Please check domain1 , theta2.offset and theta1.offset")
- }
- return( numerator$value/denominator$value )
- }
-
- ##Bayes estimator for theta2
- bayes.theta2.tilda <- function(yuima,prior2=prior2){
- Hn.tmp <- function(yuima,h,theta1,theta2){
- tmp <- numeric(length(theta2))
- for(i in 1:length(theta2)){
- Hn1 <- Hn( yuima , h , theta1 , theta2[i] )
- tmp[i] <- exp( Hn1 -Hn2 )
- }
- if( is.infinite(sum(tmp)) != 0 ){
- stop("Hn is too big.\n Please check theta2.offset and theta1.offset")
- }
- return(tmp)
- }
- ## if theta2 is a vector , prior1 returns vector compliant with each theta2 value
- prior2.tmp <- function(theta2){
- tmp <- numeric(length(theta2))
- for( i in 1:length(theta2) ){
- tmp[i] <- prior2(theta2[i],domain2)
- }
- return(tmp)
- }
- ## numerator function
- Numerator.tmp <- function(theta2){
- tmp <- theta2 * Hn.tmp(yuima=yuima,h=h,theta1=theta1.tilda,theta2) * prior2.tmp(theta2)
- return(tmp)
- }
- ## denominator function
- Denominator.tmp <- function(theta2){
- tmp <- Hn.tmp(yuima=yuima,h=h,theta1=theta1.tilda,theta2) * prior2.tmp(theta2)
- return(tmp)
- }
- numerator <- integrate(Numerator.tmp,domain2[1,1],domain2[1,2])
- denominator <- integrate(Denominator.tmp,domain2[1,1],domain2[1,2])
- if(numerator$value == 0 || denominator$value == 0){
- stop("Quasi-likelihood is too small.\n Please check domain2 , theta2.offset and theta1.offset")
- }
- return( numerator$value/denominator$value )
- }
- ## Start error Check ##
- ## check yuima
- if( missing(yuima) ){
- stop("yuima is missing.")
- }
- ## check number of theta
- if( length(yuima at model@diffusion) != 1 || length(yuima at model@drift) != 1){
- stop("This version do not support multi-dimention parameters.")
- }else{
- if( is.vector(domain2) ) domain2 <- t(as.matrix(domain2))
- if( is.vector(domain1) ) domain1 <- t(as.matrix(domain1))
- }
- ## check domain
- if( !is.matrix(domain2) || !is.matrix(domain1) ){
- stop("domain2 or domain1 is not a matrix.")
- }else if( length(yuima at model@diffusion) != nrow(domain1) || length(yuima at model@drift) != nrow(domain2)){
- stop("dimention of domain1 or domain2 is different from length with respect to each theta.")
- }else if( ncol(domain1) != 2 || ncol(domain2) != 2){
- stop("ncol(domain1) or ncol(domain2) is irrelevance.")
- }
- ## check X
- ##if(is.matrix(yuima at data@original.data) == FALSE){##original.data
- ## stop("data is not a matrix. Please check again.")
- ##}else if
- if( dim(as.matrix(onezoo(yuima)))[1] < 2 ){
- stop("length of data is too short.")
- }
- ## check h
- if( (!is.double(h) && !is.integer(h)) || length(h) != 1){
- stop("'h' is not correct type.")
- }else if( h <= 0 ){
- stop("'h' is negative or zero. Please set 'h' to be a positive value.")
- }
- ## check theta offset
- if( !is.null(theta1.offset) && !is.null(theta2.offset) ){
- if(!is.vector(theta2.offset) || !is.vector(theta1.offset)){
- stop("theta2.offset or theta1.offset is not a vector.")
- }
- if(length(yuima at model@drift) != length(theta2.offset) || length(yuima at model@diffusion) != length(theta1.offset)){
- stop("length of theta2.offset or theta1.offset is different from length with respect to each theta")
- }
- Hn2 <- Hn( yuima , h , theta1 = theta1.offset , theta2 = theta2.offset )
- }else{
- Hn2 <- 0
- }
- ## Finish error check ##
-
- ## Culculate theta1.tilda and theta2.tilda
- theta1.tilda <- bayes.theta1.tilda(yuima,prior1=prior1)
- theta2.tilda <- bayes.theta2.tilda(yuima,prior2=prior2)
- return(c(theta2.tilda,theta1.tilda))
-}
-)
+setMethod("adaBayes", "yuima",
+ function(yuima, print=FALSE, init, prior,propose,n.iter=100,lower,upper,n.burnin,method="nomcmc",mhtype="independent"){
+ if( missing(yuima)){
+ cat("\nyuima object is missing.\n")
+ return(NULL)
+ }
+
+ if(length(match(yuima at model@parameter at drift ,names(init),nomatch=0))!=length(yuima at model@parameter at drift)){
+ cat("\ndrift parameters in yuima model and init do not match.\n")
+ return(NULL)
+ }
+ if(length(match(yuima at model@parameter at diffusion ,names(init),nomatch=0))!=length(yuima at model@parameter at diffusion)){
+ cat("\ndiffusion parameters in yuima model and init do not match.\n")
+ return(NULL)
+ }
+
+ ## BEGIN burnin handling
+ if(missing(n.burnin)){
+ n.burnin <- n.iter
+ }
+ ## END burnin handling
+
+ h <- deltat(yuima at data@zoo.data[[1]])
+
+ ## BEGIN Prior construction
+ liprior<- function(prior,term){
+ mvec <- numeric(0); mdom <- numeric(0)
+ for(i in 1:length(slot(yuima at model@parameter,term))){
+ if(prior[slot(yuima at model@parameter,term)[i]][[1]]$measure.type=="density"){
+ mvec <- append(mvec,prior[slot(yuima at model@parameter,term)[i]][[1]]$density)
+ if(is.null(prior[slot(yuima at model@parameter,term)[i]][[1]]$domain)){
+ mdom <- cbind(mdom,c(-Inf,Inf))
+ }else{
+ mdom <- cbind(mdom,prior[slot(yuima at model@parameter,term)[i]][[1]]$domain)
+ }
+ }
+ }
+ mdensity <- function(param){
+ res <- 1
+ for(i in 1:length(slot(yuima at model@parameter,term))){
+ res <- res*mvec[i][[1]](param[slot(yuima at model@parameter,term)[i]])
+ }
+ return(res)
+
+ }
+ return(list(mdom=mdom,mdensity=mdensity))
+ }
+ ## END Prior construction
+
+ n <- length(yuima)[1]
+
+ ilparam <- function(liparam){
+ return(list("theta2"=unlist(liparam[slot(yuima at model@parameter,"drift")]),"theta1"=unlist(liparam[slot(yuima at model@parameter,"diffusion")])))
+ }
+ liparam <- function(ilparam){
+ return(relist(unlist(ilparam),init))
+ }
+
+ if(method=="nomcmc"){
+ require(adapt)
+ ## BEGIN numerical integration function
+ nintegral <- function(term,param=init,prior=prior,lower=lower,upper=upper,print=FALSE){
+ lip <- liprior(prior,term); mdom <- lip$mdom; mdensity <- lip$mdensity
+
+ ## BEGIN denominator calculation
+ nparam <- param
+ denominator <- function(subparam){
+ if(is.null(dim(subparam))){
+ matparam <- matrix(subparam,ncol=length(unlist(param[slot(yuima at model@parameter,term)])))
+ }else{
+ matparam <- subparam
+ }
+
+ nparam <- param
+
+ qvec <- numeric(dim(matparam)[1])
+ for(k in 1:dim(matparam)[1]){
+ for(j in 1:(length(slot(yuima at model@parameter,term)))){
+ nparam[slot(yuima at model@parameter,term)][j] <- matparam[k,][grep(slot(yuima at model@parameter,term)[j],names(unlist(param[slot(yuima at model@parameter,term)])))]
+ }
+ qvec[k] <- exp(-negquasilik(yuima,param=nparam,print=print)+log(mdensity(nparam))+negquasilik(yuima,param=param,print=print)-log(mdensity(param)))
+ }
+ return(qvec)
+ }
+ if(length(slot(yuima at model@parameter,term))==1){
+ denomvalue <- integrate(denominator,mdom[1,1],mdom[2,1])$value
+ }else{
+ denomvalue <- adapt(length(unlist(param[slot(yuima at model@parameter,term)])),lower=mdom[1,],upper=mdom[2,],functn=denominator)$value
+ }
+ ## END denominator calculation
+
+ ## BEGIN numerator calculation
+ unlistparam <- numeric(length(param[slot(yuima at model@parameter,term)]))
+ for(i in 1:length(param[slot(yuima at model@parameter,term)])){
+ numevalue <- numeric(1)
+ nparam <- param
+ numerator <- function(subparam){
+ if(is.null(dim(subparam))){
+ matparam <- matrix(subparam,ncol=length(unlist(param[slot(yuima at model@parameter,term)])))
+ }else{
+ matparam <- subparam
+ }
+ qvec <- numeric(dim(matparam)[1])
+ for(k in 1:(dim(matparam)[1])){
+ nparam[slot(yuima at model@parameter,term)] <- matparam[k,]
+ qvec[k] <- matparam[k,i]*exp(-negquasilik(yuima,param=nparam,print=print)+log(mdensity(nparam))+negquasilik(yuima,param=param,print=print)-log(mdensity(param)))
+ }
+ return(qvec)
+ }
+ if(length(slot(yuima at model@parameter,term))==1){
+ numevalue <- integrate(numerator,mdom[1,1],mdom[2,1])$value
+ }else{
+ numevalue <- adapt(length(param[slot(yuima at model@parameter,term)]),lower=mdom[1,],upper=mdom[2,],functn=numerator)$value
+ }
+ unlistparam[i] <- numevalue/denomvalue
+ }
+ ## END numerator calculation
+ newparam <- param
+ newparam[slot(yuima at model@parameter,term)] <- relist(unlistparam,param[slot(yuima at model@parameter,term)])
+
+ return(newparam)
+ }
+ ## END numerical integration function
+
+ ## BEGIN numerical integration procedure
+ param <- nintegral("diffusion",param=init,prior=prior,lower=lower,upper=upper)
+ param <- nintegral("drift",param=param,prior=prior,lower=lower,upper=upper)
+ ## END numerical integration procedure
+
+ }else if(method=="mcmc"){
+ if(mhtype=="RW"){
+ ## BEGIN propose construction
+ if(missing(propose)){
+ rpropose1 <- function(n,mean,varcov){
+ sd <- sqrt(varcov)
+ d <- length(init[slot(yuima at model@parameter,"diffusion")])
+ return(matrix(data=rnorm(d*n,mean,sd),nrow=n))
+ }
+ rpropose2 <- function(n,mean,varcov){
+ sd <- sqrt(varcov)
+ d <- length(init[slot(yuima at model@parameter,"drift")])
+ return(matrix(data=rnorm(d*n,mean,sd),nrow=n))
+ }
+
+ propose.param1 <- list(mean=0,varcov=1/n)
+ propose.param2 <- list(mean=0,varcov=1/(n*h))
+ propose1 <- list(rpropose=rpropose1,propose.param=propose.param1)
+ propose2 <- list(rpropose=rpropose2,propose.param=propose.param2)
+ propose <- list(propose2 = propose2,propose1 = propose1)
+ }
+ ## END propose construction
+
+ ## BEGIN mh RW-algorithm function
+ mhalgorithm <- function(j,length=n.iter,param=init,prior=prior,print=FALSE){
+ if(prior[[j]]$measure.type=="density"){
+ mdensity <- prior[[j]]$density
+ if(is.null(prior[[j]]$domain)){
+ mdom <- c(-Inf,Inf)
+ }else{
+ mdom <- prior[[j]]$domain
+ }
+ }
+
+ rpropose <- propose[[j]][[1]]
+ propose.param <- propose[[j]][[2]]
+ d <- length(param[[j]])
+
+ dh <- matrix(t(rpropose(n.iter,propose.param[[1]],propose.param[[2]])),nrow=d)
+
+ u <- runif(n.iter)
+
+ pparam <- param
+ pql <- -negquasilik(yuima,param=pparam,print=print)+log(mdensity(pparam[[j]]))
+ nparam <- param
+ mhestimator <- 0
+
+ for(i in 1:n.iter){
+ nparam[[j]] <- pparam[[j]] + dh[,i]
+ nql <- -negquasilik(yuima,param=nparam,print=print)+log(mdensity(nparam[[j]]))
+ alpha <- exp(nql-pql)
+ if(is.na(alpha)){alpha <- -Inf}
+ pparam[[j]] <- (alpha>u[i])*(nparam[[j]]-pparam[[j]])+pparam[[j]]
+ if(alpha>u[i]){ pql <- nql}
+ mhestimator <- mhestimator + pparam[[j]]
+ }
+
+ return(mhestimator/n.iter)
+ }
+ ## END mh RW-algorithm function
+ }
+
+ if(mhtype=="independent"){
+ require(mvtnorm)
+ ## propose distribution construction
+ if(missing(propose)){
+ rpropose1 <- function(n,mean,varcov){
+ d <- length(init[slot(yuima at model@parameter,"diffusion")])
+
+ if(d==1) varcov <- as.numeric(varcov)
+
+ return( t( as.matrix(rmnorm(n,mean,varcov)) ) )
+ }
+ rpropose2 <- function(n,mean,varcov){
+ d <- length(init[slot(yuima at model@parameter,"drift")])
+
+ if(d==1) varcov <- as.numeric(varcov)
+
+ return( t( as.matrix(rmnorm(n,mean,varcov)) ) )
+ }
+ dpropose <- function(x,mean,varcov){
+ n <- dim(as.matrix(x))[2]
+ d <- dim(as.matrix(x))[1]
+
+ if(n==1 & d==1) x <- matrix(x)
+
+ if(!is.matrix(x)){
+ x <- as.matrix(x)
+ }
+ dvector <- numeric(n)
+
+ for(i in 1:n){
+ dvector[i] <- dmnorm(x[,i],mean,varcov)
+ }
+ return(dvector)
+ }
+ }
+
+ ## BEGIN mh Independent type algorithm function
+ mhalgorithm <- function(j,length=n.iter,param=init,prior=prior,print=FALSE){
+ if(j==2){term <- "diffusion"}else{term <- "drift"}
+ lip <- liprior(prior,term); mdom <- lip$mdom; mdensity <- lip$mdensity
+ if(sum(c("theta2","theta1") %in% names(param))==2){
+ unlistparam <- param
+ }else{
+ unlistparam <- ilparam(param)
+ }
+ newton <- newton.ml.qlb(yuima, theta2=unlistparam[[1]],theta1=unlistparam[[2]], h, iteration=1)
+ if(print){ cat("newton result\n"); print(newton) }
+
+ if(j==1){##estimate theta2
+ newton.tmp <- newton.ml.qlb(yuima, theta2=unlist(newton$theta2.new), theta1=unlistparam[[2]], h, param.only=TRUE)
+ propose.param1 <- list(mean=newton.tmp$theta1.new,varcov=newton.tmp$Sigma1temp)
+ propose.param2 <- list(mean=newton.tmp$theta2.new,varcov=newton.tmp$Sigma2temp)
+ if(length(prior)!=1) prior <- prior$prior.theta2
+
+ }else if(j==2){ ##estimate theta1
+ newton.tmp <- newton.ml.qlb(yuima, theta2=unlistparam[[1]], theta1=newton$theta1.new, h, param.only=TRUE)
+ propose.param1 <- list(mean=newton.tmp$theta1.new,varcov=newton.tmp$Sigma1temp)
+ propose.param2 <- list(mean=newton.tmp$theta2.new,varcov=newton.tmp$Sigma2temp)
+ if(length(prior)!=1) prior <- prior$prior.theta1
+ }
+ propose1 <- list(rpropose=rpropose1, dpropose=dpropose, propose.param=propose.param1)
+ propose2 <- list(rpropose=rpropose2, dpropose=dpropose, propose.param=propose.param2)
+ propose <- list(propose2=propose2, propose1=propose1)
+ param <- list(theta2=newton$theta2.new, theta1=newton$theta1.new)
+
+ rpropose <- propose[[j]][[1]]
+ dpropose <- propose[[j]][[2]]
+ propose.param <- propose[[j]][[3]]
+ d <- length(param[[j]])
+
+ state <- rpropose(n.iter,propose.param[[1]],propose.param[[2]])
+ tmp.param <- param
+ weight <- numeric(n.iter)
+
+ for(i in 1:n.iter){
+ tmp.param[[j]] <- state[,i]
+ weight[i] <- exp(-negquasilik(yuima,param=liparam(tmp.param),print=print))*mdensity(liparam(tmp.param))/
+ dpropose(state[,i],propose.param[[1]],propose.param[[2]])
+ }
+
+ cstate <- param[[j]]
+ cweight <- exp(-negquasilik(yuima,param=liparam(param),print=print))*mdensity(liparam(tmp.param))/
+ dpropose(param[[j]],propose.param[[1]],as.matrix(propose.param[[2]]))
+
+ u <- runif(n.iter)
+
+ mhestimator <- 0
+
+ for(i in 1:n.iter){
+ alpha <- weight[i]/cweight
+ if(is.na(alpha)){alpha <- -Inf}
+ if(alpha>u[i]){ cweight <- weight[i]; cstate <- state[,i]}
+ mhestimator <- mhestimator + state[,i]
+ }
+ return(mhestimator/n.iter)
+ }
+ }
+
+ ## BEGIN mh procedure
+ ##theta1 estim.
+ param <- liparam(list(init[slot(yuima at model@parameter,"drift")],mhalgorithm(2,length=n.burnin,param=init,prior=prior, print=print)))
+ param <- liparam(list(param[slot(yuima at model@parameter,"drift")],mhalgorithm(2,length=n.iter,param=param,prior=prior, print=print)))
+ ##theta2 estim.
+ param <- liparam(list(list(mhalgorithm(1,length=n.burnin,param=param,prior=prior, print=print),param[slot(yuima at model@parameter,"diffusion")])))
+ param <- liparam(list(mhalgorithm(1,length=n.iter,param=param,prior=prior, print=print),param[slot(yuima at model@parameter,"diffusion")]))
+ ## END mh procedure
+ }
+ return(liparam(param))
+ })
More information about the Yuima-commits
mailing list