mloglik1b {IHSEP}

Minus loglikelihood of an IHSEP model

Description

Calculates the minus loglikelihood of an IHSEP model with given baseline inensity function \nu and excitation function g for event times jtms on interval [0,TT].

Usage

mloglik1b(jtms, TT = max(jtms), nu, g,
          Ig=function(x)sapply(x,function(y)integrate(g,0,y,
               rel.tol=1e-12,abs.tol=1e-12,subdivisions=1000)$value),
          Inu=function(x)sapply(x,function(y)integrate(nu,0,y)$value))

Arguments

jtms

A numeric vector, with values sorted in ascending order. Jump times to fit the inhomogeneous self-exciting point process model on.

TT

A scalar. The censoring time, or the terminal time for observation. Should be (slightly) greater than the maximum of jtms.

nu

A (vectorized) function. The baseline intensity function.

g

A (vectorized) function. The excitation function.

Ig

A (vectorized) function. Its value at t gives the integral of the excitation function from 0 to t.

Inu

A (vectorized) function. Its value at t gives the integral of the baseline intensity function \nu from 0 to t.

Details

This version of the mloglik function uses external C code to speedup the calculations. When given the analytical form of Inu or a quickly calculatable Inu, it should be (probably slightly) faster than mloglik1a. Otherwise it is the same as mloglik0 and mloglik1a.

Value

The value of the negative log-liklihood.

Author(s)

Feng Chen <[e-mail]>

See Also

mloglik0 and mloglik1a

Examples

## simulated data of an IHSEP on [0,1] with baseline intensity function
## nu(t)=200*(2+cos(2*pi*t)) and excitation function
## g(t)=8*exp(-16*t)
data(asep)

## get the birth times of all generations and sort in ascending order 
tms <- sort(unlist(asep))
## calculate the minus loglikelihood of an SEPP with the true parameters 
mloglik1b(tms,TT=1,nu=function(x)200*(2+cos(2*pi*x)),
          g=function(x)8*exp(-16*x),Ig=function(x)8/16*(1-exp(-16*x)))
## calculate the MLE for the parameter assuming known parametric forms
## of the baseline intensity and excitation functions  
## Not run: 
system.time(est <- optim(c(400,200,2*pi,8,16),
                         function(p){
                           mloglik1b(jtms=tms,TT=1,
                                     nu=function(x)p[1]+p[2]*cos(p[3]*x),
                                     g=function(x)p[4]*exp(-p[5]*x),
                                     Ig=function(x)p[4]/p[5]*(1-exp(-p[5]*x)))
                         },
                         hessian=TRUE,control=list(maxit=5000,trace=TRUE))
            )

## point estimate by MLE
est$par
## standard error estimates:
diag(solve(est$hessian))^0.5

## End(Not run)