ippm {spatstat} | R Documentation |
Experimental extension to ppm
.
Find optimal values of the irregular trend parameters in a
point process model using Fisher scoring algorithm.
ippm(..., iScore=NULL, start=list(), covfunargs=start, maxiter=20, tol=1e-4, progress=TRUE, stepfactor=1, dbug=FALSE)
... |
Arguments passed to |
iScore |
A named list of R functions that compute the partial derivatives
of |
start |
Named list containing initial values of the irregular parameters over which to optimise. |
covfunargs |
Argument passed to |
maxiter |
Integer. Maximum number of iterations of Fisher scoring algorithm. |
tol |
Numeric value or vector. The algorithm stops when the
difference between two successive estimates of the irregular
parameter is less than |
progress |
Logical. Whether to print progress reports. |
stepfactor |
Numeric value between 0 and 1 indicating that the change in the parameter between successive iterations is only a specified fraction of the step computed by the Newton-Raphson algorithm. |
dbug |
Logical. Whether to print debugging output. |
This function is an experimental extension to the
point process model fitting command ppm
.
The extension allows the trend of the model to include irregular parameters,
which will be maximised by a Fisher scoring method.
For the sake of explanation, consider a Poisson point process with intensity function lambda(u) at location u. Assume that
lambda(u) = \exp(alpha + beta * Z(u)) * f(u, gamma)
where alpha,beta,gamma are parameters to be estimated, Z(u) is a spatial covariate function, and f is some known function. Then the parameters alpha,beta are called regular because they appear in a loglinear form; the parameter gamma is called irregular.
To fit this model using ippm
, we specify the
intensity using the trend
formula
in the same way as usual for ppm
.
The trend formula is a representation of the log intensity.
In the above example the log intensity is
log(lambda(u)) = alpha + beta * Z(u) + log(f(u, gamma))
So the model above would be encoded with the trend formula
~Z + offset(log(f))
. Note that the irregular part of the model
is an offset term, which means that it is included in the log trend
as it is, without being multiplied by another regular parameter.
To perform Fisher scoring we also need the derivative
of log(f(u,gamma)) with
respect to gamma. We call this the
irregular score. The user must write an R function
that computes the irregular score for any value of
gamma at any location (x,y)
.
Thus, to code such a problem,
The argument trend
should define the
log intensity, with the irregular part as an offset;
The argument start
should be a list
containing initial values of each of the irregular parameters;
The argument iScore
must be a list (with one entry
for each entry of start
) of functions
with arguments x,y,...
, that evaluate the partial derivatives
of log(f(u,gamma)) with
respect to each irregular parameter.
The coded example below illustrates the model with two irregular parameters gamma,delta and irregular term
f((x,y), (gamma, delta)) = 1 + \exp(gamma - delta * x^3)
Arguments ...
passed to ppm
may
also include interaction
. In this case the model is not
a Poisson point process but a more general Gibbs point process;
the trend formula trend
determines the first-order trend
of the model (the first order component of the conditional intensity),
not the intensity.
A fitted point process model (object of class "ppm"
).
Adrian Baddeley Adrian.Baddeley@csiro.au http://www.maths.uwa.edu.au/~adrian/ and Rolf Turner r.turner@auckland.ac.nz
nd <- 32 gamma0 <- 3 delta0 <- 5 POW <- 3 # Terms in intensity Z <- function(x,y) { -2*y } f <- function(x,y,gamma,delta) { 1 + exp(gamma - delta * x^POW) } # True intensity lamb <- function(x,y,gamma,delta) { 200 * exp(Z(x,y)) * f(x,y,gamma,delta) } # Simulate realisation lmax <- max(lamb(0,0,gamma0,delta0), lamb(1,1,gamma0,delta0)) set.seed(42) X <- rpoispp(lamb, lmax=lmax, win=owin(), gamma=gamma0, delta=delta0) # Partial derivatives of log f DlogfDgamma <- function(x,y, gamma, delta) { topbit <- exp(gamma - delta * x^POW) topbit/(1 + topbit) } DlogfDdelta <- function(x,y, gamma, delta) { topbit <- exp(gamma - delta * x^POW) - (x^POW) * topbit/(1 + topbit) } # irregular score Dlogf <- list(gamma=DlogfDgamma, delta=DlogfDdelta) # fit model ippm(X, ~Z + offset(log(f)), covariates=list(Z=Z, f=f), iScore=Dlogf, start=list(gamma=1, delta=1), tol=0.01, nd=nd)