[Yuima-commits] r680 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 23 18:47:45 CEST 2018
Author: lorenzo
Date: 2018-10-23 18:47:45 +0200 (Tue, 23 Oct 2018)
New Revision: 680
Modified:
pkg/yuima/R/AuxMethodforPPR.R
pkg/yuima/R/FunctionAndOperators.R
pkg/yuima/R/lambdaPPR.R
pkg/yuima/R/simulateForPpr.R
Log:
Modified: pkg/yuima/R/AuxMethodforPPR.R
===================================================================
--- pkg/yuima/R/AuxMethodforPPR.R 2018-10-04 11:21:04 UTC (rev 679)
+++ pkg/yuima/R/AuxMethodforPPR.R 2018-10-23 16:47:45 UTC (rev 680)
@@ -5,10 +5,13 @@
Internal.LogLikPPR <- function(param,my.envd1=NULL,
my.envd2=NULL,my.envd3=NULL){
param<-unlist(param)
-
- IntLambda<-InternalConstractionIntensity2(param,my.envd1,
+ if(my.envd3$CondIntensityInKern){
+ IntLambda<- InternalConstractionIntensityFeedBackIntegrand(param,my.envd1,
+ my.envd2,my.envd3)
+ }else{
+ IntLambda<-InternalConstractionIntensity2(param,my.envd1,
my.envd2,my.envd3)
-
+ }
# IntLambda<-InternalConstractionIntensity(param,my.envd1,
# my.envd2,my.envd3)
Index<-my.envd3$gridTime
@@ -62,7 +65,7 @@
#+sum((param-oldpar)^2*param^2)/2
# line 40 necessary for the development of the cod
- #cat("\n ",logLik, param)
+ cat("\n ",logLik, param)
#assign("oldpar",param,envir = my.envd1)
@@ -308,9 +311,11 @@
assign("Univariate",Univariate,envir=my.envd3)
assign("ExistdN",ExistdN,envir=my.envd3)
assign("ExistdX",ExistdX,envir=my.envd3)
+ assign("JumpTimeLogical",c(FALSE,as.integer(diff(my.envd3$N))!=0),envir=my.envd3)
+ assign("CondIntensityInKern",
+ my.envd3$YUIMA.PPR at PPR@additional.info %in% all.vars(my.envd3$Integrand2expr),
+ envir=my.envd3)
-
-
out<-NULL
if(length(lower)==0 && length(upper)>0 && length(fixed)==0){
Modified: pkg/yuima/R/FunctionAndOperators.R
===================================================================
--- pkg/yuima/R/FunctionAndOperators.R 2018-10-04 11:21:04 UTC (rev 679)
+++ pkg/yuima/R/FunctionAndOperators.R 2018-10-23 16:47:45 UTC (rev 680)
@@ -270,6 +270,13 @@
}
paramIntegrand <- paramIntegrand[!Cond]
+ # state variable
+ Cond <- (paramIntegrand %in% mod at state.variable)
+ if(sum(Cond)==0){
+ yuima.warn("Integrand fuction does not depend on state.variable")
+ }
+
+ paramIntegrand <- paramIntegrand[!Cond]
# time variable
Cond <- (paramIntegrand %in% mod at time.variable)
paramIntegrand <- paramIntegrand[!Cond]
Modified: pkg/yuima/R/lambdaPPR.R
===================================================================
--- pkg/yuima/R/lambdaPPR.R 2018-10-04 11:21:04 UTC (rev 679)
+++ pkg/yuima/R/lambdaPPR.R 2018-10-23 16:47:45 UTC (rev 680)
@@ -85,6 +85,160 @@
}
+InternalConstractionIntensityFeedBackIntegrand<-function(param,my.envd1,
+ my.envd2,my.envd3){
+ paramPPR <- my.envd3$YUIMA.PPR at PPR@allparamPPR
+ namesparam <-my.envd3$namesparam
+
+
+ gridTime <-my.envd3$gridTime
+ Univariate <-my.envd3$Univariate
+ ExistdN <-my.envd3$ExistdN
+ ExistdX <-my.envd3$ExistdX
+
+ gfun<-my.envd3$gfun
+
+ allVarsInG<- all.vars(gfun)
+ CondIntFeedBacksToG <- my.envd3$YUIMA.PPR at PPR@additional.info %in% allVarsInG
+
+ Integrand2<-my.envd3$Integrand2
+ Integrand2expr<-my.envd3$Integrand2expr
+
+ if(ExistdN){
+ for(i in c(1:length(paramPPR))){
+ cond<-namesparam %in% paramPPR[i]
+ assign(paramPPR[i], param[cond], envir = my.envd1 )
+ }
+ }
+
+ if(ExistdX){
+ for(i in c(1:length(paramPPR))){
+ cond<-namesparam %in% paramPPR[i]
+ assign(paramPPR[i], param[cond], envir = my.envd2)
+ }
+ }
+
+ #param
+ for(i in c(1:length(paramPPR))){
+ cond<-namesparam %in% paramPPR[i]
+ assign(paramPPR[i], param[cond], envir = my.envd3)
+ }
+
+
+ # for(i in c(1:length(gridTime))){
+ # KerneldN[i] <- InternalKernelFromPPRModel(Integrand2,Integrand2expr,my.envd1=my.envd1,my.envd2=my.envd2,
+ # Univariate=Univariate, ExistdN, ExistdX, gridTime=gridTime[i])
+ # }
+
+ if(Univariate){
+ NameCol <- colnames(Integrand2)
+
+ # Kernel <- sapply(X=gridTime,FUN = InternalKernelFromPPRModel3,
+ # Integrand2=Integrand2, Integrand2expr = Integrand2expr,my.envd1=my.envd1,my.envd2=my.envd2,
+ # my.envd3=my.envd3,
+ # Univariate=Univariate, ExistdN =ExistdN, ExistdX=ExistdX,
+ # dimCol=dim(Integrand2)[2], NameCol = NameCol,
+ # JumpTimeName =paste0("JumpTime.",NameCol))
+ # Kernel <- evalKernelCpp2(Integrand2,
+ # Integrand2expr,
+ # my.envd1, my.envd2, my.envd3$YUIMA.PPR at PPR@IntensWithCount,
+ # my.envd3$YUIMA.PPR at PPR@counting.var,
+ # my.envd3$YUIMA.PPR at PPR@covariates,
+ # ExistdN, ExistdX,
+ # gridTime, dimCol = dim(Integrand2)[2], NameCol = NameCol,
+ # JumpTimeName =paste0("JumpTime.",NameCol))
+ # Evalgfun <- internalGfunFromPPRModel(gfun[i],my.envd3, univariate=TRUE)
+ # result<-Kernel+Evalgfun
+ kernel <- numeric(length = length(gridTime))
+ Intensity <- numeric(length = length(gridTime))
+ JumpTimeName <- paste0("JumpTime.", NameCol)
+ dimCol <- dim(Integrand2)[2]
+ IntensityatJumpTime<-as.list(numeric(length=length(gridTime)))
+ differentialCounting <- paste0("d",
+ c(my.envd3$YUIMA.PPR at PPR@counting.var,
+ my.envd3$YUIMA.PPR at PPR@additional.info))
+ if(!CondIntFeedBacksToG){
+ EvalGFUN <- eval(gfun,envir=my.envd3)
+ Intensity[1]<-EvalGFUN[1]
+ }
+ aaaa<- length(gridTime)
+ CondallJump <- rep(FALSE,aaaa+1)
+
+ for(t in c(2:aaaa)){
+ assign(my.envd1$t.time,gridTime[[t]][1],envir=my.envd1)
+ for(j in c(1:dimCol)){
+ if(NameCol[j] %in% differentialCounting){
+ # Intensity at Jump Time
+ assign(my.envd3$YUIMA.PPR at PPR@additional.info,
+ IntensityatJumpTime[[t-1]],
+ envir=my.envd1)
+ # Counting Var at Jump Time
+ assign(my.envd3$YUIMA.PPR at PPR@counting.var,
+ my.envd1[[my.envd1$PosListCountingVariable]][[t]],
+ envir=my.envd1)
+ # Jump time <= t
+ assign(my.envd1$var.time,my.envd1[[JumpTimeName[j]]][[t]],envir=my.envd1)
+ KerneldN <- sum(eval(Integrand2expr,envir=my.envd1)*my.envd1[[my.envd1$namedX[j]]][[t]])
+ kernel[t-1] <- KerneldN
+ }
+
+ }
+ # Evaluation gFun
+ if(!CondIntFeedBacksToG){
+ #EvalGFUN <- eval(gfun,envir=my.envd3)
+ if(t<=aaaa){
+ Intensity[t]<- EvalGFUN[t-1]+kernel[t-1]
+ }
+ }else{
+ # Here we evaluate gFun time by time
+
+ }
+
+ for(j in c(1:dimCol)){
+ if(t+1<=aaaa){
+ if(NameCol[j] %in% differentialCounting){
+ #
+ CondallJump[t] <-my.envd3$JumpTimeLogical[t]
+ IntensityatJumpTime[[t]]<- Intensity[CondallJump[-1]]
+
+ }
+ }
+ }
+ if(t==77){
+ ff<-2
+ }
+ }
+
+ }else{
+ n.row <- length(my.envd3$YUIMA.PPR at PPR@counting.var)
+ n.col <- length(gridTime)
+ result <- matrix(NA,n.row, n.col)
+ #Kernel<- numeric(length=n.col-1)
+
+ dimCol<- dim(Integrand2)[2]
+ NameCol<-colnames(Integrand2)
+ JumpTimeName <- paste0("JumpTime.",NameCol)
+
+ for(i in c(1:n.row)){
+
+ # Kernel <- sapply(X=gridTime,FUN = InternalKernelFromPPRModel2,
+ # Integrand2=t(Integrand2[i,]), Integrand2expr = Integrand2expr[[i]],my.envd1=my.envd1,my.envd2=my.envd2,
+ # Univariate=TRUE, ExistdN =ExistdN, ExistdX=ExistdX, dimCol=dimCol, NameCol = NameCol,
+ # JumpTimeName =JumpTimeName)
+ # Kernel <- evalKernelCpp2(t(Integrand2[i,]), Integrand2expr[[i]],
+ # my.envd1, my.envd2, my.envd3$YUIMA.PPR at PPR@IntensWithCount,
+ # my.envd3$YUIMA.PPR at PPR@counting.var,
+ # my.envd3$YUIMA.PPR at PPR@covariates,
+ # ExistdN, ExistdX,
+ # gridTime, dimCol = dim(Integrand2)[2], NameCol = NameCol,
+ # JumpTimeName =paste0("JumpTime.",NameCol))
+
+ # Evalgfun <- internalGfunFromPPRModel(gfun[i],my.envd3, univariate=TRUE)
+ # result[i,]<-Kernel+Evalgfun
+ }
+ }
+ return(Intensity)
+}
InternalKernelFromPPRModel2<-function(Integrand2,Integrand2expr,my.envd1=NULL,my.envd2=NULL,
Univariate=TRUE, ExistdN, ExistdX, gridTime, dimCol, NameCol,
@@ -575,8 +729,9 @@
assign("Univariate",Univariate,envir=my.envd3)
assign("ExistdN",ExistdN,envir=my.envd3)
assign("ExistdX",ExistdX,envir=my.envd3)
+
+ assign("JumpTimeLogical",c(FALSE,as.integer(diff(my.envd3$N))!=0),envir=my.envd3)
-
# end construction my.envd3
@@ -589,8 +744,13 @@
#We define the parameters value for each environment
param<-unlist(param)
- result<-InternalConstractionIntensity2(param,my.envd1,
+ if(my.envd3$YUIMA.PPR at PPR@additional.info %in% all.vars(my.envd3$Integrand2expr)){
+ result <- InternalConstractionIntensityFeedBackIntegrand(param,my.envd1,
+ my.envd2,my.envd3)
+ }else{
+ result<-InternalConstractionIntensity2(param,my.envd1,
my.envd2,my.envd3)
+ }
if(class(result)=="matrix"){
my.matr <- t(result)
colnames(my.matr) <-paste0("lambda",c(1:yuimaPPR at Kernel@Integrand at dimIntegrand[1]))
Modified: pkg/yuima/R/simulateForPpr.R
===================================================================
--- pkg/yuima/R/simulateForPpr.R 2018-10-04 11:21:04 UTC (rev 679)
+++ pkg/yuima/R/simulateForPpr.R 2018-10-23 16:47:45 UTC (rev 680)
@@ -52,7 +52,7 @@
}
)
-constHazIntPr <- function(g.Fun , Kern.Fun, covariates, counting.var){
+constHazIntPr <- function(g.Fun , Kern.Fun, covariates, counting.var,statevar=NULL){
numb.Int <- length(g.Fun)
Int.Intens <- list()
dum.g <- character(length=numb.Int)
@@ -90,6 +90,33 @@
dum.Ker <- as.character(unlist(Kern.Fun at Integrand@IntegrandList))
dum.Ker <- gsub("(", "( ", fixed=TRUE,x = dum.Ker)
dum.Ker <- gsub(")", " )", fixed=TRUE,x = dum.Ker)
+
+ dimKernel<-length(dum.Ker)
+ condIntInKer <- FALSE
+ for(i in c(1:dimKernel)){
+ if(!condIntInKer)
+ condIntInKer <- statevar%in%all.vars(parse(text=dum.Ker[i]))
+ }
+ if(condIntInKer){
+ for(i in c(1:length(statevar))){
+ my.countOld <- paste0(statevar[i] ," ")
+ my.countNew <- paste0( statevar[i] ,
+ "ForKernel[CondJumpGrid]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+ my.countOld <- paste0(statevar[i] ,"[",Kern.Fun at variable.Integral@upper.var,"]")
+ # my.countNew <- paste0( counting.var[i] ,
+ # "[ as.character( ",Kern.Fun at variable.Integral@upper.var ," ) ]")
+ my.countNew <- paste0( "tail(",statevar[i] ,"ForKernel ,n=1L) ")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+
+ my.countOld <- paste0(statevar[i] ,"[",Kern.Fun at variable.Integral@var.time,"]")
+ my.countNew <- paste0(statevar[i] ,
+ "ForKernel[CondJumpGrid]")
+ dum.Ker <- gsub(my.countOld, my.countNew, x = dum.Ker, fixed=TRUE)
+
+ }
+ }
+
if(length(counting.var)>0){
for(i in c(1:length(counting.var))){
@@ -814,7 +841,7 @@
Kern <- object at Kernel
object at sampling <- samp
randomGenerator<-object at model@measure$df
-
+ nameIntensityProc <- object at PPR@additional.info
if(missing(increment.W) | is.null(increment.W)){
if( Model at hurst!=0.5 ){
@@ -950,15 +977,19 @@
}
}
CovariateSim <-matrix(0,Model at equation.number,(samp at n+1))
-
+ #IntensityProcInter <- matrix(0,length(nameIntensityProc),(samp at n+1))
ExprHaz <- constHazIntPr(g.Fun = object at gFun@formula,
Kern.Fun = object at Kernel, covariates = object at PPR@covariates,
- counting.var = object at PPR@counting.var)$Intens
-
+ counting.var = object at PPR@counting.var, statevar=nameIntensityProc)$Intens
+ IntensityProcInter <- as.matrix(tryCatch(eval(object at gFun@formula,envir=myenv),error =function(){1}))
dN <- matrix(0,object at gFun@dimension[1],object at gFun@dimension[2])
+ rownames(CovariateSim)<- Model at solve.variable
+ assign(object at PPR@counting.var,CovariateSim[object at PPR@counting.var,1],envir=myenv)
grid <- samp at grid[[1]]
const <- -log(runif(gFun at dimension[1]))
condMyTR <- const<delta
+ AllnameIntensityProc <- paste0(nameIntensityProc,"ForKernel")
+ assign(AllnameIntensityProc,IntensityProcInter,envir=myenv)
while(any(condMyTR)){
if(sum(condMyTR)==0){
const <- -log(runif(length(condMyTR)))
@@ -984,7 +1015,14 @@
lambda<-compErrHazR4(samp, Kern, capitalTime=samp at grid[[1]][inter_i[j]],
Model, myenv, ExprHaz, Time=jumpT, dN, Index, j)
if(is.matrix(posLambda)){}else{
- assign(object at model@state.variable[posLambda],lambda, envir = myenv)
+ #assign(object at model@state.variable[posLambda],lambda, envir = myenv)
+ assign(nameIntensityProc,lambda[j], envir = myenv)
+ # myenv[[AllnameIntensityProc]][j,]<-cbind(myenv[[AllnameIntensityProc]][j,],
+ # lambda[j])
+ assign(AllnameIntensityProc,
+ cbind(t(myenv[[AllnameIntensityProc]][j,]),
+ lambda[j]),
+ envir=myenv)
}
@@ -1010,18 +1048,20 @@
Terminal=samp at grid[[1]][inter_i[j]],n=1,
env=myenv)[,-1]
}
+ # if(inter_i[j]==66){
+ # aaaaa<-1
+ # }
-
if(inter_i[j]>=(dimGrid)){
noExit[j] <- FALSE
}
if(inter_i[j]<=dimGrid){
+ assign(object at PPR@counting.var,CovariateSim[object at PPR@counting.var[j],1:inter_i[j]],envir=myenv)
dimCov <- length(object at PPR@covariates)
-
if (dimCov>0){
- for(j in c(1:dimCov)){
- assign(object at PPR@covariates[j],
- as.numeric(CovariateSim[object at PPR@covariates[j],1:inter_i[j]]),
+ for(jj in c(1:dimCov)){
+ assign(object at PPR@covariates[jj],
+ as.numeric(CovariateSim[object at PPR@covariates[jj],1:inter_i[j]]),
envir = myenv)
}
}
@@ -1045,7 +1085,7 @@
Noise.Laux[object at PPR@counting.var,i-1] <- dumdN[i==inter_i]
dN <- cbind(dN,dumdN)
}
- cat("\n ", i, grid[i])
+ # cat("\n ", i, grid[i])
# assign("dL",Noise.Laux,myenv)
#
@@ -1825,7 +1865,7 @@
}
jumpT<-c(jumpT,next_jump)
- cat("\n ", c(next_jump, checkFunDum,count))
+ # cat("\n ", c(next_jump, checkFunDum,count))
u_bar <- tail(jumpT,1L)
posInitial<- sum(grid<=next_jump)
More information about the Yuima-commits
mailing list