# Quadrat-based statistics in `pointpats`

This notebook demonstrates quadrat-based methods for analyzing planar point patterns using the current `pointpats.quadrat_statistics` implementation.

The workflow is:

1. Represent the point pattern as a NumPy array of coordinates with shape `(n, 2)`.

2. Construct a rectangular or hexagonal quadrat tessellation over the minimum bounding box (MBB) of the points.

3. Compute a Pearson chi-square statistic from per-quadrat counts.

4. Optionally, obtain a Monte Carlo p-value by simulating CSR realizations within the same study window, with reproducibility controlled by `rng`.

Key implementation facts (relevant for interpretation):

- The default chi-square statistic uses equal expected counts across included cells.

- A `window` may be supplied; if omitted, the MBB rectangle is used.

- Monte Carlo inference makes the reference distribution window-aware, but it does not change how expected counts are defined.

## 1. Quadrat statistics and CSR

Quadrat analysis partitions space into a set of non-overlapping cells and compares the observed counts per cell to what is expected under complete spatial randomness (CSR).

Let \(O_i\) be the count in cell (i), and let (k) be the number of cells included in the analysis. The implemented test statistic is the Pearson chi-square statistic with equal expected counts across included cells:

\[X^2 = \sum_{i=1}^k \frac{(O_i - \bar{O})^2}{\bar{O}}, \qquad \bar{O} = \frac{1}{k} \sum_{i=1}^k O_i.\]

Analytical inference compares \(X^2\) to a chi-square distribution with \((k-1)\) degrees of freedom. Simulation-based inference (Monte Carlo) instead compares the observed \(X^2\) to the sampling distribution obtained from CSR realizations generated in the same window.

## 2. Example: the juvenile delinquency point pattern

We will illustrate the quadrat workflow using the `juvenile.shp` example dataset. The dataset is read as a GeoDataFrame and then converted to a NumPy coordinate array `xy` with shape `(n, 2)`.

These coordinates are passed directly to `QStatistic`.

    [ ]:

    import libpysal as ps
    import numpy as np
    import geopandas as gpd
    %matplotlib inline
    import matplotlib.pyplot as plt
    import warnings
    warnings.filterwarnings("ignore")

Import `pointpats.quadrat_statistics` and compute quadrat statistics.

The main entry point is:

    QStatistic(points, shape="rectangle"|"hexagon", realizations=0, window=None,
               drop_nonintersecting=False, rng=None)

- `points` is always a NumPy array of shape `(n, 2)`.

- `shape` selects the tessellation type.

- `realizations > 0` enables Monte Carlo inference.

- `rng` controls reproducibility of simulations.

    [ ]:

    from pointpats.quadrat_statistics import QStatistic

Open the point shapefile and extract coordinates.

Here we use `GeoDataFrame.get_coordinates()` to obtain an `(n, 2)` NumPy array, which is the required input type for `QStatistic`.

    [ ]:

    juv = gpd.read_file((ps.examples.get_path("juvenile.shp")))

    [ ]:

    juv.shape

    [ ]:

    xy = juv.get_coordinates()

    [ ]:

    xy

    [ ]:

    resRect = QStatistic(xy)

    [ ]:

    resRect.chi2_pvalue

    [ ]:

    resRect.plot()

    [ ]:

    resHex = QStatistic(xy, shape='hexagon')

    [ ]:

    resHex.plot()

    [ ]:

    resHex.chi2_pvalue

## 3. Simulation-based inference (Monte Carlo)

When `realizations > 0`, CSR realizations are generated within the same window and the chi-square statistic is recomputed for each realization. The Monte Carlo p-value uses the standard +1 correction:

\[p = \frac{\#\{X^2_{sim} \ge X^2_{obs}\} + 1}{R + 1},\]

where \(R\) is the number of realizations.

### Reproducibility

All simulation-based results can be made reproducible by supplying `rng`:

- `None` (default): create a new `numpy.random.default_rng()`

- integer seed

- `numpy.random.Generator`

- legacy `numpy.random.RandomState` (wrapped)

The same generator is passed through the CSR simulation code, ensuring repeatable Monte Carlo p-values.

    [ ]:

    from numpy.random import default_rng

    [ ]:

    rng = default_rng(42)

    [ ]:

    resHex_sim = QStatistic(xy, shape='hexagon', realizations=99, rng=rng)

    [ ]:

    resHex_sim.chi2_r_pvalue

    [ ]:

    resHex_sim.chi2_realizations

    [ ]:

    resHex_sim.chi2

## 4. Irregular windows

Quadrat grids are constructed over the MBB of the points. If the study window is irregular (e.g., a polygon), some grid cells may lie entirely outside the window.

- If such outside cells are included, they contribute guaranteed zeros and can distort the test.

### How Monte Carlo interacts with this

Comparing the observed chi-square statistic to simulated chi-square values is a valid Monte Carlo test for the statistic being used, because the same grid, window, and inclusion rule are applied to both observed and simulated patterns.

Monte Carlo inference makes the reference distribution window-aware, but it does not convert the test into an area-proportional quadrat test when the window is irregular.

    [ ]:

    from shapely.geometry import Polygon

    [ ]:

    window = Polygon([ [0, 0], [1, 0], [0, 1] ])

    [ ]:

    from pointpats.random import poisson

    [ ]:

    pp = poisson(window, 100/1, rng=42)

    [ ]:

    pp

    [ ]:

    qs = QStatistic(pp, nx=5, ny=5)

    [ ]:

    qs.plot()

    [ ]:

    qs.chi2_pvalue

    [ ]:

    qsw = QStatistic(pp, nx=5, ny=5, window=window, realizations=99, rng=42)

    [ ]:

    qsw.chi2_r_pvalue

    [ ]:

    qsw.plot()

## 5. Recap

In this notebook, we:

- represented a point pattern as an `(n, 2)` NumPy coordinate array,

- computed quadrat counts using rectangular and hexagonal tessellations over the MBB,

- used a Pearson chi-square statistic with equal expected counts across included cells,

- obtained either an analytical p-value or a Monte Carlo p-value from CSR realizations,

- clarified the issues associated with irregular windows, and

- ensured reproducibility of simulations via an explicit `rng`.