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.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.density.fd
# 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 <- 10
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)
Run the code above in your browser using DataLab