[Yuima-commits] r290 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 19 17:33:42 CET 2014
Author: iacus
Date: 2014-03-19 17:33:41 +0100 (Wed, 19 Mar 2014)
New Revision: 290
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/asymptotic_term_third_function.R
pkg/yuima/R/yuima.R
pkg/yuima/R/yuima.model.R
Log:
up
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-03-19 14:02:57 UTC (rev 289)
+++ pkg/yuima/DESCRIPTION 2014-03-19 16:33:41 UTC (rev 290)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 1.0.0
+Version: 1.0.1
Date: 2014-03-19
Depends: methods, zoo, stats4, utils, expm
Suggests: cubature, mvtnorm
Modified: pkg/yuima/R/asymptotic_term_third_function.R
===================================================================
--- pkg/yuima/R/asymptotic_term_third_function.R 2014-03-19 14:02:57 UTC (rev 289)
+++ pkg/yuima/R/asymptotic_term_third_function.R 2014-03-19 16:33:41 UTC (rev 290)
@@ -9019,179 +9019,180 @@
}
- p2 <- function(z,get_F_tilde1__2,get_F_tilde2,
- get_F_tilde1H1,get_H2_1,get_H2_2,get_H2_3,env){
+# p2 <- function(z,get_F_tilde1__2,get_F_tilde2,
+# get_F_tilde1H1,get_H2_1,get_H2_2,get_H2_3,env){
+#
+# k.size <- env$k.size
+# delta <- env$delta
+# mu <- env$mu
+# Sigma <- env$Sigma
+#
+# z.tilde <- z - mu
+#
+# first <- 0
+# second <- 0
+# third <- 0 #added
+#
+# if(k.size == 1){
+#
+# first <- (F_tilde1__2_x(1,1,z.tilde+2*delta,get_F_tilde1__2,env) *
+# dnorm(z+2*delta,mean=mu,sd=sqrt(Sigma)) -
+# 2 * F_tilde1__2_x(1,1,z.tilde+delta,get_F_tilde1__2,env) *
+# dnorm(z+delta,mean=mu,sd=sqrt(Sigma)) +
+# F_tilde1__2_x(1,1,z.tilde,get_F_tilde1__2,env) *
+# dnorm(z,mean=mu,sd=sqrt(Sigma)))/(delta)^2
+#
+# first <- first/2
+#
+# second <- (F_tilde2_x(1,z.tilde+delta,get_F_tilde2,env) *
+# dnorm(z+delta,mean=mu,sd=sqrt(Sigma)) -
+# F_tilde2_x(1,z.tilde,get_F_tilde2,env) *
+# dnorm(z,mean=mu,sd=sqrt(Sigma)))/delta
+#
+##added:start
+# obj1 <- F_tilde1H1_x(1,z.tilde+delta,get_F_tilde1H1,env)
+#
+# obj2 <- F_tilde1H1_x(1,z.tilde,get_F_tilde1H1,env)
+#
+# third <- (obj1 * dnorm(z+delta,mean=mu,sd=sqrt(Sigma)) -
+# obj2 * dnorm(z,mean=mu,sd=sqrt(Sigma)))/delta
+#
+# forth <- H2_x(z.tilde,get_H2_1,get_H2_2,get_H2_3,env) * #dnorm(z,mean=mu,sd=sqrt(Sigma))
+#adde#d:end
+#
+# }else{
+#
+# tmp1 <- matrix(0,k.size,k.size)
+#
+# for(l1 in 1:k.size){
+# for(l2 in 1:k.size){
+#
+# dif1 <- numeric(k.size)
+# dif2 <- numeric(k.size)
+# dif1[l1] <- dif1[l1] + delta
+# dif2[l2] <- dif2[l2] + delta
+#
+# tmp1[l1,l2] <- (F_tilde1__2_x(l1,l2,z.tilde+dif1+dif2,get_F_tilde1__2,env) *
+# ft_norm(z+dif1+dif2,env) -
+# F_tilde1__2_x(l1,l2,z.tilde+dif1,get_F_tilde1__2,env) *
+# ft_norm(z+dif1,env) -
+# F_tilde1__2_x(l1,l2,z.tilde+dif2,get_F_tilde1__2,env) *
+# ft_norm(z+dif2,env) +
+# F_tilde1__2_x(l1,l2,z.tilde,get_F_tilde1__2,env) *
+# ft_norm(z,env))/(delta)^2
+#
+# }
+# }
+#
+# first <- sum(tmp1)/2
+#
+# tmp2 <- double(k.size)
+# tmp3 <- double(k.size) #added
+#
+# for(l in 1:k.size){
+#
+# dif <- numeric(k.size)
+# dif[l] <- dif[l] + delta
+#
+# tmp2[l] <- (F_tilde2_x(l,z.tilde+dif,get_F_tilde2,env) *
+# ft_norm(z+dif,env) -
+# F_tilde2_x(l,z.tilde,get_F_tilde2,env) *
+# ft_norm(z,env))/delta
+#
+##added:start
+# obj1 <- F_tilde1H1_x(l,z.tilde+dif,get_F_tilde1H1,env)
+#
+# obj2 <- F_tilde1H1_x(l,z.tilde,get_F_tilde1H1,env)
+#
+# tmp3[l] <- (obj1 * ft_norm(z+dif,env) -
+# obj2 * ft_norm(z,env))/delta
+##added:end
+#
+# }
+#
+# second <- sum(tmp2)
+# third <- sum(tmp3) #added
+# forth <- H2_x(z.tilde,get_H2_1,get_H2_2,get_H2_3,env)*ft_norm(z,env) #added
+# }
+#
+# result <- first - second -
+# third + forth #added
+# return(result)
+# }
+#
+#
+# p2_z <- function(z){
+#
+# if(k.size == 1){
+# zlen <- length(z)
+# result <- c()
+#
+# for(i in 1:zlen){
+# result[i] <- p2(z[i],get_F_tilde1__2,get_F_tilde2,
+# get_F_tilde1H1,get_H2_1,get_H2_2,get_H2_3,env)
+# }
+# }else{
+# result <- p2(z,get_F_tilde1__2,get_F_tilde2,
+# get_F_tilde1H1,get_H2_1,get_H2_2,get_H2_3,env)
+# }
+#
+# return(result)
+# }
- k.size <- env$k.size
- delta <- env$delta
- mu <- env$mu
- Sigma <- env$Sigma
- z.tilde <- z - mu
+# d2.term <- function(){
+#
+# gz_p2 <- function(z){
+#
+# result <- H0 * G(z) * p2_z(z)
+# return(result)
+# }
+#
+# if(k.size == 1){
+#
+# ztmp <- seq(mu-7*sqrt(Sigma),mu+7*sqrt(Sigma),length=1000)
+# dt <- ztmp[2] - ztmp[1]
+#
+# p2tmp <- gz_p2(ztmp)
+#
+# result <- sum(p2tmp) * dt
+#
+# }else if(2 <= k.size || k.size <= 20){
+#
+# lambda <- eigen(Sigma)$values
+# matA <- eigen(Sigma)$vector
+#
+# gz_p2 <- function(z){
+# tmpz <- matA %*% z
+# tmp <- H0 * G(tmpz) * p2_z(tmpz) #det(matA) = 1
+# return( tmp )
+# }
+#
+# my.x <- matrix(0,k.size,20^k.size)
+# dt <- 1
+#
+# for(k in 1:k.size){
+# max <- 5 * sqrt(lambda[k])
+# min <- -5 * sqrt(lambda[k])
+# tmp.x <- seq(min,max,length=20)
+# dt <- dt * (tmp.x[2] - tmp.x[1])
+# my.x[k,] <- rep(tmp.x,each=20^(k.size-k),times=20^(k-1))
+# }
+#
+# tmp <- 0
+#
+# for(i in 1:20^k.size){
+# tmp <- tmp + gz_pi1(my.x[,i])
+# }
+#
+# tmp <- tmp * dt
+#
+# }else{
+# stop("length k is too long.")
+# }
+#
+# return(result)
+# }
- first <- 0
- second <- 0
- third <- 0 #added
- if(k.size == 1){
-
- first <- (F_tilde1__2_x(1,1,z.tilde+2*delta,get_F_tilde1__2,env) *
- dnorm(z+2*delta,mean=mu,sd=sqrt(Sigma)) -
- 2 * F_tilde1__2_x(1,1,z.tilde+delta,get_F_tilde1__2,env) *
- dnorm(z+delta,mean=mu,sd=sqrt(Sigma)) +
- F_tilde1__2_x(1,1,z.tilde,get_F_tilde1__2,env) *
- dnorm(z,mean=mu,sd=sqrt(Sigma)))/(delta)^2
-
- first <- first/2
-
- second <- (F_tilde2_x(1,z.tilde+delta,get_F_tilde2,env) *
- dnorm(z+delta,mean=mu,sd=sqrt(Sigma)) -
- F_tilde2_x(1,z.tilde,get_F_tilde2,env) *
- dnorm(z,mean=mu,sd=sqrt(Sigma)))/delta
-
-#added:start
- obj1 <- F_tilde1H1_x(1,z.tilde+delta,get_F_tilde1H1,env)
-
- obj2 <- F_tilde1H1_x(1,z.tilde,get_F_tilde1H1,env)
-
- third <- (obj1 * dnorm(z+delta,mean=mu,sd=sqrt(Sigma)) -
- obj2 * dnorm(z,mean=mu,sd=sqrt(Sigma)))/delta
-
- forth <- H2_x(z.tilde,get_H2_1,get_H2_2,get_H2_3,env) * dnorm(z,mean=mu,sd=sqrt(Sigma))
-#added:end
-
- }else{
-
- tmp1 <- matrix(0,k.size,k.size)
-
- for(l1 in 1:k.size){
- for(l2 in 1:k.size){
-
- dif1 <- numeric(k.size)
- dif2 <- numeric(k.size)
- dif1[l1] <- dif1[l1] + delta
- dif2[l2] <- dif2[l2] + delta
-
- tmp1[l1,l2] <- (F_tilde1__2_x(l1,l2,z.tilde+dif1+dif2,get_F_tilde1__2,env) *
- ft_norm(z+dif1+dif2,env) -
- F_tilde1__2_x(l1,l2,z.tilde+dif1,get_F_tilde1__2,env) *
- ft_norm(z+dif1,env) -
- F_tilde1__2_x(l1,l2,z.tilde+dif2,get_F_tilde1__2,env) *
- ft_norm(z+dif2,env) +
- F_tilde1__2_x(l1,l2,z.tilde,get_F_tilde1__2,env) *
- ft_norm(z,env))/(delta)^2
-
- }
- }
-
- first <- sum(tmp1)/2
-
- tmp2 <- double(k.size)
- tmp3 <- double(k.size) #added
-
- for(l in 1:k.size){
-
- dif <- numeric(k.size)
- dif[l] <- dif[l] + delta
-
- tmp2[l] <- (F_tilde2_x(l,z.tilde+dif,get_F_tilde2,env) *
- ft_norm(z+dif,env) -
- F_tilde2_x(l,z.tilde,get_F_tilde2,env) *
- ft_norm(z,env))/delta
-
-#added:start
- obj1 <- F_tilde1H1_x(l,z.tilde+dif,get_F_tilde1H1,env)
-
- obj2 <- F_tilde1H1_x(l,z.tilde,get_F_tilde1H1,env)
-
- tmp3[l] <- (obj1 * ft_norm(z+dif,env) -
- obj2 * ft_norm(z,env))/delta
-#added:end
-
- }
-
- second <- sum(tmp2)
- third <- sum(tmp3) #added
- forth <- H2_x(z.tilde,get_H2_1,get_H2_2,get_H2_3,env)*ft_norm(z,env) #added
- }
-
- result <- first - second -
- third + forth #added
- return(result)
- }
-
-
- p2_z <- function(z){
-
- if(k.size == 1){
- zlen <- length(z)
- result <- c()
-
- for(i in 1:zlen){
- result[i] <- p2(z[i],get_F_tilde1__2,get_F_tilde2,
- get_F_tilde1H1,get_H2_1,get_H2_2,get_H2_3,env)
- }
- }else{
- result <- p2(z,get_F_tilde1__2,get_F_tilde2,
- get_F_tilde1H1,get_H2_1,get_H2_2,get_H2_3,env)
- }
-
- return(result)
- }
-
-
- d2.term <- function(){
-
- gz_p2 <- function(z){
-
- result <- H0 * G(z) * p2_z(z)
- return(result)
- }
-
- if(k.size == 1){
-
- ztmp <- seq(mu-7*sqrt(Sigma),mu+7*sqrt(Sigma),length=1000)
- dt <- ztmp[2] - ztmp[1]
-
- p2tmp <- gz_p2(ztmp)
-
- result <- sum(p2tmp) * dt
-
- }else if(2 <= k.size || k.size <= 20){
-
- lambda <- eigen(Sigma)$values
- matA <- eigen(Sigma)$vector
-
- gz_p2 <- function(z){
- tmpz <- matA %*% z
- tmp <- H0 * G(tmpz) * p2_z(tmpz) #det(matA) = 1
- return( tmp )
- }
-
- my.x <- matrix(0,k.size,20^k.size)
- dt <- 1
-
- for(k in 1:k.size){
- max <- 5 * sqrt(lambda[k])
- min <- -5 * sqrt(lambda[k])
- tmp.x <- seq(min,max,length=20)
- dt <- dt * (tmp.x[2] - tmp.x[1])
- my.x[k,] <- rep(tmp.x,each=20^(k.size-k),times=20^(k-1))
- }
-
- tmp <- 0
-
- for(i in 1:20^k.size){
- tmp <- tmp + gz_pi1(my.x[,i])
- }
-
- tmp <- tmp * dt
-
- }else{
- stop("length k is too long.")
- }
-
- return(result)
- }
-
-
+#
\ No newline at end of file
Modified: pkg/yuima/R/yuima.R
===================================================================
--- pkg/yuima/R/yuima.R 2014-03-19 14:02:57 UTC (rev 289)
+++ pkg/yuima/R/yuima.R 2014-03-19 16:33:41 UTC (rev 290)
@@ -15,154 +15,168 @@
# This isn't a serious attempt at simplification code. It just does some
# obvious things like 0 + x => x. It was written to support Deriv.R.
-yuima.Simplify<- function(expr)
+yuima.Simplify <- function(expr, yuima.env){
+
+ ###
+
+
+ Simplify_ <- function(expr)
+ {
+ if (is.symbol(expr)) {
+ expr
+ } else if (is.language(expr) && is.symbol(expr[[1]])) {
+ # is there a rule in the table?
+ sym.name <- as.character(expr[[1]])
+ if (class(try(Simplify.rule <-
+ get(sym.name, envir=yuima.env,
+ inherits=FALSE), silent=TRUE))
+ != "try-error")
+ return(Simplify.rule(expr))
+ }
+ expr
+ }
+
+ Simplify.function <- function(f, x=names(formals(f)), env=parent.frame())
+ {
+ stopifnot(is.function(f))
+ as.function(c(as.list(formals(f)),
+ Simplify_(body(f))),
+ envir=env)
+ }
+
+ `Simplify.+` <- function(expr)
+ {
+ if (length(expr) == 2)
+ {
+ if (is.numeric(expr[[2]]))
+ return(+expr[[2]])
+ return(expr)
+ }
+
+ a <- Simplify_(expr[[2]])
+ b <- Simplify_(expr[[3]])
+
+ if (is.numeric(a) && all(a == 0)) {
+ b
+ } else if (is.numeric(b) && all(b == 0)) {
+ a
+ } else if (is.numeric(a) && is.numeric(b)) {
+ a + b
+ } else {
+ expr[[2]] <- a
+ expr[[3]] <- b
+ expr
+ }
+ }
+
+ `Simplify.-` <- function(expr)
+ {
+ if (length(expr) == 2)
+ {
+ if (is.numeric(expr[[2]]))
+ return(-expr[[2]])
+ return(expr)
+ }
+
+ a <- Simplify_(expr[[2]])
+ b <- Simplify_(expr[[3]])
+
+ if (is.numeric(a) && all(a == 0)) {
+ -b
+ } else if (is.numeric(b) && all(b == 0)) {
+ a
+ } else if (is.numeric(a) && is.numeric(b)) {
+ a - b
+ } else {
+ expr[[2]] <- a
+ expr[[3]] <- b
+ expr
+ }
+ }
+
+ `Simplify.(` <- function(expr)
+ expr[[2]]
+
+ `Simplify.*` <- function(expr)
+ {
+ a <- Simplify_(expr[[2]])
+ b <- Simplify_(expr[[3]])
+
+ if (is.numeric(a) && all(a == 0)) {
+ 0
+ } else if (is.numeric(b) && all(b == 0)) {
+ 0
+ } else if (is.numeric(a) && all(a == 1)) {
+ b
+ } else if (is.numeric(b) && all(b == 1)) {
+ a
+ } else if (is.numeric(a) && is.numeric(b)) {
+ a * b
+ } else {
+ expr[[2]] <- a
+ expr[[3]] <- b
+ expr
+ }
+ }
+
+ `Simplify.^` <- function(expr)
+ {
+ a <- Simplify_(expr[[2]])
+ b <- Simplify_(expr[[3]])
+
+ if (is.numeric(a) && all(a == 0)) {
+ 0
+ } else if (is.numeric(b) && all(b == 0)) {
+ 1
+ } else if (is.numeric(a) && all(a == 1)) {
+ 1
+ } else if (is.numeric(b) && all(b == 1)) {
+ a
+ } else if (is.numeric(a) && is.numeric(b)) {
+ a ^ b
+ } else {
+ expr[[2]] <- a
+ expr[[3]] <- b
+ expr
+ }
+ }
+
+ `Simplify.c` <- function(expr)
+ {
+ args <- expr[-1]
+ args.simplified <- lapply(args, Simplify_)
+ if (all(lapply(args.simplified, is.numeric))) {
+ as.numeric(args.simplified)
+ } else {
+ for (i in 1:length(args))
+ expr[[i + 1]] <- args.simplified[[i]]
+ expr
+ }
+ }
+
+
+ assign("+", `Simplify.+`, envir=yuima.env)
+ assign("-", `Simplify.-`, envir=yuima.env)
+ assign("*", `Simplify.*`, envir=yuima.env)
+ assign("(", `Simplify.(`, envir=yuima.env)
+ assign("c", `Simplify.c`, envir=yuima.env)
+ assign("^", `Simplify.^`, envir=yuima.env)
+
+
+ ###
+
+
+
+
+
+
+
+
as.expression(Simplify_(expr[[1]]))
-
-Simplify_ <- function(expr)
-{
- if (is.symbol(expr)) {
- expr
- } else if (is.language(expr) && is.symbol(expr[[1]])) {
- # is there a rule in the table?
- sym.name <- as.character(expr[[1]])
- if (class(try(Simplify.rule <-
- get(sym.name, envir=simplifications,
- inherits=FALSE), silent=TRUE))
- != "try-error")
- return(Simplify.rule(expr))
- }
- expr
}
-Simplify.function <- function(f, x=names(formals(f)), env=parent.frame())
-{
- stopifnot(is.function(f))
- as.function(c(as.list(formals(f)),
- Simplify_(body(f))),
- envir=env)
-}
-`Simplify.+` <- function(expr)
-{
- if (length(expr) == 2)
- {
- if (is.numeric(expr[[2]]))
- return(+expr[[2]])
- return(expr)
- }
-
- a <- Simplify_(expr[[2]])
- b <- Simplify_(expr[[3]])
-
- if (is.numeric(a) && all(a == 0)) {
- b
- } else if (is.numeric(b) && all(b == 0)) {
- a
- } else if (is.numeric(a) && is.numeric(b)) {
- a + b
- } else {
- expr[[2]] <- a
- expr[[3]] <- b
- expr
- }
-}
-
-`Simplify.-` <- function(expr)
-{
- if (length(expr) == 2)
- {
- if (is.numeric(expr[[2]]))
- return(-expr[[2]])
- return(expr)
- }
-
- a <- Simplify_(expr[[2]])
- b <- Simplify_(expr[[3]])
-
- if (is.numeric(a) && all(a == 0)) {
- -b
- } else if (is.numeric(b) && all(b == 0)) {
- a
- } else if (is.numeric(a) && is.numeric(b)) {
- a - b
- } else {
- expr[[2]] <- a
- expr[[3]] <- b
- expr
- }
-}
-
-`Simplify.(` <- function(expr)
- expr[[2]]
-
-`Simplify.*` <- function(expr)
-{
- a <- Simplify_(expr[[2]])
- b <- Simplify_(expr[[3]])
-
- if (is.numeric(a) && all(a == 0)) {
- 0
- } else if (is.numeric(b) && all(b == 0)) {
- 0
- } else if (is.numeric(a) && all(a == 1)) {
- b
- } else if (is.numeric(b) && all(b == 1)) {
- a
- } else if (is.numeric(a) && is.numeric(b)) {
- a * b
- } else {
- expr[[2]] <- a
- expr[[3]] <- b
- expr
- }
-}
-
-`Simplify.^` <- function(expr)
-{
- a <- Simplify_(expr[[2]])
- b <- Simplify_(expr[[3]])
-
- if (is.numeric(a) && all(a == 0)) {
- 0
- } else if (is.numeric(b) && all(b == 0)) {
- 1
- } else if (is.numeric(a) && all(a == 1)) {
- 1
- } else if (is.numeric(b) && all(b == 1)) {
- a
- } else if (is.numeric(a) && is.numeric(b)) {
- a ^ b
- } else {
- expr[[2]] <- a
- expr[[3]] <- b
- expr
- }
-}
-
-`Simplify.c` <- function(expr)
-{
- args <- expr[-1]
- args.simplified <- lapply(args, Simplify_)
- if (all(lapply(args.simplified, is.numeric))) {
- as.numeric(args.simplified)
- } else {
- for (i in 1:length(args))
- expr[[i + 1]] <- args.simplified[[i]]
- expr
- }
-}
-
-assign("simplifications", new.env(), envir=globalenv())
-
-assign("+", `Simplify.+`, envir=simplifications)
-assign("-", `Simplify.-`, envir=simplifications)
-assign("*", `Simplify.*`, envir=simplifications)
-assign("(", `Simplify.(`, envir=simplifications)
-assign("c", `Simplify.c`, envir=simplifications)
-assign("^", `Simplify.^`, envir=simplifications)
-
-
## Constructor and Initializer of class 'yuima'
# we convert objects to "zoo" internally
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2014-03-19 14:02:57 UTC (rev 289)
+++ pkg/yuima/R/yuima.model.R 2014-03-19 16:33:41 UTC (rev 290)
@@ -128,10 +128,14 @@
time.variable="t",
solve.variable,
xinit=NULL){
+ ## we need a temp env for simplifications
+
+ yuimaENV <- new.env()
##::measure and jump term #####################################
##::initialize objects ########
MEASURE <- list()
+
##::end initialize objects ########
##::make type function list ####
@@ -388,8 +392,8 @@
}
##22/11:: function to simplify expression in drift, diffusion, jump and xinit characters
yuima.Simplifyobj<-function(x){
- dummy<-yuima.Simplify(x)
- dummy1<-yuima.Simplify(dummy)
+ dummy<-yuima.Simplify(x, yuima.env=yuimaENV)
+ dummy1<-yuima.Simplify(dummy, yuima.env=yuimaENV)
dummy2<-as.character(dummy1)
res<-parse(text=paste0("(",dummy2,")",collapse=NULL))
return(res)
More information about the Yuima-commits
mailing list