You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
where $GP(\mu(s), K(s,s'))$ denotes a Gaussian process with mean function $\mu(s)$ and covariance kernel $K(s,s')$ for a location $s \in X$. This is one of the simplest models of point patterns of $n$ events recorded as locations $s_1,...,s_n$ in an arbitrary metric space. In conjunction with a Bayesian analysis, this model can be used to answering questions of interest such as:
34
+
where $GP(\mu(s), K(s,s'))$ denotes a Gaussian process with mean function $\mu(s)$ and covariance kernel $K(s,s')$ for a location $s \in X$. This is one of the simplest models of point patterns of $n$ events recorded as locations $s_1,...,s_n$ in an arbitrary metric space. In conjunction with a Bayesian analysis, this model can be used to answer questions of interest such as:
35
35
* Does an observed point pattern imply a statistically significant shift in spatial intensities?
36
36
* What would randomly sampled patterns with the same statistical properties look like?
37
37
* Is there a statistical correlation between the *frequency* and *magnitude* of point events?
38
38
39
-
In this notebook, we'll use a grid-based approximation to the full LGCP with PyMC to fit a model and analyze its posterior summaries. We will also explore the usage of a marked Poisson process, an extension of this model to account for the distribution of *marks* associated with each data point.
39
+
In this notebook, we'll use a grid-based approximation to the full LGCP with PyMC to fit a model and analyze its posterior summaries. We will also explore the use of a marked Poisson process, an extension of this model to account for the distribution of *marks* associated with each data point.
40
40
41
41
+++
42
42
43
43
## Data
44
44
45
45
+++
46
46
47
-
Our observational data concerns 231 sea anemones whose sizes and locations on the French coast were recorded. This data was taken from the [`spatstat` spatial modeling package in R](https://github.com/spatstat/spatstat) which is designed to address models like the LGCP and its subsequent refinements. The original source of this data is the textbook *Spatial data analysis by example* by Upton and Fingleton (1985) and a longer description of the data can be found there.
47
+
Our observational data concerns 231 sea anemones whose sizes and locations on the French coast were recorded. These data were taken from the [`spatstat` spatial modeling package in R](https://github.com/spatstat/spatstat), which is designed to address models like the LGCP and its subsequent refinements. The original source of this data is the textbook *Spatial data analysis by example* by Upton and Fingleton {cite}`upton1985spatial` and a longer description of the data can be found there.
48
48
49
49
```{code-cell} ipython3
50
+
import os
50
51
import warnings
51
52
52
53
from itertools import product
@@ -57,19 +58,19 @@ import numpy as np
57
58
import pandas as pd
58
59
import pymc as pm
59
60
60
-
from matplotlib import MatplotlibDeprecationWarning
@@ -91,7 +92,7 @@ The 'marks' column indicates the size of each anemone. If we were to model both
91
92
92
93
+++
93
94
94
-
While there are multiple ways to conduct inference, perhaps the simplest way is to slice up our domain $X$ into many small pieces $A_1, A_2,...,A_M$ and fix the intensity field to be constant within each subset. Then, we will treat the number of points within each $A_j$ as a Poisson random variable such that $Y_j \sim Poisson(\lambda_j)$. and we also consider the $\log{\lambda_1}...,\log{\lambda_M}$ variables as a single draw from a Gaussian process.
95
+
While there are multiple ways to conduct inference, perhaps the simplest way is to slice up our domain $X$ into many small pieces $A_1, A_2,...,A_M$ and fix the intensity field to be constant within each subset. Then, we will treat the number of points within each $A_j$ as a Poisson random variable such that $Y_j \sim Poisson(\lambda_j)$, and we also consider the $\log{\lambda_1}...,\log{\lambda_M}$ variables as a single draw from a Gaussian process.
95
96
96
97
+++
97
98
@@ -103,7 +104,6 @@ xy = data[["x", "y"]].values
103
104
# Jitter the data slightly so that none of the points fall exactly
@@ -154,7 +154,9 @@ We can see that all of the counts are fairly low and range from zero to five. Wi
154
154
Our first step is to place prior distributions over the high-level parameters for the Gaussian process. This includes the length scale $\rho$ for the covariance function and a constant mean $\mu$ for the GP.
With the model fully specified, we can start sampling from the posterior using the default NUTS sampler. I'll also tweak the target acceptance rate to reduce the number of divergences.
Posterior inference on the length_scale parameter is useful for understanding whether or not there are long-range correlations in the data. We can also examine the mean of the log-intensity field, but since it is on the log scale it is hard to directly interpret.
190
192
191
193
```{code-cell} ipython3
192
-
az.summary(trace, var_names=["mu", "rho"])
194
+
az.summary(idata, var_names=["mu", "rho"])
193
195
```
194
196
195
197
We are also interested in looking at the value of the intensity field at a large number of new points in space. We can accommodate this within our model by including a new random variable for the latent Gaussian process evaluated at a denser set of points. Using `sample_posterior_predictive`, we generate posterior predictions on new data points contained in the variable `intensity_new`.
Let's take a look at a few realizations of $\lambda(s)$. Since the samples are on the log scale, we'll need to exponentiate them to obtain the spatial intensity field of our 2D Poisson process. In the plot below, the observed point pattern is overlaid.
0 commit comments