From 23814c29f382b0b2acd9b4b966c715506d4acfbf Mon Sep 17 00:00:00 2001
From: Dmitry Shkolnik <shkolnikd@gmail.com>
Date: Mon, 27 Jul 2026 21:16:45 -0700
Subject: [PATCH 1/2] Add homogeneity_test() for regime and spatial-lag
 transition inference

Tests whether transition dynamics differ across regimes. For
grd_spatial_markov objects this is the Rey, Kang & Wolf (2016) question:
do the lag-conditioned matrices differ from the pooled matrix? Reports:

- Pearson chi-squared Q and likelihood-ratio tests with the
  Bickenbach & Bald (2003) degrees-of-freedom adjustment
- the Kullback, Kupperman & Ku (1962) information test

Returns a tidy one-row-per-test tibble. A list method accepts arbitrary
regime count matrices (e.g. LISA Markov quadrant regimes or time
regimes), so the machinery is not tied to spatial conditioning.

Validated against PySAL giddy 2.3.9: exact agreement with
Spatial_Markov Q/LR/dof/p-values and giddy.markov.kullback() on the
usjoin spatial count matrices (fixture generated by
tools/oracle/pysal_giddy_inference.py), plus the published Kullback
(1962) Table 9.2 example as a shipped reference test.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 NAMESPACE                                     |   3 +
 NEWS.md                                       |  10 ++
 R/inference.R                                 | 144 ++++++++++++++++++
 _pkgdown.yml                                  |   3 +
 man/homogeneity_test.Rd                       |  55 +++++++
 .../pysal/pysal_spatial_homogeneity.csv       |   2 +
 tests/testthat/test-inference.R               |  96 ++++++++++++
 tests/testthat/test-oracles-optional.R        |  25 +++
 tools/oracle/pysal_giddy_inference.py         |  49 ++++++
 9 files changed, 387 insertions(+)
 create mode 100644 R/inference.R
 create mode 100644 man/homogeneity_test.Rd
 create mode 100644 tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv
 create mode 100644 tests/testthat/test-inference.R
 create mode 100644 tools/oracle/pysal_giddy_inference.py

diff --git a/NAMESPACE b/NAMESPACE
index fe42776..7cf69ed 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,6 +7,8 @@ S3method(as.data.frame,grd_spatial_markov)
 S3method(class_intervals,grd_classes)
 S3method(class_intervals,grd_markov)
 S3method(class_intervals,grd_spatial_markov)
+S3method(homogeneity_test,grd_spatial_markov)
+S3method(homogeneity_test,list)
 S3method(lag_intervals,grd_spatial_markov)
 S3method(print,grd_classes)
 S3method(print,grd_markov)
@@ -18,6 +20,7 @@ S3method(transition_matrix,grd_markov)
 S3method(transition_matrix,grd_spatial_markov)
 export(class_intervals)
 export(classify_dynamics)
+export(homogeneity_test)
 export(lag_intervals)
 export(markov_dynamics)
 export(plot_rank_mobility)
diff --git a/NEWS.md b/NEWS.md
index fda9478..8c64754 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,15 @@
 # griddy (development version)
 
+## New features
+
+- `homogeneity_test()` tests whether transition dynamics differ across
+  regimes -- for `grd_spatial_markov` objects, whether the matrices
+  conditioned on spatial-lag class differ from the pooled matrix (Rey, Kang,
+  and Wolf 2016). Reports the Pearson chi-squared Q and likelihood-ratio
+  tests with the Bickenbach--Bald (2003) degrees-of-freedom adjustment and
+  the Kullback (1962) information test, all validated against PySAL `giddy`
+  via static fixtures. A list method accepts arbitrary regime count matrices.
+
 # griddy 0.1.1
 
 ## Performance
diff --git a/R/inference.R b/R/inference.R
new file mode 100644
index 0000000..d65c9da
--- /dev/null
+++ b/R/inference.R
@@ -0,0 +1,144 @@
+#' Test Homogeneity Of Transition Matrices Across Regimes
+#'
+#' Tests whether transition dynamics are homogeneous across regimes -- for a
+#' spatial Markov object, whether the transition matrices conditioned on
+#' spatial-lag class genuinely differ from the pooled matrix (Rey, Kang, and
+#' Wolf 2016). Three statistics are reported: a Pearson chi-squared Q test and
+#' a likelihood-ratio test with the degrees-of-freedom adjustment of
+#' Bickenbach and Bald (2003), and the Kullback information test (Kullback,
+#' Kupperman, and Ku 1962). All mirror their PySAL `giddy` counterparts
+#' (`Spatial_Markov` homogeneity attributes and `giddy.markov.kullback()`).
+#'
+#' @param x A `grd_spatial_markov` object, or a list of regime transition
+#'   count matrices of identical dimensions (counts, not probabilities).
+#' @param ... Reserved for future methods.
+#'
+#' @return A tibble with one row per test (`chi2`, `likelihood_ratio`,
+#'   `kullback`) and columns `test`, `statistic`, `df`, and `p_value`. The
+#'   null hypothesis is that all regimes share the pooled transition matrix;
+#'   small p-values indicate the conditional dynamics differ.
+#'
+#' @examplesIf identical(Sys.getenv("IN_PKGDOWN"), "true")
+#' panel <- data.frame(
+#'   id = rep(1:9, each = 3),
+#'   year = rep(2020:2022, times = 9),
+#'   value = c(
+#'     1, 2, 3, 4, 3, 5, 6, 7, 6, 8, 9, 10, 11, 10,
+#'     12, 13, 14, 13, 15, 16, 17, 18, 17, 19, 20, 21, 20
+#'   )
+#' )
+#' grid <- sf::st_sf(
+#'   id = 1:9,
+#'   geometry = sf::st_make_grid(
+#'     sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 3, ymax = 3)),
+#'     n = c(3, 3)
+#'   )
+#' ) |>
+#'   dplyr::mutate(
+#'     nb = sfdep::st_contiguity(geometry),
+#'     wt = sfdep::st_weights(nb)
+#'   )
+#' spatial <- spatial_markov(panel, id, year, value, geometry = grid, k = 3)
+#' homogeneity_test(spatial)
+#' @export
+homogeneity_test <- function(x, ...) {
+  UseMethod("homogeneity_test")
+}
+
+#' @export
+homogeneity_test.grd_spatial_markov <- function(x, ...) {
+  mats <- lapply(x$lag_states, function(lag_class) {
+    transition_matrix(x, type = "count", lag_class = lag_class)
+  })
+  names(mats) <- x$lag_states
+  homogeneity_test(mats)
+}
+
+#' @export
+homogeneity_test.list <- function(x, ...) {
+  if (length(x) < 2) {
+    stop("Homogeneity tests need at least two regimes.", call. = FALSE)
+  }
+  mats <- lapply(x, function(m) {
+    m <- as.matrix(m)
+    storage.mode(m) <- "double"
+    m
+  })
+  dims <- unique(lapply(mats, dim))
+  if (length(dims) != 1 || dims[[1]][1] != dims[[1]][2]) {
+    stop(
+      "All regime matrices must be square and share the same dimensions.",
+      call. = FALSE
+    )
+  }
+  if (any(vapply(mats, function(m) anyNA(m) || any(m < 0), logical(1)))) {
+    stop(
+      "Regime matrices must contain non-negative transition counts.",
+      call. = FALSE
+    )
+  }
+  bb <- .grd_homogeneity(mats)
+  kb <- .grd_kullback(mats)
+  tibble::tibble(
+    test = c("chi2", "likelihood_ratio", "kullback"),
+    statistic = c(bb$Q, bb$LR, kb$statistic),
+    df = c(bb$dof, bb$dof, kb$dof),
+    p_value = c(bb$Q_p_value, bb$LR_p_value, kb$p_value)
+  )
+}
+
+# Pearson Q and likelihood-ratio homogeneity statistics with the
+# Bickenbach & Bald (2003) degrees-of-freedom adjustment, mirroring
+# giddy.markov.Homogeneity_Results.
+.grd_homogeneity <- function(mats) {
+  r <- nrow(mats[[1]])
+  total <- Reduce(`+`, mats)
+  n_i <- rowSums(total)
+  A_i <- rowSums(total > 0)
+  p_ij <- total / ifelse(n_i == 0, 1, n_i)
+  den <- p_ij + (p_ij == 0)
+
+  b_i <- numeric(r)
+  Q <- 0
+  LR <- 0
+  for (nijm in mats) {
+    nim <- rowSums(nijm)
+    b_i <- b_i + (nim > 0)
+    p_ijm <- nijm / ifelse(nim == 0, 1, nim)
+    Q <- Q + sum(nim * (p_ijm - p_ij)^2 / den)
+    mask <- (nijm > 0) & (p_ij > 0)
+    ratio <- ifelse(mask, p_ijm / p_ij, 1)
+    LR <- LR + sum(nijm * log(ratio))
+  }
+  LR <- 2 * LR
+  dof <- as.integer(sum((b_i - 1) * (A_i - 1)))
+  list(
+    Q = Q,
+    LR = LR,
+    dof = dof,
+    Q_p_value = stats::pchisq(Q, dof, lower.tail = FALSE),
+    LR_p_value = stats::pchisq(LR, dof, lower.tail = FALSE)
+  )
+}
+
+# Kullback information test of conditional homogeneity, mirroring
+# giddy.markov.kullback().
+.grd_kullback <- function(mats) {
+  s <- length(mats)
+  r <- nrow(mats[[1]])
+  xlogx <- function(v) sum(v * log(ifelse(v == 0, 1, v)))
+
+  t1 <- 2 * sum(vapply(mats, xlogx, numeric(1)))
+  pooled <- Reduce(`+`, mats)
+  t2 <- 2 * xlogx(pooled)
+  t3 <- 2 * xlogx(rowSums(pooled))
+  t4 <- 2 * sum(vapply(mats, function(m) xlogx(colSums(m)), numeric(1)))
+
+  statistic <- t1 - t4 - t2 + t3
+  dof <- as.integer(r * (s - 1) * (r - 1))
+  list(
+    statistic = statistic,
+    dof = dof,
+    p_value = stats::pchisq(statistic, dof, lower.tail = FALSE)
+  )
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index a99c477..2436860 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -49,6 +49,9 @@ reference:
       - spatial_markov
       - lag_intervals
       - rank_mobility
+  - title: Inference
+    contents:
+      - homogeneity_test
   - title: Plotting
     contents:
       - plot_transition_matrix
diff --git a/man/homogeneity_test.Rd b/man/homogeneity_test.Rd
new file mode 100644
index 0000000..991934d
--- /dev/null
+++ b/man/homogeneity_test.Rd
@@ -0,0 +1,55 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/inference.R
+\name{homogeneity_test}
+\alias{homogeneity_test}
+\title{Test Homogeneity Of Transition Matrices Across Regimes}
+\usage{
+homogeneity_test(x, ...)
+}
+\arguments{
+\item{x}{A \code{grd_spatial_markov} object, or a list of regime transition
+count matrices of identical dimensions (counts, not probabilities).}
+
+\item{...}{Reserved for future methods.}
+}
+\value{
+A tibble with one row per test (\code{chi2}, \code{likelihood_ratio},
+\code{kullback}) and columns \code{test}, \code{statistic}, \code{df}, and \code{p_value}. The
+null hypothesis is that all regimes share the pooled transition matrix;
+small p-values indicate the conditional dynamics differ.
+}
+\description{
+Tests whether transition dynamics are homogeneous across regimes -- for a
+spatial Markov object, whether the transition matrices conditioned on
+spatial-lag class genuinely differ from the pooled matrix (Rey, Kang, and
+Wolf 2016). Three statistics are reported: a Pearson chi-squared Q test and
+a likelihood-ratio test with the degrees-of-freedom adjustment of
+Bickenbach and Bald (2003), and the Kullback information test (Kullback,
+Kupperman, and Ku 1962). All mirror their PySAL \code{giddy} counterparts
+(\code{Spatial_Markov} homogeneity attributes and \code{giddy.markov.kullback()}).
+}
+\examples{
+\dontshow{if (identical(Sys.getenv("IN_PKGDOWN"), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+panel <- data.frame(
+  id = rep(1:9, each = 3),
+  year = rep(2020:2022, times = 9),
+  value = c(
+    1, 2, 3, 4, 3, 5, 6, 7, 6, 8, 9, 10, 11, 10,
+    12, 13, 14, 13, 15, 16, 17, 18, 17, 19, 20, 21, 20
+  )
+)
+grid <- sf::st_sf(
+  id = 1:9,
+  geometry = sf::st_make_grid(
+    sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 3, ymax = 3)),
+    n = c(3, 3)
+  )
+) |>
+  dplyr::mutate(
+    nb = sfdep::st_contiguity(geometry),
+    wt = sfdep::st_weights(nb)
+  )
+spatial <- spatial_markov(panel, id, year, value, geometry = grid, k = 3)
+homogeneity_test(spatial)
+\dontshow{\}) # examplesIf}
+}
diff --git a/tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv b/tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv
new file mode 100644
index 0000000..7067c81
--- /dev/null
+++ b/tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv
@@ -0,0 +1,2 @@
+Q,Q_p_value,LR,LR_p_value,dof,kullback,kullback_dof,kullback_p_value
+200.62372763076806,0.0,170.6588284308426,1.492250767398673e-12,60,230.0266246375395,80,2.220446049250313e-16
diff --git a/tests/testthat/test-inference.R b/tests/testthat/test-inference.R
new file mode 100644
index 0000000..8e45143
--- /dev/null
+++ b/tests/testthat/test-inference.R
@@ -0,0 +1,96 @@
+nine_cell_spatial <- function() {
+  panel <- data.frame(
+    id = rep(1:9, each = 3),
+    year = rep(2020:2022, times = 9),
+    value = c(
+      1, 2, 3, 4, 3, 5, 6, 7, 6, 8, 9, 10, 11, 10,
+      12, 13, 14, 13, 15, 16, 17, 18, 17, 19, 20, 21, 20
+    )
+  )
+  grid <- sf::st_sf(
+    id = 1:9,
+    geometry = sf::st_make_grid(
+      sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 3, ymax = 3)),
+      n = c(3, 3)
+    )
+  ) |>
+    dplyr::mutate(
+      nb = sfdep::st_contiguity(geometry),
+      wt = sfdep::st_weights(nb)
+    )
+  spatial_markov(panel, id, year, value, geometry = grid, k = 3)
+}
+
+test_that("kullback statistic reproduces the published Kullback (1962) example", {
+  s1 <- matrix(
+    c(
+      22, 11, 24, 2, 2, 7,
+      5, 23, 15, 3, 42, 6,
+      4, 21, 190, 25, 20, 34,
+      0, 2, 14, 56, 14, 28,
+      32, 15, 20, 10, 56, 14,
+      5, 22, 31, 18, 13, 134
+    ),
+    nrow = 6, byrow = TRUE
+  )
+  s2 <- matrix(
+    c(
+      3, 6, 9, 3, 0, 8,
+      1, 9, 3, 12, 27, 5,
+      2, 9, 208, 32, 5, 18,
+      0, 14, 32, 108, 40, 40,
+      22, 14, 9, 26, 224, 14,
+      1, 5, 13, 53, 13, 116
+    ),
+    nrow = 6, byrow = TRUE
+  )
+
+  out <- homogeneity_test(list(s1, s2))
+  kb <- out[out$test == "kullback", ]
+
+  expect_equal(kb$statistic, 160.961, tolerance = 1e-5)
+  expect_equal(kb$df, 30L)
+  expect_lt(kb$p_value, 1e-6)
+})
+
+test_that("identical regimes are homogeneous under the chi2 and LR tests", {
+  m <- matrix(c(10, 5, 3, 12), nrow = 2)
+  out <- homogeneity_test(list(m, m))
+  bb <- out[out$test != "kullback", ]
+
+  expect_equal(bb$statistic, rep(0, 2), tolerance = 1e-10)
+  expect_equal(bb$p_value, rep(1, 2))
+  # The Kullback statistic is not exactly zero for identical regimes (its
+  # marginal term uses terminal-state totals); giddy agrees on this value.
+  expect_equal(
+    out$statistic[out$test == "kullback"],
+    1.0698498515574215,
+    tolerance = 1e-10
+  )
+})
+
+test_that("homogeneity_test returns a tidy three-row tibble", {
+  spatial <- nine_cell_spatial()
+  out <- homogeneity_test(spatial)
+
+  expect_s3_class(out, "tbl_df")
+  expect_equal(out$test, c("chi2", "likelihood_ratio", "kullback"))
+  expect_equal(names(out), c("test", "statistic", "df", "p_value"))
+  expect_true(all(out$statistic >= 0))
+  expect_true(all(out$df >= 0))
+  expect_true(all(out$p_value >= 0 & out$p_value <= 1))
+})
+
+test_that("invalid inputs are rejected", {
+  m <- matrix(c(10, 5, 3, 12), nrow = 2)
+
+  expect_error(homogeneity_test(list(m)), "at least two regimes")
+  expect_error(
+    homogeneity_test(list(m, matrix(0, 3, 3))),
+    "same dimensions"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(-1, 2, 3, 4), 2))),
+    "non-negative"
+  )
+})
diff --git a/tests/testthat/test-oracles-optional.R b/tests/testthat/test-oracles-optional.R
index ddc7aaf..0e0c591 100644
--- a/tests/testthat/test-oracles-optional.R
+++ b/tests/testthat/test-oracles-optional.R
@@ -95,6 +95,31 @@ test_that("estdaR spatial Markov agrees on bundled state example when installed"
   expect_equal(unname(grd_counts), unname(prior$Transitions))
 })
 
+test_that("homogeneity_test agrees with the PySAL giddy oracle fixtures", {
+  mats <- lapply(1:5, function(idx) {
+    as.matrix(utils::read.csv(
+      test_path("fixtures", "pysal", sprintf("pysal_spatial_counts_lag_%d.csv", idx))
+    ))
+  })
+  fixture <- utils::read.csv(
+    test_path("fixtures", "pysal", "pysal_spatial_homogeneity.csv")
+  )
+
+  out <- homogeneity_test(mats)
+
+  expect_equal(out$statistic[out$test == "chi2"], fixture$Q)
+  expect_equal(out$p_value[out$test == "chi2"], fixture$Q_p_value)
+  expect_equal(out$statistic[out$test == "likelihood_ratio"], fixture$LR)
+  expect_equal(out$p_value[out$test == "likelihood_ratio"], fixture$LR_p_value)
+  expect_equal(out$df[out$test != "kullback"], rep(fixture$dof, 2))
+  expect_equal(out$statistic[out$test == "kullback"], fixture$kullback)
+  expect_equal(out$df[out$test == "kullback"], fixture$kullback_dof)
+  expect_lt(
+    abs(out$p_value[out$test == "kullback"] - fixture$kullback_p_value),
+    1e-15
+  )
+})
+
 test_that("PySAL giddy oracle comparison is optional", {
   skip_if(Sys.which("python3") == "", "python3 not available")
   code <- "import importlib.util, sys; sys.exit(0 if importlib.util.find_spec('giddy') else 1)"
diff --git a/tools/oracle/pysal_giddy_inference.py b/tools/oracle/pysal_giddy_inference.py
new file mode 100644
index 0000000..4ff71dc
--- /dev/null
+++ b/tools/oracle/pysal_giddy_inference.py
@@ -0,0 +1,49 @@
+"""Optional PySAL giddy inference oracle generation.
+
+Development infrastructure only; excluded from the package build. Run in an
+environment with giddy and libpysal installed to generate static homogeneity
+test fixtures for R tests. Uses the same usjoin setup as
+pysal_giddy_oracle.py so the fixtures share the spatial count matrices.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+
+def main() -> None:
+    try:
+        import libpysal
+        from giddy.markov import Spatial_Markov, kullback
+    except Exception as exc:  # pragma: no cover - development helper
+        raise SystemExit(f"Install libpysal and giddy first: {exc}") from exc
+
+    f = libpysal.io.open(libpysal.examples.get_path("usjoin.csv"))
+    pci = np.array([f.by_col[str(y)] for y in range(1929, 2010)]).T
+    rpci = pci / pci.mean(axis=0)
+    w = libpysal.io.open(libpysal.examples.get_path("states48.gal")).read()
+    w.transform = "r"
+
+    spatial = Spatial_Markov(rpci, w, fixed=True, k=5, m=5)
+    kb = kullback(spatial.T)
+
+    out_dir = Path("tests/testthat/fixtures/pysal")
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    pd.DataFrame(
+        {
+            "Q": [spatial.Q],
+            "Q_p_value": [spatial.Q_p_value],
+            "LR": [spatial.LR],
+            "LR_p_value": [spatial.LR_p_value],
+            "dof": [spatial.dof_hom],
+            "kullback": [kb["Conditional homogeneity"]],
+            "kullback_dof": [kb["Conditional homogeneity dof"]],
+            "kullback_p_value": [kb["Conditional homogeneity pvalue"]],
+        }
+    ).to_csv(out_dir / "pysal_spatial_homogeneity.csv", index=False)
+
+
+if __name__ == "__main__":
+    main()

From f0100e7a929f402b3820de5ed2774ee2dcee7645 Mon Sep 17 00:00:00 2001
From: Dmitry Shkolnik <shkolnikd@gmail.com>
Date: Sat, 1 Aug 2026 17:26:52 -0700
Subject: [PATCH 2/2] Correct Kullback homogeneity margins

---
 NEWS.md                                       |  8 +++-
 R/inference.R                                 | 43 ++++++++++++++++---
 man/homogeneity_test.Rd                       |  7 ++-
 .../pysal/pysal_spatial_homogeneity.csv       |  2 +-
 tests/testthat/test-inference.R               | 39 +++++++++++------
 tests/testthat/test-oracles-optional.R        |  4 +-
 tools/oracle/pysal_giddy_inference.py         | 31 ++++++++++---
 7 files changed, 104 insertions(+), 30 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 8c64754..3f0a915 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,8 +7,12 @@
   conditioned on spatial-lag class differ from the pooled matrix (Rey, Kang,
   and Wolf 2016). Reports the Pearson chi-squared Q and likelihood-ratio
   tests with the Bickenbach--Bald (2003) degrees-of-freedom adjustment and
-  the Kullback (1962) information test, all validated against PySAL `giddy`
-  via static fixtures. A list method accepts arbitrary regime count matrices.
+  the Kullback (1962) information test. The first two are validated against
+  PySAL `giddy` via static fixtures. The Kullback statistic uses the
+  origin-state margins required by the conditional likelihood, deliberately
+  differing from `giddy.markov.kullback()` through giddy 2.3.9 when transition
+  tables are not flow-balanced. A list method accepts arbitrary regime count
+  matrices.
 
 # griddy 0.1.1
 
diff --git a/R/inference.R b/R/inference.R
index d65c9da..d6b9a6d 100644
--- a/R/inference.R
+++ b/R/inference.R
@@ -6,8 +6,11 @@
 #' Wolf 2016). Three statistics are reported: a Pearson chi-squared Q test and
 #' a likelihood-ratio test with the degrees-of-freedom adjustment of
 #' Bickenbach and Bald (2003), and the Kullback information test (Kullback,
-#' Kupperman, and Ku 1962). All mirror their PySAL `giddy` counterparts
-#' (`Spatial_Markov` homogeneity attributes and `giddy.markov.kullback()`).
+#' Kupperman, and Ku 1962). The Pearson and likelihood-ratio statistics mirror
+#' the PySAL `giddy` `Spatial_Markov` homogeneity attributes. The Kullback
+#' statistic uses origin-state (row) margins; this deliberately differs from
+#' `giddy.markov.kullback()` through giddy 2.3.9, which uses terminal-state
+#' (column) margins for the stratum-specific term.
 #'
 #' @param x A `grd_spatial_markov` object, or a list of regime transition
 #'   count matrices of identical dimensions (counts, not probabilities).
@@ -61,6 +64,12 @@ homogeneity_test.list <- function(x, ...) {
   }
   mats <- lapply(x, function(m) {
     m <- as.matrix(m)
+    if (!is.numeric(m)) {
+      stop(
+        "Regime matrices must contain numeric transition counts.",
+        call. = FALSE
+      )
+    }
     storage.mode(m) <- "double"
     m
   })
@@ -71,12 +80,30 @@ homogeneity_test.list <- function(x, ...) {
       call. = FALSE
     )
   }
-  if (any(vapply(mats, function(m) anyNA(m) || any(m < 0), logical(1)))) {
+  if (any(vapply(mats, function(m) any(!is.finite(m)), logical(1)))) {
+    stop(
+      "Regime matrices must contain finite transition counts.",
+      call. = FALSE
+    )
+  }
+  if (any(vapply(mats, function(m) any(m < 0), logical(1)))) {
     stop(
       "Regime matrices must contain non-negative transition counts.",
       call. = FALSE
     )
   }
+  if (any(vapply(mats, function(m) any(m != floor(m)), logical(1)))) {
+    stop(
+      "Regime matrices must contain whole-number transition counts.",
+      call. = FALSE
+    )
+  }
+  if (any(vapply(mats, function(m) sum(m) == 0, logical(1)))) {
+    stop(
+      "Each regime matrix must contain at least one transition.",
+      call. = FALSE
+    )
+  }
   bb <- .grd_homogeneity(mats)
   kb <- .grd_kullback(mats)
   tibble::tibble(
@@ -121,8 +148,10 @@ homogeneity_test.list <- function(x, ...) {
   )
 }
 
-# Kullback information test of conditional homogeneity, mirroring
-# giddy.markov.kullback().
+# Kullback information test of conditional homogeneity. The stratum-specific
+# nuisance margins are origin-state (row) totals. giddy <= 2.3.9 instead sums
+# terminal-state columns here; the distinction vanishes for flow-balanced
+# transition tables but not for aggregated panels with endpoint imbalance.
 .grd_kullback <- function(mats) {
   s <- length(mats)
   r <- nrow(mats[[1]])
@@ -132,9 +161,9 @@ homogeneity_test.list <- function(x, ...) {
   pooled <- Reduce(`+`, mats)
   t2 <- 2 * xlogx(pooled)
   t3 <- 2 * xlogx(rowSums(pooled))
-  t4 <- 2 * sum(vapply(mats, function(m) xlogx(colSums(m)), numeric(1)))
+  t4 <- 2 * sum(vapply(mats, function(m) xlogx(rowSums(m)), numeric(1)))
 
-  statistic <- t1 - t4 - t2 + t3
+  statistic <- max(0, t1 - t4 - t2 + t3)
   dof <- as.integer(r * (s - 1) * (r - 1))
   list(
     statistic = statistic,
diff --git a/man/homogeneity_test.Rd b/man/homogeneity_test.Rd
index 991934d..930d488 100644
--- a/man/homogeneity_test.Rd
+++ b/man/homogeneity_test.Rd
@@ -25,8 +25,11 @@ spatial-lag class genuinely differ from the pooled matrix (Rey, Kang, and
 Wolf 2016). Three statistics are reported: a Pearson chi-squared Q test and
 a likelihood-ratio test with the degrees-of-freedom adjustment of
 Bickenbach and Bald (2003), and the Kullback information test (Kullback,
-Kupperman, and Ku 1962). All mirror their PySAL \code{giddy} counterparts
-(\code{Spatial_Markov} homogeneity attributes and \code{giddy.markov.kullback()}).
+Kupperman, and Ku 1962). The Pearson and likelihood-ratio statistics mirror
+the PySAL \code{giddy} \code{Spatial_Markov} homogeneity attributes. The Kullback
+statistic uses origin-state (row) margins; this deliberately differs from
+\code{giddy.markov.kullback()} through giddy 2.3.9, which uses terminal-state
+(column) margins for the stratum-specific term.
 }
 \examples{
 \dontshow{if (identical(Sys.getenv("IN_PKGDOWN"), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
diff --git a/tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv b/tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv
index 7067c81..a095eca 100644
--- a/tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv
+++ b/tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv
@@ -1,2 +1,2 @@
 Q,Q_p_value,LR,LR_p_value,dof,kullback,kullback_dof,kullback_p_value
-200.62372763076806,0.0,170.6588284308426,1.492250767398673e-12,60,230.0266246375395,80,2.220446049250313e-16
+200.62372763076806,0.0,170.6588284308426,1.492250767398673e-12,60,170.65882843085274,80,1.5963658564297747e-08
diff --git a/tests/testthat/test-inference.R b/tests/testthat/test-inference.R
index 8e45143..def559d 100644
--- a/tests/testthat/test-inference.R
+++ b/tests/testthat/test-inference.R
@@ -53,20 +53,15 @@ test_that("kullback statistic reproduces the published Kullback (1962) example",
   expect_lt(kb$p_value, 1e-6)
 })
 
-test_that("identical regimes are homogeneous under the chi2 and LR tests", {
+test_that("identical regimes are homogeneous under every test", {
   m <- matrix(c(10, 5, 3, 12), nrow = 2)
-  out <- homogeneity_test(list(m, m))
-  bb <- out[out$test != "kullback", ]
 
-  expect_equal(bb$statistic, rep(0, 2), tolerance = 1e-10)
-  expect_equal(bb$p_value, rep(1, 2))
-  # The Kullback statistic is not exactly zero for identical regimes (its
-  # marginal term uses terminal-state totals); giddy agrees on this value.
-  expect_equal(
-    out$statistic[out$test == "kullback"],
-    1.0698498515574215,
-    tolerance = 1e-10
-  )
+  for (scale in c(1, 10, 100, 1000)) {
+    out <- homogeneity_test(list(m * scale, m * scale))
+
+    expect_equal(out$statistic, rep(0, 3), tolerance = 1e-9)
+    expect_equal(out$p_value, rep(1, 3))
+  }
 })
 
 test_that("homogeneity_test returns a tidy three-row tibble", {
@@ -93,4 +88,24 @@ test_that("invalid inputs are rejected", {
     homogeneity_test(list(m, matrix(c(-1, 2, 3, 4), 2))),
     "non-negative"
   )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(1, 2, 3, Inf), 2))),
+    "finite"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(1, 2, 3, NA), 2))),
+    "finite"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(1, 2, 3, 4.5), 2))),
+    "whole-number"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(0, 2, 2))),
+    "at least one transition"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(letters[1:4], 2))),
+    "numeric"
+  )
 })
diff --git a/tests/testthat/test-oracles-optional.R b/tests/testthat/test-oracles-optional.R
index 0e0c591..39558c4 100644
--- a/tests/testthat/test-oracles-optional.R
+++ b/tests/testthat/test-oracles-optional.R
@@ -95,7 +95,7 @@ test_that("estdaR spatial Markov agrees on bundled state example when installed"
   expect_equal(unname(grd_counts), unname(prior$Transitions))
 })
 
-test_that("homogeneity_test agrees with the PySAL giddy oracle fixtures", {
+test_that("homogeneity_test agrees with independently generated fixtures", {
   mats <- lapply(1:5, function(idx) {
     as.matrix(utils::read.csv(
       test_path("fixtures", "pysal", sprintf("pysal_spatial_counts_lag_%d.csv", idx))
@@ -112,6 +112,8 @@ test_that("homogeneity_test agrees with the PySAL giddy oracle fixtures", {
   expect_equal(out$statistic[out$test == "likelihood_ratio"], fixture$LR)
   expect_equal(out$p_value[out$test == "likelihood_ratio"], fixture$LR_p_value)
   expect_equal(out$df[out$test != "kullback"], rep(fixture$dof, 2))
+  # Kullback uses an independent origin-margin calculation rather than
+  # giddy.markov.kullback(), whose stratum term uses column margins in 2.3.9.
   expect_equal(out$statistic[out$test == "kullback"], fixture$kullback)
   expect_equal(out$df[out$test == "kullback"], fixture$kullback_dof)
   expect_lt(
diff --git a/tools/oracle/pysal_giddy_inference.py b/tools/oracle/pysal_giddy_inference.py
index 4ff71dc..625b8a1 100644
--- a/tools/oracle/pysal_giddy_inference.py
+++ b/tools/oracle/pysal_giddy_inference.py
@@ -10,12 +10,33 @@
 
 import numpy as np
 import pandas as pd
+from scipy.stats import chi2
+
+
+def kullback_with_origin_margins(F: np.ndarray) -> tuple[float, int, float]:
+    """Kullback conditional-homogeneity test using initial-state totals."""
+
+    def xlogx(values: np.ndarray) -> float:
+        positive = values > 0
+        return float(np.sum(values[positive] * np.log(values[positive])))
+
+    pooled = F.sum(axis=0)
+    statistic = 2 * (
+        sum(xlogx(regime) for regime in F)
+        - sum(xlogx(regime.sum(axis=1)) for regime in F)
+        - xlogx(pooled)
+        + xlogx(pooled.sum(axis=1))
+    )
+    statistic = max(0.0, statistic)
+    strata, states, _ = F.shape
+    dof = states * (strata - 1) * (states - 1)
+    return statistic, dof, float(chi2.sf(statistic, dof))
 
 
 def main() -> None:
     try:
         import libpysal
-        from giddy.markov import Spatial_Markov, kullback
+        from giddy.markov import Spatial_Markov
     except Exception as exc:  # pragma: no cover - development helper
         raise SystemExit(f"Install libpysal and giddy first: {exc}") from exc
 
@@ -26,7 +47,7 @@ def main() -> None:
     w.transform = "r"
 
     spatial = Spatial_Markov(rpci, w, fixed=True, k=5, m=5)
-    kb = kullback(spatial.T)
+    kb_statistic, kb_dof, kb_p_value = kullback_with_origin_margins(spatial.T)
 
     out_dir = Path("tests/testthat/fixtures/pysal")
     out_dir.mkdir(parents=True, exist_ok=True)
@@ -38,9 +59,9 @@ def main() -> None:
             "LR": [spatial.LR],
             "LR_p_value": [spatial.LR_p_value],
             "dof": [spatial.dof_hom],
-            "kullback": [kb["Conditional homogeneity"]],
-            "kullback_dof": [kb["Conditional homogeneity dof"]],
-            "kullback_p_value": [kb["Conditional homogeneity pvalue"]],
+            "kullback": [kb_statistic],
+            "kullback_dof": [kb_dof],
+            "kullback_p_value": [kb_p_value],
         }
     ).to_csv(out_dir / "pysal_spatial_homogeneity.csv", index=False)