\newcommand{\x}{\mathbf{x}}
This doc gives a short example of how to simulate shot-noise and product-noise Cox patterns.
We will use spatstat
to simulate the generator pattern and for some other tools, so load it as well.
library(rcox) library(spatstat)
The first step is to generate the generator pattern. Let's create two of them for two different outcomes:
set.seed(1) x1 <- runifpoint(50) x2 <- rMaternII(120, 0.08) generators <- listof(a=x1, b=x2) plot(generators)
Given a generator pattern $\x$, sample a pattern $\mathbf{y}$ from $Poisson( F(.;\x))$, where the intensity $F$ is given by one of the two ways:
sum
-type: For $u\in W$ where $W$ is the window,
[
F(u; \x) = \alpha \sum_{x\in\x} k_\sigma(u-x)
]
where $\alpha>0$ adjusts for intensity and $k_\sigma$ is a kernel (prob. density in 2D) with parameter $\sigma$.
This is called shot-noise-field. Some details of this model:
product
-type: For $u\in W$
[
F(u;\x) = e^{\alpha_1}\prod_{x\in \x}[1 + \alpha_2 k_\sigma(u-x)/k_\sigma(0)]
]
where $\alpha_1$ again adjusts the intensity and $\alpha_2\in[-1,\infty)$ is an interaction parameter. This is called product-shot-noise, following Jalilian et al. 2015. Some details when $\x$ is from $Poisson(\lambda_x)$:
We set up the models like this:
lam1 <- lambda(x1, kernel = "gauss", sigma = 0.05, alpha = 3, type = "sum") lam2 <- lambda(x2, kernel = "step", sigma = 0.1, alpha = c(log(100), -0.5), type = "prod")
The documentation lists the options.
We can "evaluate" the $F$ on a grid to see how the surface looks like.
field1 <- coxintensity2matrix(lam1, W = square(1)) field2 <- coxintensity2matrix(lam2, bbox = cbind(0:1, 0:1), nx = 128)
Note how the bounding box can be given in two ways.
par(mfrow=c(1,2), mar=c(1,1,1,1)) plot(field1, axes =F) points(x1) plot(field2, axes =F, scale=log) points(x2)
Red means low intensity (not preferred), white means high intensity (preferred).
The simulation uses either a thinning approach (free point count) or Metropolis-Hastings (fixed point count). The calls:
s1 <- rcox(lambda = lam1, bbox = cbind(0:1, 0:1)) s2 <- rcox(lambda = lam2, n = 100, bbox = cbind(0:1, 0:1), iter = 1e4, verb=T) # turn to ppp-objects y1 <- cox2ppp(s1) y2 <- cox2ppp(s2)
par(mfrow=c(1,2), mar=c(1,1,1,1)) plot(y1) points(x1, pch=3, col = 2) plot(y2) points(x2, pch=3, col = 2)
There we have it.
EOF
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.