[Yuima-commits] r633 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 29 18:49:46 CET 2017
Author: lorenzo
Date: 2017-11-29 18:49:46 +0100 (Wed, 29 Nov 2017)
New Revision: 633
Added:
pkg/yuima/R/DataPPR.R
pkg/yuima/man/DataPpr.Rd
pkg/yuima/man/get.counting.data.RD
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/AuxMethodforPPR.R
Log:
Upload Utilities for yuimaPPR
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2017-11-27 14:42:36 UTC (rev 632)
+++ pkg/yuima/DESCRIPTION 2017-11-29 17:49:46 UTC (rev 633)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.7.4
+Version: 1.7.5
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2)
Author: YUIMA Project Team
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2017-11-27 14:42:36 UTC (rev 632)
+++ pkg/yuima/NAMESPACE 2017-11-29 17:49:46 UTC (rev 633)
@@ -196,6 +196,7 @@
export(lse)
export(qmle)
+export(get.counting.data)
export(quasilogl)
export(phi.test)
export(lasso)
@@ -204,6 +205,7 @@
export(qmleL)
export(qmleLevy)
export(IC)
+export(DataPpr)
export(Intensity.PPR)
@@ -246,10 +248,6 @@
S3method(plot, yuima.llag) # Oct. 10, 2015
S3method(plot, yuima.mllag) # Oct. 10, 2015
-
useDynLib(yuima)
-#,.registration = TRUE)
-
-
Modified: pkg/yuima/R/AuxMethodforPPR.R
===================================================================
--- pkg/yuima/R/AuxMethodforPPR.R 2017-11-27 14:42:36 UTC (rev 632)
+++ pkg/yuima/R/AuxMethodforPPR.R 2017-11-29 17:49:46 UTC (rev 633)
@@ -21,7 +21,7 @@
if(length(my.envd3$YUIMA.PPR at Ppr@counting.var)>0){
cond1 <- my.envd3$YUIMA.PPR at model@solve.variable %in% my.envd3$YUIMA.PPR at Ppr@counting.var
cond2 <- diff(as.numeric(my.envd3$YUIMA.PPR at data@original.data[,cond1]))
- #Integr2a<- sum(log(IntLambda[-1][cond2!=0]),na.rm=TRUE)
+ #Integr2<- sum(log(IntLambda[-1][cond2!=0]),na.rm=TRUE)
Integr2 <- sum(log(IntLambda[cond2!=0]),na.rm=TRUE)
#Integr2 <- (Integr2a+Integr2b)/2
}else{
@@ -31,8 +31,17 @@
# Integr2 <- -10^6
# }
logLik <- Integr1+Integr2
+ if(is.null(my.envd1$oldpar)){
+ oldpar <- param
+ }else{
+ oldpar <- my.envd1$oldpar
+ }
+ ret <- -logLik/sum(cond2,na.rm=TRUE)#+sum((param-oldpar)^2*param^2)/2
cat("\n ",logLik, param)
- return(-logLik/sum(cond2))
+
+ #assign("oldpar",param,envir = my.envd1)
+
+ return(ret)
}
Added: pkg/yuima/R/DataPPR.R
===================================================================
--- pkg/yuima/R/DataPPR.R (rev 0)
+++ pkg/yuima/R/DataPPR.R 2017-11-29 17:49:46 UTC (rev 633)
@@ -0,0 +1,81 @@
+# The following example is necessary for the structure of
+# the PPR Data
+# library(yuima)
+# Terminal <- 10
+# samp <- setSampling(T=Terminal,n=1000)
+#
+# # Ex 1. (Simple homogeneous Poisson process)
+# mod1 <- setPoisson(intensity="lambda", df=list("dconst(z,1)"))
+# set.seed(123)
+# y1 <- simulate(mod1, true.par=list(lambda=1),sampling=samp)
+# y1 at data@original.data
+# simprvKern->yuimaPPR
+get.counting.data<-function(yuimaPPR,type="zoo"){
+ count <- yuimaPPR at Ppr@counting.var
+ dimCount <- length(count)
+ Time_Arrivals_Grid <- index(yuimaPPR at data@zoo.data[[1]])
+ if(dimCount==1){
+ cond <- yuimaPPR at model@solve.variable %in% count
+ Count.Var <- yuimaPPR at data@original.data[,cond]
+ }
+
+ Time_Arrivals <- t(Time_Arrivals_Grid[1])
+ colnames(Time_Arrivals) <- "Time"
+ Obser <- t(yuimaPPR at data@original.data[1,])
+ colnames(Obser) <- yuimaPPR at model@solve.variable
+ cond <- diff(Count.Var)==1
+ Time_Arrivals <- rbind(Time_Arrivals,as.matrix(Time_Arrivals_Grid[-1][cond]))
+ if(is(yuimaPPR at data@original.data,"matrix")){
+ Obser <- rbind(Obser,yuimaPPR at data@original.data[-1,][cond,])
+ }else{
+ Obser <- rbind(Obser,as.matrix(yuimaPPR at data@original.data[-1][cond]))
+ }
+
+ Data <- zoo(x = Obser,order.by = Time_Arrivals)
+ plot(Data)
+ if(type=="zoo"){
+ Data <- Data
+ return(Data)
+ }
+ if(type=="yuima.Ppr"){
+ yuimaPPR at data@original.data <- Data
+ return(yuimaPPR)
+ }
+ if(type=="matrix"){
+ Data <- cbind(Time_Arrivals, Obser)
+ return(Data)
+ }
+ yuima.stop("type is not supported.")
+}
+
+DataPpr <- function(CountVar, yuimaPPR, samp){
+ if(!is(yuimaPPR,"yuima.Ppr")){
+ yuima.stop("...")
+ }
+ if(!is(samp,"yuima.sampling")){
+ yuima.stop("...")
+ }
+ if(!is(CountVar,"zoo")){
+ yuima.stop("...")
+ }
+ names <- colnames(CountVar)
+ if(all(names %in% yuimaPPR at model@solve.variable) && all(yuimaPPR at model@solve.variable %in% names)){
+ TimeArr <- index(CountVar)
+ grid <- samp at grid[[1]]
+ dimData <- dim(CountVar)[2]
+ DataOr <- NULL
+ for(i in c(1:dimData)){
+ prv <- approx(x=TimeArr, y = CountVar[,i], xout = grid, method = "constant")$y
+ DataOr <- cbind(DataOr,prv)
+ }
+ DataFin <- zoo(DataOr, order.by = grid)
+ colnames(DataFin) <- names
+ myData <- setData(DataFin)
+ yuimaPPR at data <- myData
+ yuimaPPR at sampling <- samp
+ return(yuimaPPR)
+ }else{
+ dummy <- paste0("The labels of variables should be ", paste0(yuimaPPR at model@solve.variable,collapse= ", "),collapse = " ")
+ yuima.stop(dummy)
+ }
+}
\ No newline at end of file
Added: pkg/yuima/man/DataPpr.Rd
===================================================================
--- pkg/yuima/man/DataPpr.Rd (rev 0)
+++ pkg/yuima/man/DataPpr.Rd 2017-11-29 17:49:46 UTC (rev 633)
@@ -0,0 +1,93 @@
+\name{DataPpr}
+\alias{DataPpr}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{From \code{zoo} data to \code{yuima.Ppr}.}
+\description{The function converts an object of class \code{zoo} to an object of class \code{\link{yuima.Ppr}}.}
+\usage{
+DataPpr(CountVar, yuimaPPR, samp)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{CountVar}{An object of class \code{zoo} that contains counting variables and covariates. \code{index(CountVar)} returns the arrival times.}
+ \item{yuimaPPR}{An object of class \code{\link{yuima.Ppr}} that contains a mathematical description of the point process regression model assumed to be the generator of the observed data.}
+ \item{samp}{An object of class \code{yuima.sampling}.}
+}
+\value{The function returns an object of class \code{\link{yuima.Ppr}} where the slot \code{model} contains the Point process described in \code{yuimaPPR at model}, the slot \code{data} contains the counting variables and the covariates observed on the grid in \code{samp}.}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\examples{
+\dontrun{
+# In this example we generate a dataset contains the Counting Variable N
+# and the Covariate X.
+# The covariate X is an OU driven by a Gamma process.
+
+# Values of parameters.
+mu <- 2
+alpha <- 4
+beta <-5
+
+# Law definition
+my.rKern <- function(n,t){
+ res0 <- t(t(rgamma(n, 0.1*t)))
+ res1 <- t(t(rep(1,n)))
+ res <- cbind(res0,res1)
+ return(res)
+}
+
+Law.PprKern <- setLaw(rng = my.rKern)
+
+# Point Process definition
+modKern <- setModel(drift = c("0.4*(0.1-X)","0"),
+ diffusion = c("0","0"),
+ jump.coeff = matrix(c("1","0","0","1"),2,2),
+ measure = list(df = Law.PprKern),
+ measure.type = c("code","code"),
+ solve.variable = c("X","N"),
+ xinit=c("0.25","0"))
+
+gFun <- "exp(mu*log(1+X))"
+#
+Kernel <- "alpha*exp(-beta*(t-s))"
+
+prvKern <- setPpr(yuima = modKern,
+ counting.var="N", gFun=gFun,
+ Kernel = as.matrix(Kernel),
+ lambda.var = "lambda", var.dx = "N",
+ lower.var="0", upper.var = "t")
+
+# Simulation
+
+Term<-200
+seed<-1
+n<-20000
+
+true.parKern <- list(mu=mu, alpha=alpha, beta=beta)
+
+
+set.seed(seed)
+# set.seed(1)
+
+time.simKern <-system.time(
+ simprvKern <- simulate(object = prvKern, true.parameter = true.parKern,
+ sampling = setSampling(Terminal =Term, n=n))
+)
+
+
+plot(simprvKern,main ="Counting Process with covariates" ,cex.main=0.9)
+
+# Using the function get.counting.data we extract from an object of class
+# yuima.Ppr the counting process N and the covariate X at the arrival times.
+
+CountVar <- get.counting.data(simprvKern)
+
+plot(CountVar)
+
+# We convert the zoo object in the yuima.Ppr object.
+
+sim2 <- DataPpr(CountVar, yuimaPPR=simprvKern, samp=simprvKern at sampling)
+
+}
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
Added: pkg/yuima/man/get.counting.data.RD
===================================================================
--- pkg/yuima/man/get.counting.data.RD (rev 0)
+++ pkg/yuima/man/get.counting.data.RD 2017-11-29 17:49:46 UTC (rev 633)
@@ -0,0 +1,154 @@
+\name{get.counting.data}
+\alias{get.counting.data}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+Extract arrival times from an object of class \code{yuima.Ppr}
+}
+\description{This function extracts arrival times from an object of class \code{\link{yuima.Ppr}}.}
+\usage{
+get.counting.data(yuimaPPR,type="zoo")
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{yuimaPPR}{An object of class \code{\link{yuima.Ppr}}.}
+ \item{type}{By default \code{type="zoo"} the function returns an object of class \code{zoo}. Other values are \code{yuima.Ppr} and \code{matrix}.}
+}
+%\details{
+%% ~~ If necessary, more details than the description above ~~
+%}
+\value{By default the function returns an object of class zoo. The arrival times can be extracted by applying the method \code{index} to the output}
+%\references{
+%% ~put references to the literature/web site here ~
+%}
+%\author{
+%% ~~who you are~~
+%}
+%\note{
+%% ~~further notes~~
+%}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+%\seealso{
+%% ~~objects to See Also as \code{\link{help}}, ~~~
+%}
+\examples{
+\dontrun{
+##################
+# Hawkes Process #
+##################
+
+# Values of parameters.
+mu <- 2
+alpha <- 4
+beta <-5
+
+# Law definition
+
+my.rHawkes <- function(n){
+ res <- t(t(rep(1,n)))
+ return(res)
+}
+
+Law.Hawkes <- setLaw(rng = my.rHawkes)
+
+# Point Process Definition
+
+gFun <- "mu"
+Kernel <- "alpha*exp(-beta*(t-s))"
+
+modHawkes <- setModel(drift = c("0"), diffusion = matrix("0",1,1),
+ jump.coeff = matrix(c("1"),1,1), measure = list(df = Law.Hawkes),
+ measure.type = "code", solve.variable = c("N"),
+ xinit=c("0"))
+
+prvHawkes <- setPpr(yuima = modHawkes, counting.var="N", gFun=gFun,
+ Kernel = as.matrix(Kernel), lambda.var = "lambda",
+ var.dx = "N", lower.var="0", upper.var = "t")
+
+true.par <- list(mu=mu, alpha=alpha, beta=beta)
+
+set.seed(1)
+
+Term<-70
+n<-7000
+
+# Simulation trajectory
+
+time.Hawkes <-system.time(
+ simHawkes <- simulate(object = prvHawkes, true.parameter = true.par,
+ sampling = setSampling(Terminal =Term, n=n))
+)
+
+# Arrival times of the Counting process.
+
+DataHawkes <- get.counting.data(simHawkes)
+TimeArr <- index(DataHawkes)
+
+##################################
+# Point Process Regression Model #
+##################################
+
+# Values of parameters.
+mu <- 2
+alpha <- 4
+beta <-5
+
+# Law definition
+my.rKern <- function(n,t){
+ res0 <- t(t(rgamma(n, 0.1*t)))
+ res1 <- t(t(rep(1,n)))
+ res <- cbind(res0,res1)
+ return(res)
+}
+
+Law.PprKern <- setLaw(rng = my.rKern)
+
+# Point Process definition
+modKern <- setModel(drift = c("0.4*(0.1-X)","0"),
+ diffusion = c("0","0"),
+ jump.coeff = matrix(c("1","0","0","1"),2,2),
+ measure = list(df = Law.PprKern),
+ measure.type = c("code","code"),
+ solve.variable = c("X","N"),
+ xinit=c("0.25","0"))
+
+gFun <- "exp(mu*log(1+X))"
+#
+Kernel <- "alpha*exp(-beta*(t-s))"
+
+prvKern <- setPpr(yuima = modKern,
+ counting.var="N", gFun=gFun,
+ Kernel = as.matrix(Kernel),
+ lambda.var = "lambda", var.dx = "N",
+ lower.var="0", upper.var = "t")
+
+# Simulation
+
+Term<-100
+seed<-1
+n<-10000
+
+true.parKern <- list(mu=mu, alpha=alpha, beta=beta)
+
+
+set.seed(seed)
+# set.seed(1)
+
+time.simKern <-system.time(
+ simprvKern <- simulate(object = prvKern, true.parameter = true.parKern,
+ sampling = setSampling(Terminal =Term, n=n))
+)
+
+
+plot(simprvKern,main ="Counting Process with covariates" ,cex.main=0.9)
+
+# Arrival Times
+CountVar <- get.counting.data(simprvKern)
+TimeArr <- index(CountVar)
+
+
+}
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
More information about the Yuima-commits
mailing list