[Yuima-commits] r114 - in pkg/yuima: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 20 16:54:54 CEST 2010
Author: iacus
Date: 2010-07-20 16:54:54 +0200 (Tue, 20 Jul 2010)
New Revision: 114
Modified:
pkg/yuima/R/CPoint.R
pkg/yuima/R/subsampling.R
pkg/yuima/R/yuima.sampling.R
pkg/yuima/man/CPoint.Rd
Log:
qmleR and qmleL make use of qmle
Modified: pkg/yuima/R/CPoint.R
===================================================================
--- pkg/yuima/R/CPoint.R 2010-07-20 12:29:29 UTC (rev 113)
+++ pkg/yuima/R/CPoint.R 2010-07-20 14:54:54 UTC (rev 114)
@@ -1,208 +1,38 @@
-qmleL <- function(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...){
+qmleL <- function(yuima, t, ...){
- 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
-
- if(!is.list(start))
- yuima.stop("Argument 'start' must be of list type.")
-
- fullcoef <- diff.par
- npar <- length(fullcoef)
- 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)
-
- 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]
-
times <- time(yuima at data@zoo.data[[1]])
minT <- as.numeric(times[1])
maxT <- as.numeric(times[length(times)])
- if(missing(t))
+ if(missing(t) )
t <- mean(c(minT,maxT))
- k <- max(which(times <=t),na.rm=TRUE)[1]
-
- if(k<2)
- k <- 2
- if(k>n-2)
- k <- n-2
- env <- new.env()
-# assign("X", as.matrix(onezoo(yuima)), env=env)
-# assign("deltaX", matrix(0, n-1, d.size), env=env)
-# for(t in 1:(n-1))
-# env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
-
- assign("X", as.matrix(onezoo(yuima)[1:k,]), env=env)
- assign("deltaX", matrix(0, k-1, d.size), env=env)
- for(t in 1:(k-1))
- env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
-
-
- assign("h", deltat(yuima at data@zoo.data[[1]]), env=env)
-
- f <- function(p) {
- mycoef <- as.list(p)
- names(mycoef) <- nm
- sum(pminusquasilogl(yuima=yuima, param=mycoef, print=print, env))
- }
-
-
- oout <- NULL
-
-
- if(length(start)>1){ # multidimensional optim
- oout <- optim(start, f, method = method, hessian = FALSE, lower=lower, upper=upper)
-
- } else { ### one dimensional optim
- opt <- optimize(f, ...) ## an interval should be provided
- oout <- list(par = opt$minimum, value = opt$objective)
- } ### endif( length(start)>1 )
-
-
- oout <- oout[c("par","value")]
- return(oout)
+ if(t<minT || t>maxT)
+ syuima.stop("time 't' out of bounds")
+ grid <- times[which(times<=t)]
+ qmle(subsampling(yuima, grid=grid), ...)
}
-
-qmleR <- function(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...){
+qmleR <- function(yuima, t, ...){
- 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
-
- if(!is.list(start))
- yuima.stop("Argument 'start' must be of list type.")
-
- fullcoef <- diff.par
- npar <- length(fullcoef)
-
- 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)
-
- 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]
-
times <- time(yuima at data@zoo.data[[1]])
minT <- as.numeric(times[1])
maxT <- as.numeric(times[length(times)])
- if(missing(t))
+ if(missing(t) )
t <- mean(c(minT,maxT))
- k <- max(which(times <=t),na.rm=TRUE)[1]
- if(k<2)
- k <- 2
- if(k>n-2)
- k <- n-2
+ if(t<minT || t>maxT)
+ syuima.stop("time 't' out of bounds")
+ grid <- times[which(times>=t)]
+ qmle(subsampling(yuima, grid=grid), ...)
+}
- env <- new.env()
- assign("X", as.matrix(onezoo(yuima)[-(1:k),]), env=env)
- assign("deltaX", matrix(0, dim(env$X)[1]-1, d.size), env=env)
- for(t in 1:(dim(env$X)[1]-1))
- env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
-
- assign("h", deltat(yuima at data@zoo.data[[1]]), env=env)
-
- f <- function(p) {
- mycoef <- as.list(p)
- names(mycoef) <- nm
- sum(pminusquasilogl(yuima=yuima, param=mycoef, print=print, env))
- }
-
-
- oout <- NULL
-
-
- if(length(start)>1){ # multidimensional optim
- oout <- optim(start, f, method = method, hessian = FALSE, lower=lower, upper=upper)
-
- } else { ### one dimensional optim
- opt <- optimize(f, ...) ## an interval should be provided
- oout <- list(par = opt$minimum, value = opt$objective)
- } ### endif( length(start)>1 )
-
-
- oout <- oout[c("par","value")]
- return(oout)
-}
CPoint <- function(yuima, param1, param2, print=FALSE, plot=FALSE){
d.size <- yuima at model@equation.number
Modified: pkg/yuima/R/subsampling.R
===================================================================
--- pkg/yuima/R/subsampling.R 2010-07-20 12:29:29 UTC (rev 113)
+++ pkg/yuima/R/subsampling.R 2010-07-20 14:54:54 UTC (rev 114)
@@ -35,11 +35,12 @@
tmpsamp <- sampling
#}
-
+
Data <- get.zoo.data(x)
n.data <- length(Data)
tmpgrid <- vector(n.data, mode="list")
-
+ tmpsamp at grid <- rep(tmpsamp at grid, n.data)[1:n.data]
+
# prepares a grid of times
if(is.logical(tmpsamp at random)){
@@ -48,7 +49,7 @@
if(length(tmpsamp at delta) < n.data)
tmpsamp at delta <- rep( tmpsamp at delta, n.data)[1:n.data]
for(i in 1:n.data){
- tmpgrid[[i]] <- seq(start(Data[[i]]), end(Data[[i]]), by=tmpsamp at delta[i])
+ tmpgrid[[i]] <- tmpsamp at grid[[i]] #seq(start(Data[[i]]), end(Data[[i]]), by=tmpsamp at delta[i])
tmpsamp at regular[i] <- TRUE
tmpsamp at random[i] <- FALSE
}
Modified: pkg/yuima/R/yuima.sampling.R
===================================================================
--- pkg/yuima/R/yuima.sampling.R 2010-07-20 12:29:29 UTC (rev 113)
+++ pkg/yuima/R/yuima.sampling.R 2010-07-20 14:54:54 UTC (rev 114)
@@ -20,7 +20,7 @@
# which.delta: check is grid is regular. If regular returns the delta, otherwise NA
-which.delta <- function(x) ifelse(length(unique(diff(x)))==1, diff(x)[1], NA)
+which.delta <- function(x) ifelse(length(unique(round(diff(x),7)))==1, diff(x)[1], NA)
setMethod("initialize", "yuima.sampling",
function(.Object, Initial, Terminal, n, delta, grid, random,
@@ -42,7 +42,7 @@
.Object at Terminal <- sapply(grid, max)
.Object at n <- sapply(grid, function(x) length(x))
.Object at random <- FALSE
- .Object at delta <- sapply(grid, which.delta)
+ .Object at delta <- as.numeric(sapply(grid, which.delta))
.Object at regular <- !any(is.na( .Object at delta ) )
return(.Object)
}
Modified: pkg/yuima/man/CPoint.Rd
===================================================================
--- pkg/yuima/man/CPoint.Rd 2010-07-20 12:29:29 UTC (rev 113)
+++ pkg/yuima/man/CPoint.Rd 2010-07-20 14:54:54 UTC (rev 114)
@@ -6,21 +6,14 @@
\description{Adaptive LASSO estimation for stochastic differential equations.}
\usage{
CPoint(yuima, param1, param2, print=FALSE, plot=FALSE)
-qmleL(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...)
-qmleR(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...)
+qmleL(yuima, t, ...)
+qmleR(yuima, t, ...)
}
\arguments{
\item{yuima}{a yuima object.}
\item{param1}{parameter values before the change point t}
\item{param2}{parameter values after the change point t}
- \item{print}{produces some output.}
- \item{t}{right end point for \code{qmleL}, left end point for \code{qmleR}. See Details.}
- \item{plot}{plot change point test statistics? Default is \code{FALSE}.}
- \item{...}{passed to \code{\link{optim}} method. See Examples.}
- \item{lower}{a named list for specifying lower bounds of parameters}
- \item{upper}{a named list for specifying upper bounds of parameters}
- \item{start}{initial values to be passed to the optimizer.}
- \item{method}{see Details.}
+ \item{...}{passed to \code{\link{qmle}} method. See Examples.}
}
\details{
@@ -94,20 +87,93 @@
# first stage changepoint estimator
-t.est2 <- CPoint(yuima,param1=tmp1$par,param2=tmp2$par)
+t.est2 <- CPoint(yuima,param1=coef(tmp1),param2=coef(tmp2))
t.est2$tau
# second stage estimation of parameters given first stage
# change point estimator
-tmp11 <- qmleL(yuima,start=as.list(tmp1$par), t=t.est2$tau-0.1,lower=low, method="L-BFGS-B")
+tmp11 <- qmleL(yuima,start=as.list(coef(tmp1)), t=t.est2$tau-0.1,lower=low, method="L-BFGS-B")
tmp11
-tmp21 <- qmleR(yuima,start=as.list(tmp2$par), t=t.est2$tau+0.1,lower=low, method="L-BFGS-B")
+tmp21 <- qmleR(yuima,start=as.list(coef(tmp2)), t=t.est2$tau+0.1,lower=low, method="L-BFGS-B")
tmp21
# second stage estimator of the change point
-CPoint(yuima,param1=tmp11$par,param2=tmp21$par)
+CPoint(yuima,param1=coef(tmp11),param2=coef(tmp21))
+
+## One dimensional example
+diff.matrix <- matrix("(1+x1^2)^theta1", 1, 1)
+drift.c <- c("x1")
+
+ymodel <- setModel(drift=drift.c, diffusion=diff.matrix, time.variable="t",
+state.variable=c("x1"), solve.variable=c("x1"))
+n <- 500
+
+set.seed(123)
+
+y0 <- 5 # initial value
+theta00 <- 1/5
+gamma <- 1/4
+
+
+theta01 <- theta0+n^(-gamma)
+
+
+t1 <- list(theta1= theta00)
+t2 <- list(theta1= theta01)
+
+tau <- 0.4
+ysamp1 <- setSampling(n=tau*n, Initial=0, delta=1/n)
+yuima1 <- setYuima(model=ymodel, sampling=ysamp1)
+yuima1 <- simulate(yuima1, xinit=c(5), true.parameter=t1)
+x1 <- yuima1 at data@zoo.data[[1]]
+x1 <- as.numeric(x1[length(x1)])
+
+ysamp2 <- setSampling(Initial=tau, n=n*(1-tau), delta=1/n)
+yuima2 <- setYuima(model=ymodel, sampling=ysamp2)
+
+yuima2 <- simulate(yuima2, xinit=c(x1), true.parameter=t2)
+
+
+yuima <- yuima1
+yuima at data@zoo.data[[1]] <- c(yuima1 at data@zoo.data[[1]], yuima2 at data@zoo.data[[1]][-1])
+
+
+plot(yuima)
+
+
+t.est <- CPoint(yuima,param1=t1,param2=t2)
+t.est$tau
+
+low <- list(theta1=0)
+upp <- list(theta1=1)
+
+# first state estimate of parameters using small
+# portion of data in the tails
+tmp1 <- qmleL(yuima,start=list(theta1=0.5),t=.15,lower=low, upper=upp,method="L-BFGS-B")
+tmp1
+tmp2 <- qmleR(yuima,start=list(theta1=0.5), t=.85,lower=low, upper=upp,method="L-BFGS-B")
+tmp2
+
+
+
+# first stage changepoint estimator
+t.est2 <- CPoint(yuima,param1=coef(tmp1),param2=coef(tmp2))
+t.est2$tau
+
+
+# second stage estimation of parameters given first stage
+# change point estimator
+tmp11 <- qmleL(yuima,start=as.list(coef(tmp1)), t=t.est2$tau-0.1, lower=low, upper=upp,method="L-BFGS-B")
+tmp11
+
+tmp21 <- qmleR(yuima,start=as.list(coef(tmp2)), t=t.est2$tau+0.1,lower=low, upper=upp,method="L-BFGS-B")
+tmp21
+
+
+# second stage estimator of the change point
+CPoint(yuima,param1=coef(tmp11),param2=coef(tmp21),plot=TRUE)
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
More information about the Yuima-commits
mailing list