[Yuima-commits] r461 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 2 19:40:52 CEST 2016
Author: lorenzo
Date: 2016-09-02 19:40:52 +0200 (Fri, 02 Sep 2016)
New Revision: 461
Added:
pkg/yuima/R/AuxMethodforPPR.R
pkg/yuima/man/lambdaFromData.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/qmle.R
Log:
Added qmle for Point Process Regression Models
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2016-07-27 17:30:44 UTC (rev 460)
+++ pkg/yuima/DESCRIPTION 2016-09-02 17:40:52 UTC (rev 461)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.1.4
+Version: 1.1.5
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
@@ -9,4 +9,4 @@
Description: Simulation and Inference for Stochastic Differential Equations.
License: GPL-2
URL: http://R-Forge.R-project.org/projects/yuima/
-LinkingTo: Rcpp
\ No newline at end of file
+LinkingTo: Rcpp
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2016-07-27 17:30:44 UTC (rev 460)
+++ pkg/yuima/NAMESPACE 2016-09-02 17:40:52 UTC (rev 461)
@@ -166,6 +166,7 @@
export(rIG, rNIG, rbgamma, rvgamma, rstable,rpts,rnts) ## random number generators in rng.R
export(dIG, dNIG, dbgamma, dvgamma) ## pdfs in rng.R
export(limiting.gamma)
+export(lambdaFromData)
export(setFunctional)
Added: pkg/yuima/R/AuxMethodforPPR.R
===================================================================
--- pkg/yuima/R/AuxMethodforPPR.R (rev 0)
+++ pkg/yuima/R/AuxMethodforPPR.R 2016-09-02 17:40:52 UTC (rev 461)
@@ -0,0 +1,255 @@
+## Here we write all auxiliar functions for the Point Process
+## Regression Model
+is.Ppr <- function(yuimaPpr){is(yuimaPpr,"yuima.Ppr")}
+
+quasiLogLik.Ppr <- function(yuimaPpr, parLambda=list(), method=method, fixed = list(),
+ lower, upper, call, ...){
+ PprData<-yuimaPpr at data
+ Time <- index(yuimaPpr at data@zoo.data[[1]])
+ envPpr <- list()
+ dY <- paste0("d",yuimaPpr at Ppr@var.dx)
+ for(i in (c(1:(length(Time)-1)))){
+ envPpr[[i]]<-new.env()
+ assign(yuimaPpr at gFun@param at time.var,rep(Time[i+1],i),envir=envPpr[[i]])
+ assign(yuimaPpr at Ppr@var.dt,Time[1:i],envir=envPpr[[i]])
+ if(length(yuimaPpr at Ppr@covariates)>0){
+ for(j in c(1:length(yuimaPpr at Ppr@covariates))){
+ cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@covariates[[j]]
+ assign(yuimaPpr at Ppr@covariates[[j]],
+ as.numeric(yuimaPpr at data@original.data[1:(i+1),cond]),
+ envir=envPpr[[i]])
+ }
+ }
+ for(j in c(1:length(yuimaPpr at Ppr@counting.var))){
+ cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@counting.var[[j]]
+ assign(yuimaPpr at Ppr@counting.var[[j]],
+ as.numeric(yuimaPpr at data@original.data[1:(i+1),cond]),
+ envir=envPpr[[i]])
+ }
+
+ for(j in c(1:length(yuimaPpr at Ppr@var.dx))){
+ cond<-c(colnames(yuimaPpr at data@original.data),yuimaPpr at Ppr@var.dt)%in%c(yuimaPpr at Ppr@var.dx,yuimaPpr at Ppr@var.dt)[[j]]
+ if(any(cond[-length(cond)])){
+ assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
+ diff(as.numeric(yuimaPpr at data@original.data[1:(i+1),cond[-length(cond)]])),
+ envir=envPpr[[i]])
+ }
+ if(tail(cond,n=1L)){
+ assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
+ as.numeric(diff(Time[1:(i+1),cond[-length(cond)]])),
+ envir=envPpr[[i]])
+ }
+ }
+ }
+ IntKernExpr<- function(Kern,dy){
+ dum<-paste(Kern,dy,sep="*")
+ dum<-paste0(dum, collapse = " + ")
+ dum <- paste0("sum( ", dum, " )")
+ return(parse(text = dum))
+ }
+
+ IntegKern <- lapply(yuimaPpr at Kernel@Integrand at IntegrandList,IntKernExpr,dY)
+ Integrator <- t(as.matrix(eval(parse(text=dY[1]),envir=envPpr[[length(envPpr)]])))
+ if(length(dY)>1){
+ for(i in c(2:length(dY))){
+ Integrator <- rbind(Integrator,
+ t(as.matrix(eval(parse(text=dY[1]),envir=envPpr[[length(envPpr)]]))))
+ }
+ }
+ assign("Integrator",Integrator,envir=envPpr[[length(envPpr)]])
+ assign("Nlamb",length(yuimaPpr at Ppr@counting.var),envir=envPpr[[length(envPpr)]])
+
+ out<-NULL
+ param1<-unlist(parLambda)
+ my.env <- envPpr
+ if(length(lower)==0 && length(upper)>0 && length(fixed)==0){
+ out <- optim(par=param1, fn=aux.lambdaFromData,
+ method = method, upper=upper, envPpr = my.env, gFun=yuimaPpr at gFun,
+ Kern =IntegKern, intensityParm = yuimaPpr at Ppr@allparamPpr,
+ logLikelihood=TRUE, ...)
+
+ }
+
+ if(length(lower)==0 && length(upper)==0 && length(fixed)>0){
+ out <- optim(par = param1, fn = aux.lambdaFromData,
+ method = method, fixed = fixed, envPpr = my.env,
+ gFun = yuimaPpr at gFun, Kern = IntegKern,
+ intensityParm = yuimaPpr at Ppr@allparamPpr,
+ logLikelihood = TRUE, ...)
+
+ }
+
+
+ if(length(lower)>0 && length(upper)==0 && length(fixed)==0){
+ out <- optim(par = param1, fn = aux.lambdaFromData,
+ method = method, lower=lower, envPpr = my.env,
+ gFun = yuimaPpr at gFun, Kern = IntegKern,
+ intensityParm = yuimaPpr at Ppr@allparamPpr,
+ logLikelihood = TRUE, ...)
+ }
+
+ if(length(lower)>0 && length(upper)>0 && length(fixed)==0){
+ out <- optim(par = param1, fn = aux.lambdaFromData,
+ method = method, upper = upper,
+ lower=lower, envPpr = my.env,
+ gFun = yuimaPpr at gFun, Kern = IntegKern,
+ intensityParm = yuimaPpr at Ppr@allparamPpr,
+ logLikelihood = TRUE, ...)
+ }
+
+
+ if(length(lower)==0 && length(upper)>0 && length(fixed)>0){
+ out <- optim(par = param1, fn = aux.lambdaFromData,
+ method = method, upper = upper,
+ fixed = fixed, envPpr = my.env,
+ gFun = yuimaPpr at gFun, Kern = IntegKern,
+ intensityParm = yuimaPpr at Ppr@allparamPpr,
+ logLikelihood = TRUE, ...)
+ }
+
+ if(length(lower)>0 && length(upper)==0 && length(fixed)>0){
+ out <- optim(par=param1, fn=aux.lambdaFromData,
+ method = method, lower = lower,
+ fixed = fixed, envPpr = my.env,
+ gFun = yuimaPpr at gFun, Kern = IntegKern,
+ intensityParm = yuimaPpr at Ppr@allparamPpr,
+ logLikelihood = TRUE, ...)
+ }
+
+
+ if(length(lower)>0 && length(upper)>0 && length(fixed)>0){
+ out <- optim(par = param1, fn = aux.lambdaFromData,
+ method = method, lower = lower, fixed = fixed, upper = upper,
+ envPpr = my.env, gFun = yuimaPpr at gFun, Kern = IntegKern,
+ intensityParm = yuimaPpr at Ppr@allparamPpr, logLikelihood = TRUE, ...)
+ }
+
+
+ if(is.null(out)){
+ out <- optim(par = param1, fn = aux.lambdaFromData,
+ method = method, envPpr = my.env, gFun = yuimaPpr at gFun,
+ Kern = IntegKern, intensityParm = yuimaPpr at Ppr@allparamPpr,
+ logLikelihood = TRUE, ...)
+ }
+
+
+ return(out)
+
+}
+
+lambdaFromData <- function(yuimaPpr, PprData=NULL, parLambda=list()){
+ if(is.null(PprData)){
+ PprData<-yuimaPpr at data
+ }else{
+ # checklambdaFromData(yuimaPpr,PprData)
+ }
+ if(!any(names(parLambda) %in% yuimaPpr at Ppr@allparamPpr)){yuima.stop("1 ...")}
+ if(!any(yuimaPpr at Ppr@allparamPpr %in% names(parLambda))){yuima.stop("2 ...")}
+ Time <- index(yuimaPpr at data@zoo.data[[1]])
+ envPpr <- list()
+ dY <- paste0("d",yuimaPpr at Ppr@var.dx)
+ for(i in (c(1:(length(Time)-1)))){
+ envPpr[[i]]<-new.env()
+ assign(yuimaPpr at gFun@param at time.var,rep(Time[i+1],i),envir=envPpr[[i]])
+ assign(yuimaPpr at Ppr@var.dt,Time[1:i],envir=envPpr[[i]])
+ if(length(yuimaPpr at Ppr@covariates)>0){
+ for(j in c(1:length(yuimaPpr at Ppr@covariates))){
+ cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@covariates[[j]]
+ assign(yuimaPpr at Ppr@covariates[[j]],
+ as.numeric(yuimaPpr at data@original.data[1:(i+1),cond]),
+ envir=envPpr[[i]])
+ }
+ }
+ for(j in c(1:length(yuimaPpr at Ppr@counting.var))){
+ cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@counting.var[[j]]
+ assign(yuimaPpr at Ppr@counting.var[[j]],
+ as.numeric(yuimaPpr at data@original.data[1:(i+1),cond]),
+ envir=envPpr[[i]])
+ }
+
+ for(j in c(1:length(yuimaPpr at Ppr@var.dx))){
+ cond<-c(colnames(yuimaPpr at data@original.data),yuimaPpr at Ppr@var.dt)%in%c(yuimaPpr at Ppr@var.dx,yuimaPpr at Ppr@var.dt)[[j]]
+ if(any(cond[-length(cond)])){
+ assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
+ diff(as.numeric(yuimaPpr at data@original.data[1:(i+1),cond[-length(cond)]])),
+ envir=envPpr[[i]])
+ }
+ if(tail(cond,n=1L)){
+ assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
+ as.numeric(diff(Time[1:(i+1),cond[-length(cond)]])),
+ envir=envPpr[[i]])
+ }
+ }
+ }
+ IntKernExpr<- function(Kern,dy){
+ dum<-paste(Kern,dy,sep="*")
+ dum<-paste0(dum, collapse = " + ")
+ dum <- paste0("sum( ", dum, " )")
+ return(parse(text = dum))
+ }
+
+ IntegKern <- lapply(yuimaPpr at Kernel@Integrand at IntegrandList,IntKernExpr,dY)
+ Integrator <- t(as.matrix(eval(parse(text=dY[1]),envir=envPpr[[length(envPpr)]])))
+ if(length(dY)>1){
+ for(i in c(2:length(dY))){
+ Integrator <- rbind(Integrator,
+ t(as.matrix(eval(parse(text=dY[1]),envir=envPpr[[length(envPpr)]]))))
+ }
+ }
+ assign("Integrator",Integrator,envir=envPpr[[length(envPpr)]])
+ assign("Nlamb",length(yuimaPpr at Ppr@counting.var),envir=envPpr[[length(envPpr)]])
+ res<-aux.lambdaFromData(param = unlist(parLambda), gFun=yuimaPpr at gFun,
+ Kern =IntegKern, intensityParm = yuimaPpr at Ppr@allparamPpr,
+ envPpr)
+ return(res)
+}
+# my.lapply <- function (X, FUN, ...){
+# # FUN <- match.fun(FUN)
+# .Internal(lapply(X, FUN))
+# }
+myfun3<-function(X,Kern){
+ t(simplify2array(lapply(Kern, FUN=DumFun,
+ Y = X), higher = (TRUE == "array")))}
+dumFun2<-function(X,Y){list2env(Y,envir=X)}
+
+aux.lambdaFromData <-function(param, gFun, Kern, intensityParm, envPpr,logLikelihood = FALSE){
+ lapply(envPpr,FUN=dumFun2,Y=as.list(param))
+ lastEnv <- tail(envPpr,n=1L)[[1]]
+ # gFunVect<- t(simplify2array(my.lapply(gFun at formula, FUN=DumFun,
+ # Y = lastEnv), higher = (TRUE == "array")))
+
+ gFunVect<- matrix(unlist(lapply(gFun at formula, FUN=DumFun,
+ Y = lastEnv)),nrow=lastEnv$Nlamb,byrow=TRUE)
+ # IntKer<- simplify2array(my.lapply(envPpr,function(my.env){
+ # t(simplify2array(my.lapply(Kern, FUN=DumFun,
+ # Y = my.env), higher = (TRUE == "array")))}),
+ # higher = (TRUE == "array")
+ # )
+ # IntKer<- simplify2array(my.lapply(envPpr,myfun3,Kern=Kern),
+ # higher = (TRUE == "array")
+ # )
+
+ # IntKer<- matrix(unlist(my.lapply(envPpr,myfun3,Kern=Kern)),
+ # nrow=1,byrow=TRUE)
+ IntKer<- matrix(unlist(lapply(envPpr,myfun3,Kern=Kern)),
+ nrow=lastEnv$Nlamb)
+ lambda <- gFunVect+cbind(0,IntKer)
+ time <- (c(lastEnv$s,lastEnv$t[1]))
+ if(!logLikelihood){
+ Intensity <- zoo(t(lambda), order.by = time)
+ return(Intensity)
+ }
+ dn <- dim(lambda)
+ if(dn[1]==1){
+ logLiklihood2 <- -sum(lambda[,-1]*diff(time)[1])
+ logLiklihood1 <- sum(log(lambda[,-1])*lastEnv$Integrator)
+ }else{
+ logLiklihood2 <- -rowSums(lambda[,-1]*diff(time)[1])
+ logLiklihood1 <- rowSums(log(lambda[,-1])*lastEnv$Integrator)
+ }
+ minusLoglik <- -sum(logLiklihood2+logLiklihood1)
+ #cat(sprintf("\n%.5f",minusLoglik))
+ return(minusLoglik)
+}
+
+DumFun<- function(X,Y){eval(X,envir=Y)}
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2016-07-27 17:30:44 UTC (rev 460)
+++ pkg/yuima/R/qmle.R 2016-09-02 17:40:52 UTC (rev 461)
@@ -153,6 +153,25 @@
return(res)
}
+ if(is.Ppr(yuima)){
+ if(missing(lower))
+ lower <- list()
+
+ if(missing(upper))
+ upper <- list()
+
+ # res <- NULL
+ # if("grideq" %in% names(as.list(call)[-(1:2)])){
+ res <- quasiLogLik.Ppr(yuimaPpr = yuima, parLambda = start, method=method, fixed = list(),
+ lower, upper, call, ...)
+ # }else{
+ # res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
+ # lower, upper, Est.Incr, call, grideq = FALSE, aggregation = aggregation,...)
+ # }
+
+ return(res)
+ }
+
orig.fixed <- fixed
orig.fixed.par <- names(orig.fixed)
if(is.Poisson(yuima))
@@ -1810,9 +1829,9 @@
####r <- length(a[[1]])
r <- yuima at model@noise.number
xdim <- length(yuima at model@state.variable)
-
+
#if(thetadim!=length(initial)) stop("check dim of initial") #error check
-
+
d_b <- NULL
for(i in 1:d){
check_x <- NULL
@@ -1825,7 +1844,7 @@
}
}
#d_b <- c(d_b,b[[i]])
-
+
v_a<-matrix(list(NULL),d,r)
for(i in 1:d){
for(j in 1:r){
@@ -1839,7 +1858,7 @@
}
}
}
-
+
for(i in 1:d) assign(yuima at model@state.variable[i], data[-length(data[,1]),i])
dx <- as.matrix((data-rbind(numeric(d),as.matrix(data[-length(data[,1]),])))[-1,])
drift <- diffusion <- NULL
Added: pkg/yuima/man/lambdaFromData.Rd
===================================================================
--- pkg/yuima/man/lambdaFromData.Rd (rev 0)
+++ pkg/yuima/man/lambdaFromData.Rd 2016-09-02 17:40:52 UTC (rev 461)
@@ -0,0 +1,33 @@
+\name{lambdaFromData}
+\alias{lambdaFromData}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+Intensity of a Point Process Regression Model
+}
+\description{This function returns the intensity process of a Ppr model when covariates and counting processes are obsered on discrete time}
+\usage{
+lambdaFromData(yuimaPpr, PprData = NULL, parLambda = list())
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{yuimaPpr}{ Mathematical Description of Ppr model
+}
+ \item{PprData}{Observed data}
+ \item{parLambda}{Values of intesity parameters}
+}
+\details{ ... }
+\value{...}
+\references{
+...}
+\author{YUIMA TEAM}
+\note{...}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\seealso{
+...}
+%\examples{}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+%\keyword{ ~kwd1 }% use one of RShowDoc("KEYWORDS")
+%\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
More information about the Yuima-commits
mailing list