[Yuima-commits] r150 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 25 06:32:03 CEST 2011
Author: hinohide
Date: 2011-04-25 06:31:57 +0200 (Mon, 25 Apr 2011)
New Revision: 150
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/simFunctional.R
Log:
asymptotic_term bug fixed
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2011-04-08 06:46:54 UTC (rev 149)
+++ pkg/yuima/DESCRIPTION 2011-04-25 04:31:57 UTC (rev 150)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.185
-Date: 2011-04-08
+Version: 0.1.186
+Date: 2011-04-25
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/R/simFunctional.R
===================================================================
--- pkg/yuima/R/simFunctional.R 2011-04-08 06:46:54 UTC (rev 149)
+++ pkg/yuima/R/simFunctional.R 2011-04-25 04:31:57 UTC (rev 150)
@@ -1,216 +1,218 @@
-# funcF
-# function to calculate F in (13.2)
-funcF <- function(yuima,X,e,expand.var="e"){
- ##:: fix bug 07/23
- assign(expand.var, e)
-
- division <- yuima at sampling@n[1] #number of observed time
- F <- getF(yuima at functional)
- d.size <- yuima at model@equation.number
- k.size <- length(F)
- modelstate <- yuima at model@state.variable
- XT <- X[division,] #X observed last. size:vector[d.size]
- resF <- numeric(k.size) #values of F. size:vector[k.size]
-
- for(d in 1:d.size){
- assign(modelstate[d],XT[d]) #input XT in state("x1","x2") to use eval function to F
- }
- for(k in 1:k.size){
- resF[k] <- eval(F[k]) #calculate F with XT
- }
- return(resF)
- }
-
-# funcf
-# function to calculate fa in (13.2)
-funcf <- function(yuima,X,e, expand.var){
- ##:: fix bug 07/23
- assign(expand.var, e)
-
- ## division <- yuima at sampling@n #number of observed time
- division <- yuima at sampling@n[1]
- F <- getF(yuima at functional)
- f <- getf(yuima at functional)
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- k.size <- length(F)
- modelstate <- yuima at model@state.variable
- resf <- array(0,dim=c(k.size,division,(r.size+1))) #value of f0,f1,~,fr. size:array[k.size,division,r.size+1]
-
-
- for(r in 1:(r.size+1)){
- for(t in 1:division){
- Xt <- X[t,] #Xt is data X observed on time t. size:vector[d.size]
- for(d in 1:d.size){
- assign(modelstate[d],Xt[d]) #input Xt in state to use eval function to f
- }
- for(k in 1:k.size){
- resf[k,t,r] <- eval(f[[r]][k]) #calculate k th expression of fr.
- }
- }
- }
- return(resf)
- }
-
-# funcFe.
-# function to calculate Fe in (13.2).
-# core function of 'simFunctional'
-funcFe. <- function(yuima,X,e,expand.var="e"){
-
- ##:: fix bug 07/23
- assign(expand.var, e) ## expand.var
- F <- getF(yuima at functional)
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- k.size <- length(F)
- modelstate <- yuima at model@state.variable
- ## division <- yuima at sampling@n
- division <- yuima at sampling@n[1] ## 2010/11/13 modified. Currently, division must be the same for each path
- Terminal <- yuima at sampling@Terminal
- delta <- Terminal/division #length between observed times
- dw <- matrix(0,division-1,r.size+1) #Wr(t) size:matrix[division,r.size+1]
-
- dw[,1] <- rep(Terminal/(division-1),length=division-1) #W0(t)=t
-
- for(r in 2:(r.size+1)){
- tmp <- rnorm(division,0,sqrt(delta)) #calculate Wr(t)
- dw[,r] <- tmp[2:division]-tmp[1:(division-1)] #calculate dWr(t)
- }
-
- resF <- funcF(yuima,X,e,expand.var=expand.var) #calculate F with X,e. size:vector[k.size]
- resf <- funcf(yuima,X,e,expand.var=expand.var) #calculate f with X,e. size:array[k.size,division,r.size+1] ## ¤¤¤Þ, ¤³¤³¤Ç¥¨¥é¡¼ 2010/11/13
- Fe <- numeric(k.size)
- for(k in 1:k.size){
- Fe[k] <- sum(resf[k,1:division-1,]*dw)+resF[k] #calculate Fe using resF and resf as (13.2).
- }
- return(Fe)
-}
-
-# Get.X.t0
-# function to calculate X(t=t0) using runge kutta method
-Get.X.t0 <- function(yuima, expand.var="e"){
- ##:: fix bug 07/23
- assign(expand.var,0) ## epsilon=0
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- k.size <- length(getF(yuima at functional))
- modelstate <- yuima at model@state.variable
- V0 <- yuima at model@drift
- V <- yuima at model@diffusion
-
- Terminal <- yuima at sampling@Terminal[1]
-## division <- yuima at sampling@n
- division <- yuima at sampling@n[1] ## 2010/11/13
-
- ##:: fix bug 07/23
- #pars <- #yuima at model@parameter at all[1] #epsilon
-
- xinit <- getxinit(yuima at functional)
-
- ## init
- dt <- Terminal/division
- X <- xinit
- Xt <- xinit
-
- k <- numeric(d.size)
- k1 <- numeric(d.size)
- k2 <- numeric(d.size)
- k3 <- numeric(d.size)
- k4 <- numeric(d.size)
-
- ## runge kutta
- for(t in 1:division){
- ## k1
- for(i in 1:d.size){
- assign(modelstate[i],Xt[i])
- }
- for(i in 1:d.size){
- k1[i] <- dt*eval(V0[i])
- }
- ## k2
- for(i in 1:d.size){
- assign(modelstate[i],Xt[i]+k1[i]/2)
- }
- for(i in 1:d.size){
- k2[i] <- dt*eval(V0[i])
- }
- ## k3
- for(i in 1:d.size){
- assign(modelstate[i],Xt[i]+k2[i]/2)
- }
- for(i in 1:d.size){
- k3[i] <- dt*eval(V0[i])
- }
- ## k4
- for(i in 1:d.size){
- assign(modelstate[i],Xt[i]+k3[i])
- }
- for(i in 1:d.size){
- k4[i] <- dt*eval(V0[i])
- }
- ## V0(X(t+dt))
- k <- (k1+k2+k2+k3+k3+k4)/6
- Xt <- Xt+k
- X <- rbind(X,Xt)
- }
-
- ## return matrix : (division+1)*d
- rownames(X) <- NULL
- colnames(X) <- modelstate
- return(ts(X,deltat=dt[1],start=0))
- }
-
-# simFunctional
-# public function to calculate Fe in (13.2).
-setGeneric("simFunctional",
- function(yuima, expand.var="e")
- standardGeneric("simFunctional")
- )
-setMethod("simFunctional", signature(yuima="yuima"),
- function(yuima, expand.var="e"){
- Xlen <- length(yuima at data)
- if(sum(Xlen != mean(Xlen)) != 0) {
- yuima.warn("All length must be same yet.")
- return(NULL)
- }
- if( (Xlen[1]-1) != yuima at sampling@n){
- yuima.warn("Length of time series and division do not much.")
- return(NULL)
- }
-
- e <- gete(yuima at functional)
-
- ##:: fix bug 07/23
- Fe <- funcFe.(yuima,as.matrix(onezoo(yuima)),e,expand.var=expand.var)
- return(Fe)
- })
-
-# F0
-# public function to calculate Fe at e=0
-setGeneric("F0",
- function(yuima, expand.var="e")
- standardGeneric("F0")
- )
-setMethod("F0", signature(yuima="yuima"),
- function(yuima, expand.var="e"){
-
- X.t0 <- Get.X.t0(yuima, expand.var=expand.var)
- F0 <- funcFe.(yuima, X.t0, 0, expand.var=expand.var)
- return(F0)
- })
-
-# Fnorm
-# public function to calculate (Fe-F0)/e
-setGeneric("Fnorm",
- function(yuima, expand.var="e")
- standardGeneric("Fnorm")
- )
-setMethod("Fnorm", signature(yuima="yuima"),
- function(yuima, expand.var="e"){
- e <- gete(yuima at functional)
-
- Fe <- simFunctional(yuima, expand.var=expand.var)
- F0 <- F0(yuima, expand.var=expand.var)
-
- return((Fe-F0)/e)
- })
+# funcF
+# function to calculate F in (13.2)
+funcF <- function(yuima,X,e,expand.var="e"){
+ ##:: fix bug 07/23
+ assign(expand.var, e)
+
+ division <- yuima at sampling@n[1] #number of observed time
+ F <- getF(yuima at functional)
+ d.size <- yuima at model@equation.number
+ k.size <- length(F)
+ modelstate <- yuima at model@state.variable
+ XT <- X[division+1,] #X observed last. size:vector[d.size]
+ resF <- numeric(k.size) #values of F. size:vector[k.size]
+
+ for(d in 1:d.size){
+ assign(modelstate[d],XT[d]) #input XT in state("x1","x2") to use eval function to F
+ }
+ for(k in 1:k.size){
+ resF[k] <- eval(F[k]) #calculate F with XT
+ }
+ return(resF)
+ }
+
+# funcf
+# function to calculate fa in (13.2)
+funcf <- function(yuima,X,e, expand.var){
+ ##:: fix bug 07/23
+ assign(expand.var, e)
+
+ ## division <- yuima at sampling@n #number of observed time
+ division <- yuima at sampling@n[1]
+ F <- getF(yuima at functional)
+ f <- getf(yuima at functional)
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ k.size <- length(F)
+ modelstate <- yuima at model@state.variable
+ resf <- array(0,dim=c(k.size,division,(r.size+1))) #value of f0,f1,~,fr. size:array[k.size,division,r.size+1]
+
+
+ for(r in 1:(r.size+1)){
+ for(t in 1:division){
+ Xt <- X[t,] #Xt is data X observed on time t. size:vector[d.size]
+ for(d in 1:d.size){
+ assign(modelstate[d],Xt[d]) #input Xt in state to use eval function to f
+ }
+ for(k in 1:k.size){
+ resf[k,t,r] <- eval(f[[r]][k]) #calculate k th expression of fr.
+ }
+ }
+ }
+ return(resf)
+ }
+
+# funcFe.
+# function to calculate Fe in (13.2).
+# core function of 'simFunctional'
+funcFe. <- function(yuima,X,e,expand.var="e"){
+
+ ##:: fix bug 07/23
+ assign(expand.var, e) ## expand.var
+ F <- getF(yuima at functional)
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ k.size <- length(F)
+ modelstate <- yuima at model@state.variable
+ ## division <- yuima at sampling@n
+ division <- yuima at sampling@n[1] ## 2010/11/13 modified. Currently, division must be the same for each path
+ Terminal <- yuima at sampling@Terminal
+ delta <- Terminal/division #length between observed times
+ dw <- matrix(0,division,r.size+1) #Wr(t) size:matrix[division,r.size+1]
+ dw[,1] <- rep(delta,length=division) #W0(t)=t
+
+ for(r in 2:(r.size+1)){
+ tmp <- rnorm(division,0,sqrt(delta)) #calculate Wr(t)
+ tmp <- c(0,tmp)
+## dw[,r] <- tmp[-1]-tmp[-(division+1)] #calculate dWr(t)
+ dw[,r] <- diff(tmp) #calculate dWr(t)
+
+ }
+
+ resF <- funcF(yuima,X,e,expand.var=expand.var) #calculate F with X,e. size:vector[k.size]
+ resf <- funcf(yuima,X,e,expand.var=expand.var) #calculate f with X,e. size:array[k.size,division,r.size+1] ## ¤¤¤Þ, ¤³¤³¤Ç¥¨¥é¡¼ 2010/11/13
+ Fe <- numeric(k.size)
+ for(k in 1:k.size){
+ Fe[k] <- sum(resf[k,1:division,]*dw)+resF[k] #calculate Fe using resF and resf as (13.2).
+ }
+ return(Fe)
+}
+
+# Get.X.t0
+# function to calculate X(t=t0) using runge kutta method
+Get.X.t0 <- function(yuima, expand.var="e"){
+ ##:: fix bug 07/23
+ assign(expand.var,0) ## epsilon=0
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ k.size <- length(getF(yuima at functional))
+ modelstate <- yuima at model@state.variable
+ V0 <- yuima at model@drift
+ V <- yuima at model@diffusion
+
+ Terminal <- yuima at sampling@Terminal[1]
+## division <- yuima at sampling@n
+ division <- yuima at sampling@n[1] ## 2010/11/13
+
+ ##:: fix bug 07/23
+ #pars <- #yuima at model@parameter at all[1] #epsilon
+
+ xinit <- getxinit(yuima at functional)
+
+ ## init
+ dt <- Terminal/division
+ X <- xinit
+ Xt <- xinit
+
+ k <- numeric(d.size)
+ k1 <- numeric(d.size)
+ k2 <- numeric(d.size)
+ k3 <- numeric(d.size)
+ k4 <- numeric(d.size)
+
+ ## runge kutta
+ for(t in 1:division){
+ ## k1
+ for(i in 1:d.size){
+ assign(modelstate[i],Xt[i])
+ }
+ for(i in 1:d.size){
+ k1[i] <- dt*eval(V0[i])
+ }
+ ## k2
+ for(i in 1:d.size){
+ assign(modelstate[i],Xt[i]+k1[i]/2)
+ }
+ for(i in 1:d.size){
+ k2[i] <- dt*eval(V0[i])
+ }
+ ## k3
+ for(i in 1:d.size){
+ assign(modelstate[i],Xt[i]+k2[i]/2)
+ }
+ for(i in 1:d.size){
+ k3[i] <- dt*eval(V0[i])
+ }
+ ## k4
+ for(i in 1:d.size){
+ assign(modelstate[i],Xt[i]+k3[i])
+ }
+ for(i in 1:d.size){
+ k4[i] <- dt*eval(V0[i])
+ }
+ ## V0(X(t+dt))
+ k <- (k1+k2+k2+k3+k3+k4)/6
+ Xt <- Xt+k
+ X <- rbind(X,Xt)
+ }
+
+ ## return matrix : (division+1)*d
+ rownames(X) <- NULL
+ colnames(X) <- modelstate
+ return(ts(X,deltat=dt[1],start=0))
+ }
+
+# simFunctional
+# public function to calculate Fe in (13.2).
+setGeneric("simFunctional",
+ function(yuima, expand.var="e")
+ standardGeneric("simFunctional")
+ )
+setMethod("simFunctional", signature(yuima="yuima"),
+ function(yuima, expand.var="e"){
+ Xlen <- length(yuima at data)
+ if(sum(Xlen != mean(Xlen)) != 0) {
+ yuima.warn("All length must be same yet.")
+ return(NULL)
+ }
+ if( (Xlen[1]-1) != yuima at sampling@n){
+ yuima.warn("Length of time series and division do not much.")
+ return(NULL)
+ }
+
+ e <- gete(yuima at functional)
+
+ ##:: fix bug 07/23
+ Fe <- funcFe.(yuima,as.matrix(onezoo(yuima)),e,expand.var=expand.var)
+ return(Fe)
+ })
+
+# F0
+# public function to calculate Fe at e=0
+setGeneric("F0",
+ function(yuima, expand.var="e")
+ standardGeneric("F0")
+ )
+setMethod("F0", signature(yuima="yuima"),
+ function(yuima, expand.var="e"){
+
+ X.t0 <- Get.X.t0(yuima, expand.var=expand.var)
+ F0 <- funcFe.(yuima, X.t0, 0, expand.var=expand.var)
+ return(F0)
+ })
+
+# Fnorm
+# public function to calculate (Fe-F0)/e
+setGeneric("Fnorm",
+ function(yuima, expand.var="e")
+ standardGeneric("Fnorm")
+ )
+setMethod("Fnorm", signature(yuima="yuima"),
+ function(yuima, expand.var="e"){
+ e <- gete(yuima at functional)
+
+ Fe <- simFunctional(yuima, expand.var=expand.var)
+ F0 <- F0(yuima, expand.var=expand.var)
+
+ return((Fe-F0)/e)
+ })
More information about the Yuima-commits
mailing list