[Yuima-commits] r96 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 11 12:02:36 CEST 2010
Author: iacus
Date: 2010-07-11 12:02:36 +0200 (Sun, 11 Jul 2010)
New Revision: 96
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/adaBayes.R
pkg/yuima/R/asymptotic_term.R
pkg/yuima/R/qmle.R
pkg/yuima/man/quasi-likelihood.Rd
Log:
qmle with joint estimation
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-07-11 08:41:16 UTC (rev 95)
+++ pkg/yuima/DESCRIPTION 2010-07-11 10:02:36 UTC (rev 96)
@@ -1,10 +1,10 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.0.96
+Version: 0.0.97
Date: 2010-07-11
Depends: methods, zoo, stats4, utils
-Suggests: adapt, mvtnorm
+Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Description: The YUIMA Project for Simulation and Inference
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2010-07-11 08:41:16 UTC (rev 95)
+++ pkg/yuima/NEWS 2010-07-11 10:02:36 UTC (rev 96)
@@ -1,3 +1,12 @@
+Changes in Version 0.0.97
+
+ o no longer using adapt, we now use cubature which is GPL
+ and available on CRAN. Package cubature and adapt use
+ the same algorithm, with cubature being more general and
+ up to date
+
+ o adaBayes and asymptotic_term use not cubature
+
Changes in Version 0.0.94
o implemented two stage estimation in qmle
Modified: pkg/yuima/R/adaBayes.R
===================================================================
--- pkg/yuima/R/adaBayes.R 2010-07-11 08:41:16 UTC (rev 95)
+++ pkg/yuima/R/adaBayes.R 2010-07-11 10:02:36 UTC (rev 96)
@@ -447,7 +447,8 @@
}
if(method=="nomcmc"){
- require(adapt)
+ require(cubature)
+# 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
@@ -475,7 +476,9 @@
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
+
+# denomvalue <- adapt(length(unlist(param[slot(yuima at model@parameter,term)])),lower=mdom[1,],upper=mdom[2,],functn=denominator)$value
+ denomvalue <- adaptIntegrate(denominator, lower=mdom[1,],upper=mdom[2,])$integral
}
## END denominator calculation
@@ -500,7 +503,8 @@
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
+# numevalue <- adapt(length(param[slot(yuima at model@parameter,term)]),lower=mdom[1,],upper=mdom[2,],functn=numerator)$value
+ numevalue <- adaptIntegrate(numerator, lower=mdom[1,],upper=mdom[2,])$integral
}
unlistparam[i] <- numevalue/denomvalue
}
Modified: pkg/yuima/R/asymptotic_term.R
===================================================================
--- pkg/yuima/R/asymptotic_term.R 2010-07-11 08:41:16 UTC (rev 95)
+++ pkg/yuima/R/asymptotic_term.R 2010-07-11 10:02:36 UTC (rev 96)
@@ -583,13 +583,15 @@
## integrate
if( k.size ==1){ # use 'integrate' if k=1
- tmp <- integrate(gz_pi02,-Inf,Inf)
+ tmp <- integrate(gz_pi02,-Inf,Inf)$value
}else if( 2 <= k.size || k.size <= 20 ){ # use library 'adapt' to solve multi-dimentional integration
max <- 10 * lambda.max
min <- -10 * lambda.max
L <- (max - min)
- if(require(adapt)){
- tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi0)
+# if(require(adapt)){
+ if(require(cubature)){
+#tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi0)$value
+ tmp <- adaptIntegrate(gz_pi0, lower=rep(min,k.size),upper=rep(max,k.size))$integral
} else {
tmp <- NA
}
@@ -1260,12 +1262,14 @@
}
## integrate
if( k.size ==1){ # use 'integrate()'
- tmp <- integrate(gz_pi1,-Inf,Inf)
+ tmp <- integrate(gz_pi1,-Inf,Inf)$value
}else if( 2 <= k.size || k.size <= 20 ){ # use 'adapt()' to solve multi-dim integration8
max <- 10*lambda.max
min <- -10*lambda.max
- if(require(adapt)){
- tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi1)
+# if(require(adapt)){
+ if(require(cubature)){
+# tmp <- adapt(ndim=k.size,lower=rep(min,k.size),upper=rep(max,k.size),functn=gz_pi1)$value
+ tmp <- adaptIntegrate(gz_pi1,lower=rep(min,k.size),upper=rep(max,k.size))$integral
} else {
tmp <- NA
}
@@ -1479,6 +1483,6 @@
yuima.warn("Calculating d1 term ...")
d1 <- get.d1.term()
yuima.warn("Done\n")
- terms <- list(d0=d0$value, d1=d1$value)
+ terms <- list(d0=d0, d1=d1)
return(terms)
})
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2010-07-11 08:41:16 UTC (rev 95)
+++ pkg/yuima/R/qmle.R 2010-07-11 10:02:36 UTC (rev 96)
@@ -13,24 +13,19 @@
r.size <- yuima at model@noise.number
d.size <- yuima at model@equation.number
modelstate <- yuima at model@state.variable
-# modelpara <- yuima at model@parameter at drift
DRIFT <- yuima at model@drift
n <- length(yuima)[1]
drift <- matrix(0,n,d.size)
-# X <- as.matrix(onezoo(yuima))
for(i in 1:length(theta)){
assign(names(theta)[i],theta[[i]])
}
for(t in 1:n){
-# Xt <- X[t,]
- for(d in 1:d.size){
-# assign(modelstate[d],Xt[d])
+ for(d in 1:d.size)
assign(modelstate[d], env$X[t,d])
- }
- for(d in 1:d.size){
+# do not collapse the two for loops
+ for(d in 1:d.size)
drift[t,d] <- eval(DRIFT[d])
- }
}
return(drift)
}
@@ -44,20 +39,18 @@
DIFFUSION <- yuima at model@diffusion
n <- length(yuima)[1]
diff <- array(0, dim=c(d.size, r.size, n))
-# X <- as.matrix(onezoo(yuima))
for(i in 1:length(theta)){
assign(names(theta)[i],theta[[i]])
}
for(r in 1:r.size){
for(t in 1:n){
-# Xt <- X[t, ]
- for(d in 1:d.size){
+ for(d in 1:d.size)
assign(modelstate[d], env$X[t,d])
- }
- for(d in 1:d.size){
+# do not collapse the two for loops
+ for(d in 1:d.size)
diff[d, r, t] <- eval(DIFFUSION[[d]][r])
- }
+
}
}
return(diff)
@@ -65,7 +58,7 @@
-
+## take from original Hino-san code
##::calculate diffusion%*%t(diffusion) matrix
calc.B <- function(diff){
d.size <- dim(diff)[1]
@@ -89,7 +82,7 @@
qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
- lower, upper, ...){
+ lower, upper, joint=FALSE, ...){
call <- match.call()
@@ -109,7 +102,7 @@
measure.par <- yuima at model@parameter at measure
common.par <- yuima at model@parameter at common
- JointOptim <- FALSE
+ JointOptim <- joint
if(length(common.par)>0){
JointOptim <- TRUE
yuima.warn("Drift and diffusion parameters must be different. Doing
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-07-11 08:41:16 UTC (rev 95)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-07-11 10:02:36 UTC (rev 96)
@@ -17,7 +17,7 @@
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)
-qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper, ...)
+qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper, joint=FALSE, ...)
quasilogl(yuima, param, print=FALSE)
}
\arguments{
@@ -36,6 +36,7 @@
\item{upper}{a named list for specifying upper bounds of parameters}
\item{start}{initial values to be passed to the optimizer.}
\item{fixed}{for conditional (quasi)maximum likelihood estimation.}
+ \item{joint}{perform joint estimation or two stage estimation? by default \code{joint=FALSE}.}
\item{...}{passed to \code{\link{optim}} method. See Examples.}
}
\details{
@@ -93,16 +94,17 @@
print("ML estimator")
opt2 at 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)
-##)
-##print("True param")
-##print("theta2 = .3, theta1 = .1")
-##print("ML estimator")
-#opt at coef
+## perform joint estimation? Non-optimal, just for didactic purposes
+system.time(
+opt3 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=.05,theta2=.05),
+ upper=list(theta1=.5,theta2=.5), method="L-BFGS-B", joint=TRUE)
+)
+print("True param")
+print("theta2 = .3, theta1 = .1")
+print("ML estimator")
+opt3 at coef
+
system.time(
opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newtoon")
)
More information about the Yuima-commits
mailing list