damllRH | R Documentation |
Calculates an apprximation to the minus loglikelihood of a
RHawkes model with given immigration hazard function μ,
offspring birth time density function h and branching ratio
η relative to event times tms
on interval [0,cens].
damllRH(tms, cens, par, q=0.999, qe=0.999, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) }, keepB=FALSE, H.inv=function(x,p)qexp(x,rate=1/p) )
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
cens |
A numericl scalar. The censoring time. |
par |
A numeric vector containing the parameters of the model, in order of the immigration parameters, in μ(.), offspring distribution parameters, in h(.), and lastly the branching ratio η(.). |
q |
A numeric scalar in (0,1] and close to 1, which controls how far we look back when truncating the distribution of the most recent immigrant. |
qe |
A numeric scalar in (0,1] and close to 1, which controls how to truncation is used in the offspring birth time distribution. |
h.fn |
A (vectorized) function. The offspring birth time density function. |
mu.fn |
A (vectorized) function. The immigrant waiting time hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
keepB |
A boolean scalar, indicating whether the looking back
values |
H.inv |
A (vectorized) function, giving the inverse function of the integral of the excitation. |
A scalar giving the value of the (approximate) negative
log-likelihood, when keepB
is FALSE (the default); A list with
components mll
, whhich contains the value of the negative
log-likelihood, Bs
, which gives the look-back order of the
truncation of the distribution of the last immigrant, and Bes
,
which gives the look-forward order in determining how far into the
future the excitation effect is allowed to last.
Feng Chen <feng.chen@unsw.edu.au>
## Not run: ## earthquake times over 96 years data(quake); tms <- sort(quake$time); # add some random noise to the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) ## calculate the minus loglikelihood of an RHawkes with some parameters ## the default hazard function and density functions are Weibull and ## exponential respectively mllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5)) damllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5),q=1,qe=1) ## calculate the MLE for the parameter assuming known parametric forms ## of the immigrant hazard function and offspring density functions. system.time(est <- optim(c(0.5, 20, 1000, 0.5), mllRH, tms = tms, cens = 96*365.25, mu.fn=function(x,p)p[1]/p[2]*(x/p[2])^(p[1]-1), Mu.fn=function(x,p)(x/p[2])^p[1], control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) system.time(est1 <- optim(c(0.5, 20, 1000, 0.5), function(p){ if(any(p<0)||p[4]<0||p[4]>=1) return(Inf); damllRH(tms = tms, cens = 96*365.25, mu.fn=function(x,p)p[1]/p[2]*(x/p[2])^(p[1]-1), Mu.fn=function(x,p)(x/p[2])^p[1], par=p,q=0.999999,qe=0.999999) }, control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) ## point estimate by MLE est$par est1$par ## standard error estimates: diag(solve(est$hessian))^0.5 diag(solve(est1$hessian))^0.5 ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.