[Yuima-commits] r6 - in pkg/yuima: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 27 10:25:24 CET 2009
Author: hinohide
Date: 2009-10-27 10:25:19 +0100 (Tue, 27 Oct 2009)
New Revision: 6
Modified:
pkg/yuima/R/simulate.R
pkg/yuima/R/yuima.model.R
pkg/yuima/man/simulate.Rd
Log:
simulate.R{d} and yuima.model.R have changed
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2009-10-27 08:50:46 UTC (rev 5)
+++ pkg/yuima/R/simulate.R 2009-10-27 09:25:19 UTC (rev 6)
@@ -42,7 +42,7 @@
## readline()
if(missing(true.parameter)){
- true.parameter <- as.list(numeric(length(sdeModel at parameter@all)))
+ true.parameter <- numeric(length(sdeModel at parameter@all))
}
r.size <- sdeModel at noise.number
@@ -109,15 +109,15 @@
modeltime <- sdeModel at time.variable
V0 <- sdeModel at drift
V <- sdeModel at diffusion
-
+
par.len <- length(sdeModel at parameter@all)
if(par.len>0){
for(i in 1:par.len){
pars <- sdeModel at parameter@all[i]
- assign(pars, true.parameter[[i]])
+ assign(pars, true.parameter[i])
}
}
-
+
##:: Initialization
##:: set time step
delta <- Terminal/division
@@ -351,28 +351,24 @@
}else if(sdeModel at measure.type=="code"){ ##:: code type
##:: Jump terms
code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
- args <- eval(parse(text=paste("list(",suppressWarnings(sub("^.+?\\(.+?,(.+)", "\\1", sdeModel at measure$df$expr, perl=TRUE)))))
-
- # args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
- #print(args);readline()
+ args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
+ print(args);readline()
dZ <- switch(code,
- #rNIG=paste("rNIG(division, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta)"),
- rNIG=rNIG(division, args[[1]], args[[2]], args[[3]]*delta, args[[4]]*delta,args[[5]]*delta),
- rIG=rIG(division,args[[1]]*delta, args[[2]]),
- rgamma=rgamma(division,args[[1]]*delta, args[[2]]),
- rbgamma=rbgamma(division, args[[1]]*delta, args[[2]], args[[3]]*delta, args[[4]]),
+ rNIG=paste("rNIG(division, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta)"),
+ rIG=paste("rIG(division,", args[2], "*delta, ", args[3], ")"),
+ rgamma=paste("rgamma(division, ", args[2], "*delta, ", args[3], ")"),
+ rbgamma=paste("rbgamma(division, ", args[2], "*delta, ", args[3], ", ", args[4], "*delta, ", args[5], ")"),
## rngamma=paste("rngamma(division, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta, ", args[6], ")"),
- rngamma=rngamma(division, args[[1]]*delta, args[[2]], args[[3]], args[[4]]*delta),
+ rngamma=paste("rngamma(division, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta)"),
## rstable=paste("rstable(division, ", args[2], ", ", args[3], ", ", args[4], ", ", args[5], ", ", args[6], ")")
- rstable=rstable(division, args[[1]], args[[2]], args[[3]]*delta^(1/args[[1]]), args[[4]]*delta)
+ rstable=paste("rstable(division, ", args[2], ", ", args[3], ", ", args[4], "*delta^(1/",args[2],"), ", args[5], "*delta)")
)
if(is.null(dZ)){ ##:: "otherwise"
cat(paste("Code \"", code, "\" not supported yet.\n", sep=""))
return(NULL)
}
- #print(parse(text=dZ))
- #dZ <- eval(parse(text=dZ))
+ dZ <- eval(parse(text=dZ))
##:: calcurate difference eq.
for(i in 1:division){
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2009-10-27 08:50:46 UTC (rev 5)
+++ pkg/yuima/R/yuima.model.R 2009-10-27 09:25:19 UTC (rev 6)
@@ -83,10 +83,10 @@
if( !length(jump.coeff) || !length(measure) ){
cat("\nmeasure type isn't matched with jump term.\n")
return(NULL)
- }#else if(length(jump.coeff)!=1){
- # cat("\nmulti dimentional jump term is not supported yet.\n")
- # return(NULL)
- #}
+ }else if(length(jump.coeff)!=1){
+ cat("\nmulti dimentional jump term is not supported yet.\n")
+ return(NULL)
+ }
if(measure.type=="CP"){ ##::CP
if(length(measure)!=2){
Modified: pkg/yuima/man/simulate.Rd
===================================================================
--- pkg/yuima/man/simulate.Rd 2009-10-27 08:50:46 UTC (rev 5)
+++ pkg/yuima/man/simulate.Rd 2009-10-27 09:25:19 UTC (rev 6)
@@ -8,7 +8,7 @@
\arguments{
\item{yuima}{an \code{yuima} object.}
\item{xinit}{initial value vector of state variables.}
- \item{true.parameter}{LIST of initial values of parameters.}
+ \item{true.parameter}{initial value vector of parameters.}
\item{space.discretized}{flag to switch to space-discretized Euler
Maruyama method.}
\item{increment.W}{to specify Wiener increment for each time tics in advance.}
@@ -43,7 +43,6 @@
# Solve SDEs using Euler-Maruyama method.
ou <- simulate(yuima=ou, xinit=1)
-x11()
plot(ou)
# A multi-dimensional (correlated) diffusion process.
@@ -123,6 +122,7 @@
space.discretized=FALSE,
increment.W=my.dW)
if( !is.null(yuima.mod) ){
+ x11()
plot(yuima.mod)
}
@@ -166,6 +166,7 @@
## Solve SDEs using Euler-Maruyama method.
yuima.obj.path <- simulate(yuima=yuima.obj, space.discretized=FALSE, increment.W=my.dW)
if( !is.null(yuima.obj.path) ){
+ x11()
plot(yuima.obj.path)
}
@@ -188,7 +189,7 @@
obj.sampling <- setSampling(Terminal=T, division=division)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
-X <- simulate(yuima=obj.yuima, xinit=xinit, true.parameter=list(theta, sigma))
+X <- simulate(yuima=obj.yuima, xinit=xinit, true.parameter=c(theta, sigma))
plot(X)
@@ -204,14 +205,15 @@
## dX = (theta - A X)dt + dZ,
## theta=(theta_1, theta_2) = c(1,.5)
## A=[a_ij], a_11 = 2, a_12 = 1, a_21 = 1, a_22=2
- xinit <- c(1,1)
- beta <- c(.1,.1)
- Lambda <- matrix(c(1,0,0,1),2,2)
- obj.model <- setModel(drift=c("1 - 2*x1-x2",".5-x1-2*x2"), xinit=c(1,1), solve.variable=c("x1","x2"), jump.coeff=c("1","1"), measure.type="code", measure=list(df="rNIG(z, alpha=1, beta=beta, delta=1,mu=c(0,0), Lambda=Lambda)"))
- obj.sampling <- setSampling(Terminal=10, division=10000)
- obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
- result <- simulate(yuima=obj.yuima,true.param=list(0,0,beta,Lambda))
- plot(result)
+## xinit <- c(1,1)
+## beta <- c(.1,.1)
+## mu <- c(0,0)
+## Lambda <- matrix(c(1,0,0,1),2,2)
+## obj.model <- setModel(drift=c("1 - 2*x1-x2",".5-x1-2*x2"), xinit=c(1,1), solve.variable=c("x1","x2"), jump.coeff="1", measure.type="code", measure=list(df="rIG(z, alpha=1, beta=beta, mu=mu, Lambda=Lambda)"))
+## obj.sampling <- setSampling(Terminal=10, division=10000)
+## obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
+## result <- simulate(yuima=obj.yuima)
+## plot(result)
}
More information about the Yuima-commits
mailing list