In siliusmv/tma4300ex2:

\section{Problem 4: Repulsive event spatial variables}

a)

The Strauss model is a model used for creating a repulsion event RF. We specify the model on conditional form, given the number of events $k_\text{D}$ within some area $\text{D}$.

\begin{equation} \left[\mathbb{X}\text{D}^S | k\text{D} = k\right] \sim p(\bm{x}1, \bm{x}_2, ..., \bm{x}_k | k\text{D} = k) = const \times \prod_{i, j \in {1, 2, ..., k}} \text{exp} \left{-\phi (|\bm{x}_i - \bm{x}_j|) \right}. \end{equation}

$\phi(\tau_{ij}) \in \mathbb{R}\oplus;\ \tau{ij} = |\bm{x}i - \bm{x}_j|$ is a pairwise interaction function that declines for increasing values of $\tau{ij}$. In our implementation the interaction function takes the form

\begin{equation} \phi({\tau}) = \left { \begin{alignedat}{2} &\phi_0;&& 0 \leqslant \tau \leqslant \tau_0 \ &\phi_0 \text{exp} {-\phi_1[\tau - \tau_0]};\ \ && \tau > \tau_0 \end{alignedat}\right. \end{equation}

If $\phi_0$ takes a high value, the probability of finding an event within a distance of less than $\tau_0$ from another event will be very small. $\tau_0$ therefore works as a border, and we will often expect a minimal distance between events with the same magnitude as $\tau_0$. The interaction function decreases towards zero for larger values of $\tau$. The rate of the decrease is controlled by the parameter $\phi_1$.

To test whether the simulated Markov chain has reached convergence, one can generate and display several Markov chains with different initial values to see when the initial state has been forgotten. First, we make a guess of the parameters needed for simulating the biological cell data. Then we test the convergence of the algorithm using these values. The results of this can be seen in figure \@ref(fig:burn-in).

(ref:burn-in) Four different functions for evaluating convergence of the Markov chain used for simulating the Strauss model. One sample is drawn for every tenth iteration of the Markov chain.

############
## PRoblem 4

phi_0 <- 100
phi_1 <- 50
tau_0 <- 0.02
k <- length(cells_data$x)
num_iter <- 2000
area <- cells_data$area
jumps <- 10

args <- list(phi_0 = phi_0,
             phi_1 = phi_1,
             tau_0 = tau_0,
             k = k,
             area = area,
             acceptanceFunc = acceptanceRepulsive,
             num_iter = num_iter,
             jumps = jumps)


num_test <- 10

gg_1 <- testBurnIn(evalFunc = evaluateMeanMinDist,
                   args = args,
                   num_test = num_test,
                   title = "Mean min distance")

gg_2 <- testBurnIn(evalFunc = evaluateMeanMaxDist,
                   args = args,
                   num_test = num_test,
                   title = "Mean max distance")

gg_3 <- testBurnIn(evalFunc = evaluateMaxDist,
                   args = args,
                   num_test = num_test,
                   title = "Max distance")

gg_4 <- testBurnIn(evalFunc = evaluateMinDist,
                   args = args,
                   num_test = num_test,
                   title = "Min distance")

grid.arrange(gg_1, gg_2, gg_3, gg_4)

After studying figure \@ref(fig:burn-in) we assume that any Markov chain has reached convergence after 700 iterations for the given set of parameters. Now, it is possible to test the Strauss model. We simulate 100 realisations of the model, and construct 0.9-confidence intervals for the L-function. The confidence intervals are displayed along with the empirical L-function for the redwood tree data in figure \@ref(fig:repulsive-test). Due to the low computational cost, the 100 samples are generated from 100 independent Markov chains to ensure iid data.

(ref:repulsive-test) The left figure shows an example of a simulated Strauss model with parameters $\phi_0 = 100$, $\phi_1 = 50$ and $\tau_0 = 0.02$. The right figure displays the empirical L-function for the biological cell data along with 0.9-confidence intervals for the L-function. These are made from simulation of 100 Strauss event RF models with the same parameters.

# evaluation of MCMC algorithm

args$num_iter <- 700

data <- generateRepulsivePois(args)

plotProcess(data, title = "Biological cells")

testModel(data = cells_data,
           simulatePois = generateRepulsivePois,
           args = args,
           num_real = 100,
           title = "Biological cells")

We see that the Strauss model with the given parameters capture the repulsion of the cell data very well. For larger values of $t$, the confidence-interval obtains a lower value of L$(T)$ than that of the biological cell data. This is not the most important area of for simulation, though. Additionally, the relative size of the error is not too large. It can also be seen that the simulation in figure \@ref(fig:repulsive-test) is similar to the point process of biological cell data in figure \@ref(fig:processes). We can therefore conclude that the Strauss model is good enough for modeling our given data set.

siliusmv/tma4300ex2 documentation built on May 26, 2019, 4:32 a.m.