\documentclass[12pt]{article}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage[norsk]{babel}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{a4wide}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{parskip}
\usepackage{verbatim} 

\usepackage{graphicx}
\usepackage[hang]{subfigure}

\newcommand{\hs}{\hspace*}
\newcommand{\mc}{\mathcal} 
\newcommand{\mr}{\mathrm}
\newcommand{\pr}{x^{\prime}}
\newcommand{\dpr}{x^{\prime \prime}}








\setlength{\parindent}{0cm}
\setlength{\parskip}{\baselineskip}

\title{Spatial Statistics \\ Assignment }
\author{Henning Omre}
\date{September 2004}

\begin{document}

\maketitle
\section{Part I: Continous Random Fields}
The R/Spatial Library is required for solving the problems. The command \texttt{library(spatial)} loads the library.

\subsection*{Problem 1: 1D Gaussian RF}

Consider a Gaussian RF $ \{R(x);x \in \mathcal{D} \subset \mathcal{R}^{1} \}; \mathcal{D} \in [0,50]$  with 

\begin{equation}
\begin{split}
E\{R(x)\}&=0 \\ 
Cov\{R(\pr ),R(\dpr)\}&= \sigma^2 \times \exp\{-\frac{|\pr-\dpr|^{\nu}}{\beta}\} \\
\end{split}
\end{equation}

where $\sigma^2 , \beta $ and $\nu $ are model parameters.

\textbf{(a)} Plot the spatial covariance function for different values of the parameters, start with $(\sigma^2,\beta,\nu)=(1,3,1)$. Describe the influence of each of the parameters on the Gaussian RF.

The spatial variogram function is defined as

\begin{equation}
\gamma (\pr-\dpr)=\frac{1}{2} Var\{R(\pr)-R(\dpr)\}
\end{equation}

Develop the expression for the variogram function for the RF specified above, and plot it for one choice of parameter values. Explain which features are denoted sill and range.

\textbf{Hint:} Use the functions  \texttt{expcov} and \texttt{gaucov}.

\textbf{(b)} Discretize $\mathcal{D}$ into $\mathcal{L}_{\mathcal{D}} : [0,1,\dots,50]$ and define the discretized Gaussian RF $\{R(x);x \in \mathcal{L}_{\mathcal{D}} \} $. Use the definition of Gaussian RF to simulate this discretized RF for different values of the model parameters. Comment on the results.

\textbf{Hint:} Use the function  \texttt{mvrnorm} in the \texttt{MASS} library.

\textbf{(c)} Use the model parameters $ (\sigma^2,\beta,\nu) = (4,3,1)$ and simulate one realization of $ \{R(x);x \in \mathcal{L}_{\mathcal{D}}\} $. Plot the realization and compute the area under the realization and above level 2.0.

Let the realizations in locations $x=(10,15,40) $ be considered as observations $o$.

Consider the discretized conditional Gaussian RF $\{[R(x)|o]; x \in \mathcal{L}_{\mathcal{D}}\} $ and generate 100 realizations of it. Plot 10 of these realizations, and comment.

In each location $x \in [0,1, \dots ,50]$ compute the average and empirical variance over the 100 realizations. Plot the results, and comment.

Compute the discretized expected function $\{E\{R(x)|o\} ; x\in \mathcal{L}_{\mathcal{D}} \}$ and variance function $\{Var\{R(x)|o\} ; x \in \mathcal{L}_{\mathcal{D}} \}$, plot the results, and comment.

Compare the results obtained above.


\subsection*{Problem 2: Kriging Prediction}

Consider the terrain elevation in an area $\mathcal{D} : [(0,0) \times (7,7)] $, denoted $\{r(x); x \in \mathcal{D} \subset \mathcal{R}^2 \} $. The terrain elevation is observed in 52 locations, $ r^d=(r(x_1),r(x_2), \dots ,r(x_{52})) $, which are available in the \texttt{MASS} library. Type \texttt{data(topo)} to obtain the data
 


\textbf{(a)} Let the terrain elevations be modelled as an isotropic continuous RF $\{R(x); x \in \mathcal{D} \}$ with

\begin{equation}
\begin{split}
R(x) &= \sum_{r+s \leq p} a_{rs} (x^1)^r (x^2)^s + U(x) \\
Cov\{U(\pr),U(\dpr)\} &= exp\{-0.7 |\pr-\dpr|\}
\end{split}
\end{equation}

where $x=(x^1,x^2) $ and  $ a_{rs} $ are unknown coefficients.

Based on the model assumptions and the available observations, $ r^d $ , compute the Universal Kriging (UK) prediction for a suitable discrtization over $ \mathcal{D} $. Use $ p=0,2,4,6 $ in the calculations. Compute also the associated UK prediction standard deviation.

Plot the results as contour plots, gray-scale plots and perspective plots with the observations marked. Comment on the results.

\textbf{ Hint:} Use the function \texttt{surf.gls} and \texttt{semat} for UK and the functions \texttt{contour}, \texttt{image} and \texttt{persp} for plotting.



\textbf{(b)} Alternatively one may use an ordinary regression approach, termed trend surface, with

\begin{equation}
\begin{split}
R(x)&= \sum_{r+s \leq p} a_{rs} (x^1)^r (x^2)^s  + U(x) \\
Cov\{U(\pr),U(\dpr)\}&= \mathrm{I}(\pr=\dpr)\\
\end{split}
\end{equation}

where $ x=(x^1,x^2) $, $\mathrm{I}(A)$ is equal to 1 whenever $A$ is true an d zero else, and $ a_{rs} $ are unknown coefficients. Note that the residuals do not have any spatial dependence.

Based on these model assumptions and the available observations, $ r^d $, compute the trend surface for a suitable discretization over $ \mathcal{D} $. Use $ p= 0,2,4,6 $ in the calculations.

Plot the results as in point a) and comment.

\textbf{ Hint:} Use the function \texttt{surf.ls} for trend surface analysis.




\textbf{(c)}Discuss the following:

\begin{itemize}
\item similarities and differences for UK and trend surface
\item strengths and weaknesses for each approach
\item strengths and weaknesses of various display alternatives 
\end{itemize}

\section*{Part II: Event RF}
\subsection*{Problem 3: Gibbs RF/Strauss model}



Consider an Event RF $ \{X_i ; i= 1, \dots ,n; \mathcal{D} \} $ given number of 
points equal  $n$. Let it be represented by a Gibbs pdf:



\begin{equation}
f(x|n) = const \times \exp \{- \sum_i^n \sum_j^n \phi(\tau_{ij}) \}
\end{equation}



where $x=(x_1, \dots ,x_n) , \tau_{ij} = |x_i - x_j| $ and $ \phi( \tau ) $ is 
the interaction function,



\begin{equation}
 \phi( \tau ) = \begin{cases} \alpha & \tau< r \\
0 & \mbox{else}\\
\end{cases}
\end{equation}



This model is termed the Strauss model.







\textbf{(a)} Use the MCMC-MH algorithm to generate a realization of this Strauss 
model on a unit square with $ n=50 , r=0.1 $ and $\alpha = 10$. Repeat the 
simulations for $ \alpha =1 $ and $\alpha = 0 $. Justify the convergence of the 
algorithm.



Plot the results and comment.



\textbf{(b)} Estimate the L-function for each of the realizations in point a).



Simulate 100 realizations of a Binominal RF on the unit square with $ n=50 $. 
Estimate the expected and 0.8 envelope of the L-function for this Binominal RF 
based on these realizations.



Plot the L-functions estimated above and comment on the results.





\textbf{(c)} Assume that the number of points is random according to a Poisson 
distribution with parameter $ \lambda $. Specify the expression for $ f(x,n) $.

Specify a simulation algorithm for the corresponding Event RF - the algorithm 
need not be implemented.

\section*{Part III: Markov Fields}
\subsection*{Problem 5: Markov RF}
Consider a binary Markov RF $ L : \{L_x; x \in \mathcal{L}_\mathcal{D} \} $ 
where $ \mathcal{L}_\mathcal{D} $ is a lattice over $ \mathcal{D} \subset 
\mathcal{R}^2 $ and $ L_x \in [-1,1] $. The pdf is:



\begin{equation}
Prob\{L=l\} = const \times \exp\{\sum_{x \sim y} [\beta_- I(l_x=l_y=-1) + \beta_+ 
I(l_x=l_y=1)] \}
\end{equation}



where $ x \sim y $ indicates sum over closest neighbours, $ I(A) $ is one 
whenever $ A $ is true and zero otherwise, and $ (\beta_-,\beta_+) $ are model 
parameters. Note that by setting $ \beta_-=\beta_+ $ one obtains the Ising 
model.





\textbf{(a)} Specify the neighborhood and clique systems of the model.



Simulate three realizations on lattice $[32 \times 32] $ by the MCMC-MH 
algorithm with model parameters $\beta_- = \beta_+ = 0.23 $. Repeat the 
simulations for three other choices of model parameters. Justify the convergence of the algorithm.

Plot the results and comment. 

\textbf{ Hint:} Use Matlab in the implementation of the simulation algorithm. 

\textbf{(b)} Specify the pseudolikelihood function for the parameters in the 
model.

Use maximum pseudolikelihood to estimate $ (\beta_-,\beta_+) $ from each of the 
realizations generated in point a). Comment on the results.



\subsection*{Problem 6: Bayesian image analysis}



The noisy image $o:\{o_x; x \in \mathcal{L}_{\mathcal{D}} \} $ with $ 
\mathcal{L}_{\mathcal{D}} : [ 115 \times 52] $ and $o_x \in \mathcal{R}^1 $ can 
be copied from the file  noiseimage.mat. Assume that the likelihood model is:



\begin{equation}
O:\{O_x = l_x + U_x ; x \in \mathcal{L}_{\mathcal{D}} \}
\end{equation}



where $ l_x \in \{-1,1\} $ and $ U : \{U_x; x \in \mathcal{L}_{\mathcal{D}} \} $ 
are iid $ Gauss_1 (u_x;0,1.5^2) $.



Assume further that a representative prior model for $ l: \{ l_x ; x \in 
\mathcal{L}_{\mathcal{D}} \} $ is a Markov RF as 



\begin{equation}
Prob \{ L=l \} = const \times \exp \{- \frac{1}{2} \sum_{x \sim y} I(l_x=l_y) \}
\end{equation}



where $ x \sim y $ indicates summation over closest neighbours and $ I(A) $ 
equals one whenever  $ A $ is true and zero else.





\textbf{(a)} Specify the posterior model for $ [L|o] $. 



Demonstrate that the posterior model is also a Markov RF with the same 
neighborhood system as the prior model.





\textbf{(b)} Use the MCMC-Gibbs algorithm to generate realizations from the 
posterior model. Use torus border conditions in the algorithm. Justify that the 
algorithm has converged. 



Plot some of the realizations and comment. 



\textbf{(c)} Define the Locationwise maximum posterior (LMAP) predictor of the 
image.



Estimate LMAP for the image above, plot it and comment. 






\end{document}