intensity.fd | Language Reference for FDA Library
|
intensity.fd
.
intensity.fd(x, WfdParobj, conv=0.0001, iterlim=20, dbglev=1)
mu(x) = exp[W(x)]
.
iternum+1
by 5 matrix containing the iteration
history.
The intensity function I(t) is almost the same thing as a
probability density function p(t) estimated by function
densify.fd
. The only difference is the absence of
the normalizing constant C that a density function requires
in order to have a unit integral.
The goal of the function is provide a smooth intensity function
estimate that approaches some target intensity by an amount that is
controlled by the linear differential operator
Lfdobj
and
the penalty parameter in argument
WfdPar
.
For example, if the first derivative of
W(t) is penalized heavily, this will force the function to
approach a constant, which in turn will force the estimated Poisson
process itself to be nearly homogeneous.
To plot the intensity function or to evaluate it,
evaluate
Wfdobj
, exponentiate the resulting vector.
# Generate 101 Poisson-distributed event times with # intensity or rate two events per unit time N <- 101 mu <- 2 # generate 101 uniform deviates uvec <- runif(rep(0,N)) # convert to 101 exponential waiting times wvec <- -log(1-uvec)/mu # accumulate to get event times tvec <- cumsum(wvec) tmax <- max(tvec) # set up an order 4 B-spline basis over [0,tmax] with # 21 equally spaced knots tbasis <- create.bspline.basis(c(0,tmax), 23) # set up a functional parameter object for W(t), # the log intensity function. The first derivative # is penalized in order to smooth toward a constant lambda <- 1e1 WfdParobj <- fdPar(tbasis, 1, lambda) # estimate the intensity function Wfdobj <- intensity.fd(tvec, WfdParobj)$Wfdobj # get intensity function values at 0 and event times events <- c(0,tvec) intenvec <- exp(eval.fd(events,Wfdobj)) # plot intensity function plot(events, intenvec, type="b") lines(c(0,tmax),c(mu,mu),lty=4)