\documentclass[11pt]{article}
\usepackage{a4}
\title{75563 Spatial Statistics, \\ Assignment 2:Point processes}
\author{}
\date{}
\setlength{\parindent}{0mm}
\setlength{\parskip}{1.5mm}

\begin{document}

\maketitle

\section*{Exercise 1: Inhomogeneous Poisson process}
We will study an inhomogeneous Poisson process at the unit square.
The process has the intensity function: 

$$\lambda_{\theta}(x,y)=\mu \exp(-\alpha x -\beta y)$$
where $\theta=(\mu,\alpha,\beta)$. 

{\bf a.} The number of points in the unit square is a random
variable. Compute its mean and variance. 

{\bf b.} Simulate the inhomogeneous Poisson process
conditioned to the number of points in the unit square.
Let the intensity function be as given with 
$\theta=(200,2.0,1.0)$.

Vary the  number of points in the unit square. (f.ex: $n=(20,50,100)$).


{\bf Splus  : } Use \verb+par(pty=``s'')+ to plot at a square.

{\bf c.}  Expand the algorithm in b, to simulate unconditional to 
the number of points and do the simulation.


{\bf d.}  An alternative algorithm to sample the 
inhomogeneous Poisson process,is based on p - thinning.

To simulate at the unit square do the following.

\begin{enumerate}
\item  put $$\lambda^*=\mathop{\bf max}_{ (x,y) \in [0,1]^2 }\lambda_{\theta}(x,y) $$
\item Simulate a homogeneous point process with intensity $\lambda^*$.
\item Keep a simulated point  $(x_n,y_n)$ with probability 
$$p=\frac{\lambda_{\theta}(x_n,y_n)}{\lambda^*}$$
\end{enumerate}
Implement this algorithm as well, and do the simulation.
What is the main difference of the two algorithms.
Generally speaking, what is the pros and the cons of this 
algorithm compared to the algorithm you used in c.
                 
{\bf e.} Identify the equations used to find the maximum likelihood 
estimate of $\theta=(\mu,\alpha,\beta)$  in the model above.
The likelihood is given by:
 
$$L(x_1,y_1\ldots,x_n,y_n;\theta)=\lambda_{\theta}(x_1,y_1)\ldots
\lambda_{\theta}(x_n,y_n)\exp[-\int_0^1 \int_0^1 \lambda_{\theta}(x,y)dxdy],$$
 at the unit square.

{\bf f.} Choose a point pattern from the unconditional simulation (c
or d). Determine the maximum likelihood estimates of 
$\theta=(\mu,\alpha,\beta)$ based on that data set.
Plot the estimated intensity with the function \verb+image+. 


{\bf g.} Copy the files intest.s and integral.o. 
from \verb+/local/kurs/romstat+ . 
This is code for a non-parametric estimate of the intensity,
based on the kernel estimation, (ref. Cressie 1993, p 600)
How does the function intest.s work ?
Use it to estimate the intensity for the same set of data  you used
for the maximum likelihood estimation. Try different 
smoothing parameters. Plot the intensity, and compare the two estimates.
(Use the same color range for the two  intensities if you plot with
image.)

\section*{Exercise 2: Neyman-Scott process}

In this exercise we will study Neyman-Scott processes where the 
number of "child's" are Poisson distributed.
The distribution relative to the mother, you can choose as 
a symmetric normal distribution (Thomas process), or the uniform 
distribution for a disk with radius $R$ (Matern cluster process).
(ref: Neyman-Scott processes;  Stoyan \& Stoyan (1994)  ch.16.2)

{\bf a.} Simulate a Neyman-Scott process at the unit square
unconditional to the number of points. What would you do to 
avoid boundary effects ? Vary the parameters of the model. 
How does it affect the point patterns.


{\bf b.}A Neyman-Scott process where the number of child's are 
Poisson distributed, is equivalent with a Cox process with 
intensity function:

$$\Lambda(r)=\omega\sum_{i=1}^{\infty}f(r-s_i) \; ,$$
where $\omega$ is expected number of child's for each mother,
$f$ is the distribution of the "children's" position relative to 
the mother and $s_i$ is the location of mother number $i$.

Loosely one may say: A Cox- process is a "inhomogeneous Poisson
distribution with stochastic intensity". Which features of 
$\Lambda(r)$ are  stochastic ?
How could You simulate $\Lambda(r)$ conditioned to the number of 
child's ? How would You progress in order to simulate a 
Neyman-Scott process conditional to the number of (child) points.
You do not have to do the simulation.


{\bf c.} Calculate the  $L$-function for some of the simulated 
point patterns (ref exercise 3), and compare with the theoretical
$L$- function, given under. (NB! Use the Splus commando \verb+ppregion()+ 
to initialize, the point process region.)

For the Thomas-process we have:
$$K(t)=\pi t^2+\rho^{-1}(1-\exp(-t^2/4\sigma^2))$$
Where $\rho$ is the intensity of the mother process, and 
 $\sigma^2$ is the variance of the "children's" distribution
relative to the mother. 

For the Matern cluster-process we have
$$K(t)=\pi t^2+\rho^{-1}\times \left\{
  \begin{array}{ll}2+\pi^{-1}[(8z^2-4)\arccos z - 2\arcsin z &  \\
                   +4z\sqrt{(1-z^2)^3}-6z\sqrt{1-z^2}] & (t\leq 2R) \\
                    1 & (t>2R)
  \end{array} \right.
$$
where  $z=t/2R$ and other notation as earlier.       

\section*{Exercise 3: Strauss process}
Load the "Spatial" library using the commando \verb+library(spatial)+,
this Library contains among other things a function for calculating 
the $L$ -function. ($L(t)=\sqrt{K(t)/\pi}$) 
with "envelopes", and a function for simulating the 
Strauss process, conditional to the number of points.
Relative to the binomial process, is the unnormalized density
of a Strauss process defined as  $c^{s(x)}$ where $c>0$ and $s(x)$ is the
number of point pairs $(\mbox{\bf x}_i,\mbox{\bf x}_j)$ 
where $\| \mbox{\bf  x}_i-\mbox{\bf  x}_j \|  < r$. 


{\bf a.} Simulate the  Strauss-process with $50$ points, $r=0.1$ and $c=0.01$,
$c=0.5$, $c=1.0$. Compare the point patterns, make sure the plots are
on the same scale.
Calculate the $L$-function of your point patterns, 
Plot it together with "envelopes" based on 99 simulations
from a Binomial process. What does the plots tell.(NB! use \verb+ppregion()+ )
What happens if $c > 1$ if 1) we condition to the number of 
points, 2) we do not condition to the number of points.


{\bf b.} Describe a Metropolis-Hastings algorithm to simulate 
from a Strauss-process conditioned to the number of points.
(You do not have to do the simulation.)

{\bf c.} Describe how to simulate an unconditional Strauss process
by the use of a Metropolis-Hastings algorithm.
(You do not have to do the simulation.)

\section*{Exercise 4: Analysis of  point patterns.}
In this exercise we will use three different point patterns.
"cells.dat","redwood.dat" and "pines2.dat".
Use the command \verb+ppinit+ from the "spatial" library to 
load the files.

{\bf a.} Plot the point patterns.
Does the the Poisson assumption apply for any of the 
point patterns? Use different tests to investigate.
(f.ex $L$-function with Poisson envelopes, dispersion index, others?)


{\bf b.} Implement an estimate of the $D$-function. $D$ is the 
distribution of the distance from a point to the nearest 
neighboring point.
Plot this for the three examples. How did you solve the boundary
problems ?
Useful Splus functions \verb+dist()+,  \verb+dist2full+ as
described in help entry for  \verb+dist()+, and
 \verb+apply( , ,min)+.

{\bf c.} Choose a point pattern that is not Poisson distributed 
and suggest a better model. Try to estimate the parameters of the
model you suggest.
Check the model fit by use of the $L$ and $D$ function.
Useful Splus commands. \verb+nlmin+ ,\verb+integrate+ ,\verb+pplik+
You do not have to use any of theese.

\section*{References:}
 
Cressie; N.; 1993: Statistics for Spatial Data, Wiley \& Sons. 

Stoyan; D. and Stoyan; H.; 1994: Fractals, Random Shapes and Point Fields; Wiley \& Sons. Kap 11,12,13,15,16


\end{document}


%%% Local Variables: 
%%% TeX-master: t
%%% End: 
