[Yuima-commits] r70 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 3 02:18:46 CEST 2010
Author: hinohide
Date: 2010-05-03 02:18:46 +0200 (Mon, 03 May 2010)
New Revision: 70
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/AllClasses.R
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/R/yuima.data.R
pkg/yuima/man/adaBayes.Rd
pkg/yuima/man/quasi-likelihood.Rd
pkg/yuima/man/setData.Rd
pkg/yuima/man/yuima-class.Rd
pkg/yuima/man/yuima.data-class.Rd
Log:
changed return object of ml.ql to mle, and add LSE function for initial guess for ml.ql
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/DESCRIPTION 2010-05-03 00:18:46 UTC (rev 70)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 0.0.83
-Date: 2010-03-20
+Version: 0.0.84
+Date: 2010-05-04
Depends: methods, zoo, adapt, stats4
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/NAMESPACE 2010-05-03 00:18:46 UTC (rev 70)
@@ -24,7 +24,6 @@
"cce",
"poisson.random.sampling",
"get.zoo.data",
- "get.mle",
"initialize",
"ql",
"rql",
@@ -40,10 +39,7 @@
"simFunctional",
"F0",
"Fnorm",
- "asymptotic_term",
-
- "get.mle"
-
+ "asymptotic_term"
)
## function which we want to expose to the user
@@ -68,7 +64,6 @@
export(poisson.random.sampling)
export(get.zoo.data)
-export(get.mle)
export(ql,rql,ml.ql)
export(adaBayes)
@@ -86,4 +81,6 @@
export(Fnorm)
export(asymptotic_term)
+export(LSE)
+
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/R/AllClasses.R 2010-05-03 00:18:46 UTC (rev 70)
@@ -42,8 +42,7 @@
# classes
setClass("yuima.data", representation(original.data = "ANY",
- zoo.data = "ANY",
- mle="ANY"
+ zoo.data = "ANY"
)
)
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/R/quasi-likelihood.R 2010-05-03 00:18:46 UTC (rev 70)
@@ -955,10 +955,8 @@
vcov=vcov, min=min, details=opt, minuslogl=ql.opt,
method=method)
##END convert to mle
-
- yuima at data@mle <- opt
- return(yuima)
+ return(opt)
})
@@ -992,3 +990,83 @@
return(ret)
})
+
+##estimate theta2 by LSE
+##function name LSE
+setGeneric("LSE", function(yuima, h, theta2.init, interval, ...)
+ standardGeneric("LSE")
+ )
+setMethod("LSE", "yuima",
+ function(yuima, h, theta2.init=c(), interval=c(0,1), ...){
+ ##theta2.init : multi-dim param used by optim()
+ ##interval : 1-dim param used by optimize()
+
+ ##objective function
+ Mn <-function(theta2){
+ Mn.part <- function(deltaX, h, yuima){
+ tmp <- deltaX - h %*% eval(yuima at model@drift)
+ ret <- t(tmp) %*% tmp
+ return(ret)
+ }
+
+ ##init
+ sum.tmp <- 0
+ X <- NULL
+ X.tmp <- get.zoo.data(yuima)
+ for(i in 1:length(X.tmp)){
+ X <- cbind(X, as.matrix(X.tmp[[i]]))
+ }
+ modelpara.drift <- yuima at model@parameter at drift
+ modelstate <- yuima at model@state.variable
+ for(i in 1:length(theta2)){
+ assign(modelpara.drift[i], theta2[i])
+ }
+
+ ##sum loop
+ for(j in 2:dim(X)[1]){
+ ##get param
+ x <- X[j-1,]
+ deltaX <- X[j,] - X[j-1,]
+ for(k in 1:length(x)){
+ assign(modelstate[k], x[k])
+ }
+ ##calc
+ sum.tmp <- sum.tmp + Mn.part(deltaX, h, yuima)
+ }
+
+ #cat("theta2 value:", theta2, "\n")
+ #cat("Mn value:", sum.tmp, "\n\n")
+
+ return(sum.tmp)
+
+ }
+ ##END objective function
+
+ if(length(yuima at model@parameter at drift)==1){
+
+ if(missing(interval)){
+ stop("\ninterval missing.\n")
+ }
+ if(!missing(theta2.init)){
+ cat("\ntheta2.init is ignored.\n")
+ }
+
+ opt <- optimize(f=Mn, interval=interval, tol=1e-100, ...)
+ opt <- list(par=opt$minimum, value=opt$objective)
+
+ }else{
+
+ if(missing(theta2.init)){
+ stop("\ntheta2.init missing.\n")
+ }
+ if(!missing(interval)){
+ cat("\ninterval is ignored.\n")
+ }
+
+ opt <- optim(c(theta2.init), fn=Mn, gr=NULL, ...)
+ opt <- list(par=opt$par, value=opt$value)
+
+ }
+ return(opt)
+ })
+
Modified: pkg/yuima/R/yuima.data.R
===================================================================
--- pkg/yuima/R/yuima.data.R 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/R/yuima.data.R 2010-05-03 00:18:46 UTC (rev 70)
@@ -43,16 +43,6 @@
return(x at zoo.data)
})
-setGeneric("get.mle",
- function(x)
- standardGeneric("get.mle")
- )
-
-setMethod("get.mle", signature(x="yuima.data"),
- function(x){
- return(x at mle)
- })
-
# following funcs are basic generic funcs
setGeneric("plot",
@@ -120,12 +110,6 @@
# same methods for 'yuima'. Depend on same methods for 'data'
-
-setMethod("get.mle", "yuima",
- function(x){
- return(get.mle(x at data))
- })
-
setMethod("get.zoo.data", "yuima",
function(x){
return(get.zoo.data(x at data))
Modified: pkg/yuima/man/adaBayes.Rd
===================================================================
--- pkg/yuima/man/adaBayes.Rd 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/man/adaBayes.Rd 2010-05-03 00:18:46 UTC (rev 70)
@@ -57,7 +57,6 @@
param <- numeric(2)
for( i in 1:5){
PAR <- ml.ql(yuima.mod,theta2=rand[i,1],theta1=rand[i,2],h=(n^(-2/3)))
- PAR <- get.mle(PAR)
if(PAR at min > QL){
QL <- PAR at min
param <- PAR at coef
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-05-03 00:18:46 UTC (rev 70)
@@ -1,6 +1,7 @@
\name{ql}
\alias{ql}
\alias{ml.ql}
+\alias{LSE}
%\alias{ml.ql2}
\alias{rql}
%\alias{ql,ANY-method}
@@ -50,7 +51,7 @@
ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
yuima <- setYuima(model=ymodel, sampling=ysamp)
set.seed(123)
-yuima <- simulate(yuima, xinit=1, true.parameter=list(theta2=0.3,theta1=0.1))
+yuima <- simulate(yuima, xinit=1, true.parameter=c(0.3, 0.1))
QL <- ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)))
QL
@@ -60,30 +61,30 @@
##QL
system.time(
-yuima <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
+opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
)
print("True param")
print("theta2 = .3, theta1 = .1")
print("ML estimator")
-get.mle(yuima)@coef
+opt at coef
## another way of parameter specification
##interval <- list(theta2.lim=c(0,1), theta1.lim=c(0,1))
##system.time(
-##yuima <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
+##opt <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
##)
##print("True param")
##print("theta2 = .3, theta1 = .1")
##print("ML estimator")
-#get.mle(yuima)@coef
+#opt at coef
system.time(
-yuima <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newtoon")
+opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newtoon")
)
print("True param")
print("theta2 = .3, theta1 = .1")
print("ML estimator")
-get.mle(yuima)@coef
+opt at coef
##multidimension case
##dXt^e = - drift.matrix * Xt^e * dt + diff.matrix * dWt
@@ -100,7 +101,7 @@
set.seed(123)
##xinit=c(x1, x2) #true.parameter=c(theta2.1, theta2.2, theta1.1, theta1.2)
-yuima <- simulate(yuima, xinit=c(1, 1), true.parameter=list(theta2.1=0.5, theta2.2=0.3, theta1.1=0.6,theta1.2=0.2))
+yuima <- simulate(yuima, xinit=c(1, 1), true.parameter=c(0.5, 0.3, 0.6, 0.2))
theta2 <- c(0.8, 0.2) #c(theta2.1, theta2.2)
theta1 <- c(0.7, 0.1) #c(theta1.1, theta1.2)
@@ -120,57 +121,23 @@
theta1.lim <- t( matrix( c(theta1.1.lim, theta1.2.lim), 2, 2) )
system.time(
-yuima <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim)
+opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim)
)
-get.mle(yuima)@coef
+opt at coef
## another way of parameter specification
#interval <- list(theta2.1.lim, theta2.2.lim, theta1.1.lim, theta1.2.lim)
#system.time(
-#yuima <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
+#opt <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
#)
-#get.mle(yuima)@coef
+#opt at coef
-
-
-##dXt^e = -theta2 * Xt^e * dt + theta1 * dWt
-diff.matrix <- matrix(c("theta1"), 1, 1)
-ymodel <- setModel(drift=c("(-1)*theta2*x"), diffusion=diff.matrix,
- time.variable="t", state.variable="x", solve.variable="x")
-n <- 100
-
-ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
-yuima <- setYuima(model=ymodel, sampling=ysamp)
-set.seed(123)
-yuima <- simulate(yuima, xinit=1, true.parameter=list(theta2=0.5, theta1=0.8))
-#0.3,0.1
-h <- 1/((n)^(2/3))
-
-onezoo <- function(ydata) {
- dat <- get.zoo.data(ydata)
- dats <- dat[[1]]
- if(length(dat)>1) {
- for(i in 2:(length(dat))) {
- dats <- merge(dats,dat[[i]])
- }
- }
- return(dats)
+system.time(
+opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim, method="Newton")
+)
+opt at coef
}
-theta2 <- 0.1#0.65
-theta1 <- 0.2#0.7
-print("use Newton method")
-res1 <- ml.ql(yuima, theta2, theta1, h, method="Newton", print=TRUE)
-print(confint(res1))
-
-
-print("use optim")
-res2 <- ml.ql(yuima, theta2, theta1, h, print=FALSE)
-print(confint(res2))
-
-
-}
-
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ts}
Modified: pkg/yuima/man/setData.Rd
===================================================================
--- pkg/yuima/man/setData.Rd 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/man/setData.Rd 2010-05-03 00:18:46 UTC (rev 70)
@@ -3,7 +3,6 @@
\alias{get.zoo.data}
\alias{dim}
\alias{length}
-\alias{get.mle}
\title{
Set and access data of an object of type "yuima.data" or "yuima".
@@ -15,8 +14,6 @@
\code{\link{yuima.data-class}} object. (Note: value is a \code{list} of
\code{\link{zoo}} objects).
-\code{get.mle} return the content of mle class of ml.ql result.
-
\code{plot} plot method for object of \code{\link{yuima.data-class}} or
\code{\link{yuima-class}}.
@@ -55,11 +52,6 @@
returns the content of the slot \code{zoo.data} of \code{x} if \code{x}
is of \code{\link{yuima.data-class}} or the content of
\code{x at data@zoo.data} if \code{x} is of \code{\link{yuima-class}}.
-
-The function \code{get.mle}
-returns the content of the slot \code{mle} of \code{x} if \code{x}
-is of \code{\link{yuima.data-class}} or the content of
-\code{x at data@mle} if \code{x} is of \code{\link{yuima-class}}.
}
\value{
\item{value}{a list of object(s) of \code{\link{yuima.data-class}} for
Modified: pkg/yuima/man/yuima-class.Rd
===================================================================
--- pkg/yuima/man/yuima-class.Rd 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/man/yuima-class.Rd 2010-05-03 00:18:46 UTC (rev 70)
@@ -1,8 +1,6 @@
\name{yuima-class}
\docType{class}
\alias{yuima-class}
-
-\alias{get.mle,yuima-method}
\alias{get.zoo.data,yuima-method}
\alias{plot,yuima,ANY-method}
\alias{dim,yuima-method}
@@ -16,6 +14,7 @@
\alias{rql,yuima-method}
\alias{ml.ql,yuima-method}
\alias{limiting.gamma,yuima-method}
+\alias{LSE,yuima-method}
\title{Class for stochastic differential equations}
\description{
@@ -79,8 +78,6 @@
For more details see \code{\link{cce}}.}
\item{simulate}{simulation method. For more information see
\code{\link{simulate}}.}
- \item{get.mle}{\code{signature(x="yuima")}: return mle class object
- of mle result.}
}
}
\author{The YUIMA Project Team}
Modified: pkg/yuima/man/yuima.data-class.Rd
===================================================================
--- pkg/yuima/man/yuima.data-class.Rd 2010-03-27 04:15:24 UTC (rev 69)
+++ pkg/yuima/man/yuima.data-class.Rd 2010-05-03 00:18:46 UTC (rev 70)
@@ -14,7 +14,6 @@
%\alias{setSampling,yuima.data-method}
\alias{poisson.random.sampling,yuima.data-method}
\alias{subsampling,yuima.data-method}
-\alias{get.mle,yuima.data-method}
\title{Class "yuima.data" for the data slot of a "yuima" class object}
\description{
@@ -28,7 +27,6 @@
\describe{
\item{\code{original.data}:}{The original data.}
\item{\code{zoo.data}:}{A list of \code{zoo} format data.}
- \item{\code{mle}:}{mle class object, result of ml.ql}
}
}
\details{
@@ -87,8 +85,6 @@
\item{cce}{\code{signature(x = "yuima.data")}: calculates asyncronous
covariance estimator on the data contained in \code{x at zoo.data}.
For more details see \code{\link{cce}}.}
- \item{get.mle}{\code{signature(x="yuima.data")}: return mle class of
- ml.ql result.}
}
}
\author{The YUIMA Project Team}
More information about the Yuima-commits
mailing list