From 1dfa05918ac63102d446b91d1f7d2e9aee512922 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Fri, 14 Jul 2023 14:06:12 -0700
Subject: [PATCH 01/11] add ICARRV and ICAR

---
 pymc/distributions/multivariate.py | 211 +++++++++++++++++++++++++++++
 1 file changed, 211 insertions(+)

diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index ba613967c1..dcc163be16 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -2218,6 +2218,217 @@ def logp(value, mu, W, alpha, tau):
         )
 
 
+class ICARRV(RandomVariable):
+    name = "icar"
+    ndim_supp = 1
+    ndims_params = [2, 1, 1, 0, 0, 0]
+    dtype = "floatX"
+    _print_name = ("ICAR", "\\operatorname{ICAR}")
+
+    def __call__(
+        self, W, node1, node2, sigma=1, N=None, centering_strength=0.001, size=None, **kwargs
+    ):
+        return super().__call__(W, node1, node2, sigma, N, centering_strength, size=size, **kwargs)
+
+    @classmethod
+    def rng_fn(cls, rng, size, W, node1, node2, sigma, N, centering_strength):
+        raise NotImplementedError("Cannot sample from ICAR prior")
+
+
+icar = ICARRV()
+
+
+class ICAR(Continuous):
+    r"""
+    The intrinsic conditional autoregressive prior. It is primarily used to model
+    covariance between neighboring areas on large datasets. It is a special case
+    of the :class:`~pymc.CAR` distribution where alpha is set to 1.
+
+    The log probability density function is
+
+    .. math::
+        f(\\phi| W,\\sigma) =
+          -\frac{1}{2\\sigma^{2}} \\sum_{i\\sim j} (\\phi_{i} - \\phi_{j})^2 -
+          \frac{1}{2}*\frac{\\sum_{i}{\\phi_{i}}}{0.001N}^{2} - \\ln{\\sqrt{2\\pi}} -
+          \\ln{0.001N}
+
+    The first term represents the spatial covariance component. Each $\\phi_{i}$ is penalized
+    based on the square distance from each of its neighbors. The notation $i\\sim j$
+    indicates a sum over all the neighbors of $\\phi_{i}$. The last three terms are the
+    Normal log density function where the mean is zero and the standard deviation is
+    $N * 0.001$ (where N is the length of the vector $\\phi$). This component imposed the zero-sum
+    constraint by finding the sum of the vector $\\phi$ and penalizing based on its
+    distance from zero.
+
+    ========  ==========================================
+    Support   :math:`x \\in \\mathbb{R}^k` ?
+    Mean      :math:`0` ?
+    Variance  :math:`T^{-1}` ?
+    ========  ==========================================
+
+    Parameters
+    ----------
+    W : ndarray of int, optional
+        Symmetric adjacency matrix of 1s and 0s indicating adjacency between elements.
+        Must pass either W or both node1 and node2.
+
+    node1 : tensor_like of int, optional
+        The edgelist. Must be passed with node2 as well. See example below to
+        illustrate the relationship between the adjacency matrix and the edgelist.
+
+    node2 : tensor_like of int, optional
+        Must be passed with node1 as well.
+
+    sigma : scalar, default 1
+        Standard deviation of the vector of phi's. Putting a prior on sigma
+        will result in a centered parameterization. In most cases, it is
+        preferable to use a non-centered parameterization by using the default
+        value and multiplying the resulting phi's by sigma. See the example below.
+
+    centering_strength : scalar, default 0.001
+        Controls how strongly to enforce the zero-sum constraint. It sets the
+        standard deviation of a normal density function with mean zero.
+
+
+    Examples
+    --------
+    This example illustrates how to switch between centered and non-centered
+    parameterizations.
+
+    .. code-block:: python
+
+        # 4x4 adjacency matrix
+        # arranged in a square lattice
+
+        W = np.array([[0,1,0,1],
+                      [1,0,1,0],
+                      [0,1,0,1],
+                      [1,0,1,0]])
+
+        # centered parameterization
+
+        with pm.Model():
+            sigma = pm.Exponential('sigma',1)
+            phi = pm.ICAR('phi',W=W,sigma=sigma)
+
+            mu = phi
+
+        # non-centered parameterization
+
+        with pm.Model():
+            sigma = pm.Exponential('sigma',1)
+            phi = pm.ICAR('phi',W=W)
+
+            mu = sigma * phi
+
+    This example illustrates how to switch between adjacency matrices and
+    edge lists.
+
+    .. code-block:: python
+
+        # phi is still arranged
+        # in a square lattice
+        # as above. But now we
+        # represent the lattice
+        # through an edgelist
+
+        node1 = np.array([0,0,1,2])
+        node2 = np.array([1,3,2,3])
+
+        with pm.Model():
+            phi = pm.ICAR('phi',node1=node1,node2=node2)
+
+
+    References
+    ----------
+    ..  Mitzi, M., Wheeler-Martin, K., Simpson, D., Mooney, J. S.,
+        Gelman, A., Dimaggio, C.
+        "Bayesian hierarchical spatial models: Implementing the Besag York
+        Mollié model in stan"
+        Spatial and Spatio-temporal Epidemiology, Vol. 31, (Aug., 2019),
+        pp 1-18
+
+
+    """
+
+    rv_op = icar
+
+    @classmethod
+    def dist(
+        cls, W=None, node1=None, node2=None, sigma=1, N=None, centering_strength=0.001, **kwargs
+    ):
+        # check that they use either an adjacency matrix or edgelist
+        # specification but not both
+
+        if W is not None and (node1 is not None or node2 is not None):
+            raise ValueError(
+                "Cannot pass both an adjacency matrix (W) and an edgelist (node1, node2)"
+            )
+
+        if W is None and (node1 is None and node2 is None):
+            raise ValueError(
+                "Must pass either an adjacency matrix (W) or an edgelist (node1, node2)"
+            )
+
+        if W is None and (node1 is None or node2 is None):
+            raise ValueError("node1 must be passed with node2")
+
+        # check that adjacency matrix is two dimensional
+
+        if W is not None and not W.ndim == 2:
+            raise ValueError("W must be a matrix (ndim=2).")
+
+        # convert adjacency matrix to edgelist representation
+
+        if W is not None:
+            node1, node2 = np.where(np.tril(W) == 1)
+
+        node1 = pt.as_tensor_variable(node1, dtype=int)
+        node2 = pt.as_tensor_variable(node2, dtype=int)
+
+        # check symmetry of adjacency matrix
+
+        if W is not None:
+            W = pt.as_tensor_variable(floatX(W))
+            msg = "W must be a symmetric adjacency matrix."
+            W = Assert(msg)(W, pt.allclose(W, W.T))
+
+        # need some way of figuring out the sample size
+        # if they do not pass W. Or we only allow
+        # users to pass the adjacency matrix.
+
+        N = pt.shape(W)[0]
+
+        # check on sigma
+
+        sigma = pt.as_tensor_variable(floatX(sigma))
+        sigma = Assert("sigma > 0")(sigma, pt.gt(sigma, 0))
+
+        # check on centering_strength
+
+        centering_strength = pt.as_tensor_variable(floatX(centering_strength))
+        centering_strength = Assert("centering_strength > 0")(
+            centering_strength, pt.gt(centering_strength, 0)
+        )
+
+        return super().dist([W, node1, node2, sigma, N, centering_strength], **kwargs)
+
+    def moment(rv, size, W, node1, node2, sigma, N, centering_strength):
+        return pt.zeros(pt.shape(W)[1])
+
+    def logp(value, W, node1, node2, sigma, N, centering_strength):
+        pairwise_difference = (-1 / (2 * sigma**2)) * pt.sum(
+            pt.square(value[node1] - value[node2])
+        )
+        soft_center = (
+            -0.5 * pt.pow(pt.sum(value) / (centering_strength * N), 2)
+            - pt.log(pt.sqrt(2.0 * np.pi))
+            - pt.log(centering_strength * N)
+        )
+
+        return check_parameters(pairwise_difference + soft_center, sigma > 0, msg="sigma > 0")
+
+
 class StickBreakingWeightsRV(RandomVariable):
     name = "stick_breaking_weights"
     ndim_supp = 1

From 10fed3b7c175717f38e4e03daa42cebc43e5b770 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Mon, 17 Jul 2023 15:51:52 -0700
Subject: [PATCH 02/11] reorganized API so W is required, added tests on
 logp,rng,checks

---
 pymc/distributions/multivariate.py       | 93 +++++++-----------------
 tests/distributions/test_multivariate.py | 45 ++++++++++++
 2 files changed, 70 insertions(+), 68 deletions(-)

diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index dcc163be16..b92cf97139 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -2225,13 +2225,11 @@ class ICARRV(RandomVariable):
     dtype = "floatX"
     _print_name = ("ICAR", "\\operatorname{ICAR}")
 
-    def __call__(
-        self, W, node1, node2, sigma=1, N=None, centering_strength=0.001, size=None, **kwargs
-    ):
-        return super().__call__(W, node1, node2, sigma, N, centering_strength, size=size, **kwargs)
+    def __call__(self, W, node1, node2, N, sigma, centering_strength, size=None, **kwargs):
+        return super().__call__(W, node1, node2, N, sigma, centering_strength, size=size, **kwargs)
 
     @classmethod
-    def rng_fn(cls, rng, size, W, node1, node2, sigma, N, centering_strength):
+    def rng_fn(cls, rng, size, W, node1, node2, N, sigma, centering_strength):
         raise NotImplementedError("Cannot sample from ICAR prior")
 
 
@@ -2239,7 +2237,7 @@ def rng_fn(cls, rng, size, W, node1, node2, sigma, N, centering_strength):
 
 
 class ICAR(Continuous):
-    r"""
+    """
     The intrinsic conditional autoregressive prior. It is primarily used to model
     covariance between neighboring areas on large datasets. It is a special case
     of the :class:`~pymc.CAR` distribution where alpha is set to 1.
@@ -2268,17 +2266,10 @@ class ICAR(Continuous):
 
     Parameters
     ----------
-    W : ndarray of int, optional
+    W : ndarray of int
         Symmetric adjacency matrix of 1s and 0s indicating adjacency between elements.
         Must pass either W or both node1 and node2.
 
-    node1 : tensor_like of int, optional
-        The edgelist. Must be passed with node2 as well. See example below to
-        illustrate the relationship between the adjacency matrix and the edgelist.
-
-    node2 : tensor_like of int, optional
-        Must be passed with node1 as well.
-
     sigma : scalar, default 1
         Standard deviation of the vector of phi's. Putting a prior on sigma
         will result in a centered parameterization. In most cases, it is
@@ -2297,6 +2288,9 @@ class ICAR(Continuous):
 
     .. code-block:: python
 
+        import numpy as np
+        import pymc as pm
+
         # 4x4 adjacency matrix
         # arranged in a square lattice
 
@@ -2321,24 +2315,6 @@ class ICAR(Continuous):
 
             mu = sigma * phi
 
-    This example illustrates how to switch between adjacency matrices and
-    edge lists.
-
-    .. code-block:: python
-
-        # phi is still arranged
-        # in a square lattice
-        # as above. But now we
-        # represent the lattice
-        # through an edgelist
-
-        node1 = np.array([0,0,1,2])
-        node2 = np.array([1,3,2,3])
-
-        with pm.Model():
-            phi = pm.ICAR('phi',node1=node1,node2=node2)
-
-
     References
     ----------
     ..  Mitzi, M., Wheeler-Martin, K., Simpson, D., Mooney, J. S.,
@@ -2354,50 +2330,31 @@ class ICAR(Continuous):
     rv_op = icar
 
     @classmethod
-    def dist(
-        cls, W=None, node1=None, node2=None, sigma=1, N=None, centering_strength=0.001, **kwargs
-    ):
-        # check that they use either an adjacency matrix or edgelist
-        # specification but not both
+    def dist(cls, W, sigma=1, centering_strength=0.001, **kwargs):
+        # check that adjacency matrix is two dimensional,
+        # square
+        # and symmetrical
 
-        if W is not None and (node1 is not None or node2 is not None):
-            raise ValueError(
-                "Cannot pass both an adjacency matrix (W) and an edgelist (node1, node2)"
-            )
-
-        if W is None and (node1 is None and node2 is None):
-            raise ValueError(
-                "Must pass either an adjacency matrix (W) or an edgelist (node1, node2)"
-            )
-
-        if W is None and (node1 is None or node2 is None):
-            raise ValueError("node1 must be passed with node2")
+        if not W.ndim == 2:
+            raise ValueError("W must be matrix with ndim=2")
 
-        # check that adjacency matrix is two dimensional
+        if not W.shape[0] == W.shape[1]:
+            raise ValueError("W must be a square matrix")
 
-        if W is not None and not W.ndim == 2:
-            raise ValueError("W must be a matrix (ndim=2).")
+        if not np.allclose(W.T, W):
+            raise ValueError("W must be a symmetric matrix")
 
         # convert adjacency matrix to edgelist representation
 
-        if W is not None:
-            node1, node2 = np.where(np.tril(W) == 1)
+        node1, node2 = np.where(np.tril(W) == 1)
 
         node1 = pt.as_tensor_variable(node1, dtype=int)
         node2 = pt.as_tensor_variable(node2, dtype=int)
 
-        # check symmetry of adjacency matrix
-
-        if W is not None:
-            W = pt.as_tensor_variable(floatX(W))
-            msg = "W must be a symmetric adjacency matrix."
-            W = Assert(msg)(W, pt.allclose(W, W.T))
-
-        # need some way of figuring out the sample size
-        # if they do not pass W. Or we only allow
-        # users to pass the adjacency matrix.
+        W = pt.as_tensor_variable(W, dtype=int)
 
         N = pt.shape(W)[0]
+        N = pt.as_tensor_variable(N)
 
         # check on sigma
 
@@ -2411,12 +2368,12 @@ def dist(
             centering_strength, pt.gt(centering_strength, 0)
         )
 
-        return super().dist([W, node1, node2, sigma, N, centering_strength], **kwargs)
+        return super().dist([W, node1, node2, N, sigma, centering_strength], **kwargs)
 
-    def moment(rv, size, W, node1, node2, sigma, N, centering_strength):
-        return pt.zeros(pt.shape(W)[1])
+    def moment(rv, size, W, node1, node2, N, sigma, centering_strength):
+        return pt.zeros(N)
 
-    def logp(value, W, node1, node2, sigma, N, centering_strength):
+    def logp(value, W, node1, node2, N, sigma, centering_strength):
         pairwise_difference = (-1 / (2 * sigma**2)) * pt.sum(
             pt.square(value[node1] - value[node2])
         )
diff --git a/tests/distributions/test_multivariate.py b/tests/distributions/test_multivariate.py
index 640d7484c9..0f4d236317 100644
--- a/tests/distributions/test_multivariate.py
+++ b/tests/distributions/test_multivariate.py
@@ -2070,6 +2070,51 @@ def check_draws_match_expected(self):
         assert np.all(np.abs(draw(x, random_seed=rng) - np.array([0.5, 0, 2.0])) < 0.01)
 
 
+# checks on adjacency matrix
+
+
+@pytest.mark.parametrize(
+    "W,msg",
+    [
+        (np.array([0, 1, 0, 0]), "W must be matrix with ndim=2"),
+        (np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1]]), "W must be a square matrix"),
+        (
+            np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1], [0, 1, 1, 0]]),
+            "W must be a symmetric matrix",
+        ),
+    ],
+)
+def test_icar_matrix_checks(W, msg):
+    with pytest.raises(ValueError, match=msg):
+        pm.ICAR.dist(W=W)
+
+
+# test logp
+
+
+def test_icar_logp():
+    W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
+
+    with pm.Model() as m:
+        RV = pm.ICAR("phi", W=W)
+
+    assert pt.isclose(
+        pm.logp(RV, np.array([0.01, -0.03, 0.02, 0.00])).eval(), np.array(4.60022238)
+    ).eval(), "logp inaccuracy"
+
+
+# test draw
+
+
+def test_icar_rng_fn():
+    W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
+
+    RV = pm.ICAR.dist(W=W)
+
+    with pytest.raises(NotImplementedError, match="Cannot sample from ICAR prior"):
+        pm.draw(RV)
+
+
 @pytest.mark.parametrize("sparse", [True, False])
 def test_car_rng_fn(sparse):
     delta = 0.05  # limit for KS p-value

From 0a6e963d5812e7816161f006042ae3ecf1769b13 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Mon, 24 Jul 2023 11:22:42 -0700
Subject: [PATCH 03/11] adjusted tests and api

---
 pymc/distributions/multivariate.py       | 44 +++++++-------
 tests/distributions/test_multivariate.py | 74 +++++++++++++++++-------
 2 files changed, 75 insertions(+), 43 deletions(-)

diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index b92cf97139..636a720887 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -86,6 +86,7 @@
     "MatrixNormal",
     "KroneckerNormal",
     "CAR",
+    "ICAR",
     "StickBreakingWeights",
 ]
 
@@ -2225,11 +2226,11 @@ class ICARRV(RandomVariable):
     dtype = "floatX"
     _print_name = ("ICAR", "\\operatorname{ICAR}")
 
-    def __call__(self, W, node1, node2, N, sigma, centering_strength, size=None, **kwargs):
-        return super().__call__(W, node1, node2, N, sigma, centering_strength, size=size, **kwargs)
+    def __call__(self, W, node1, node2, N, sigma, zero_sum_strength, size=None, **kwargs):
+        return super().__call__(W, node1, node2, N, sigma, zero_sum_strength, size=size, **kwargs)
 
     @classmethod
-    def rng_fn(cls, rng, size, W, node1, node2, N, sigma, centering_strength):
+    def rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_strength):
         raise NotImplementedError("Cannot sample from ICAR prior")
 
 
@@ -2259,8 +2260,8 @@ class ICAR(Continuous):
     distance from zero.
 
     ========  ==========================================
-    Support   :math:`x \\in \\mathbb{R}^k` ?
-    Mean      :math:`0` ?
+    Support   :math:`x \\in \\mathbb{R}^k`
+    Mean      :math:`0`
     Variance  :math:`T^{-1}` ?
     ========  ==========================================
 
@@ -2276,7 +2277,7 @@ class ICAR(Continuous):
         preferable to use a non-centered parameterization by using the default
         value and multiplying the resulting phi's by sigma. See the example below.
 
-    centering_strength : scalar, default 0.001
+    zero_sum_strength : scalar, default 0.001
         Controls how strongly to enforce the zero-sum constraint. It sets the
         standard deviation of a normal density function with mean zero.
 
@@ -2323,27 +2324,28 @@ class ICAR(Continuous):
         Mollié model in stan"
         Spatial and Spatio-temporal Epidemiology, Vol. 31, (Aug., 2019),
         pp 1-18
-
+    ..  Banerjee, S., Carlin, B., Gelfand, A. Hierarchical Modeling
+        and Analysis for Spatial Data. Second edition. CRC press. (2015)
 
     """
 
     rv_op = icar
 
     @classmethod
-    def dist(cls, W, sigma=1, centering_strength=0.001, **kwargs):
+    def dist(cls, W, sigma=1, zero_sum_strength=0.001, **kwargs):
         # check that adjacency matrix is two dimensional,
-        # square
-        # and symmetrical
+        # symmetrical
+        # and composed of 1s or 0s.
 
         if not W.ndim == 2:
             raise ValueError("W must be matrix with ndim=2")
 
-        if not W.shape[0] == W.shape[1]:
-            raise ValueError("W must be a square matrix")
-
         if not np.allclose(W.T, W):
             raise ValueError("W must be a symmetric matrix")
 
+        if np.any((W != 0) & (W != 1)):
+            raise ValueError("W must be composed of only 1s and 0s")
+
         # convert adjacency matrix to edgelist representation
 
         node1, node2 = np.where(np.tril(W) == 1)
@@ -2363,24 +2365,24 @@ def dist(cls, W, sigma=1, centering_strength=0.001, **kwargs):
 
         # check on centering_strength
 
-        centering_strength = pt.as_tensor_variable(floatX(centering_strength))
-        centering_strength = Assert("centering_strength > 0")(
-            centering_strength, pt.gt(centering_strength, 0)
+        zero_sum_strength = pt.as_tensor_variable(floatX(zero_sum_strength))
+        zero_sum_strength = Assert("centering_strength > 0")(
+            zero_sum_strength, pt.gt(zero_sum_strength, 0)
         )
 
-        return super().dist([W, node1, node2, N, sigma, centering_strength], **kwargs)
+        return super().dist([W, node1, node2, N, sigma, zero_sum_strength], **kwargs)
 
-    def moment(rv, size, W, node1, node2, N, sigma, centering_strength):
+    def moment(rv, size, W, node1, node2, N, sigma, zero_sum_strength):
         return pt.zeros(N)
 
-    def logp(value, W, node1, node2, N, sigma, centering_strength):
+    def logp(value, W, node1, node2, N, sigma, zero_sum_strength):
         pairwise_difference = (-1 / (2 * sigma**2)) * pt.sum(
             pt.square(value[node1] - value[node2])
         )
         soft_center = (
-            -0.5 * pt.pow(pt.sum(value) / (centering_strength * N), 2)
+            -0.5 * pt.pow(pt.sum(value) / (zero_sum_strength * N), 2)
             - pt.log(pt.sqrt(2.0 * np.pi))
-            - pt.log(centering_strength * N)
+            - pt.log(zero_sum_strength * N)
         )
 
         return check_parameters(pairwise_difference + soft_center, sigma > 0, msg="sigma > 0")
diff --git a/tests/distributions/test_multivariate.py b/tests/distributions/test_multivariate.py
index 0f4d236317..c9bff28b66 100644
--- a/tests/distributions/test_multivariate.py
+++ b/tests/distributions/test_multivariate.py
@@ -1060,6 +1060,18 @@ def test_car_moment(self, mu, size, expected):
             pm.CAR("x", mu=mu, W=W, alpha=alpha, tau=tau, size=size)
         assert_moment_is_expected(model, expected)
 
+    @pytest.mark.parametrize(
+        "W, expected",
+        [
+            (np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]), np.array([0, 0, 0])),
+            (np.array([[0, 1], [1, 0]]), np.array([0, 0])),
+        ],
+    )
+    def test_icar_moment(self, W, expected):
+        with pm.Model() as model:
+            RV = pm.ICAR("phi", W=W)
+        assert_moment_is_expected(model, expected)
+
     @pytest.mark.parametrize(
         "nu, mu, cov, size, expected",
         [
@@ -2084,35 +2096,53 @@ def check_draws_match_expected(self):
         ),
     ],
 )
-def test_icar_matrix_checks(W, msg):
-    with pytest.raises(ValueError, match=msg):
-        pm.ICAR.dist(W=W)
-
-
-# test logp
-
-
-def test_icar_logp():
-    W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
-
-    with pm.Model() as m:
-        RV = pm.ICAR("phi", W=W)
+class TestICAR(BaseTestDistributionRandom):
+    pymc_dist = pm.ICAR
+    pymc_dist_params = {"W": np.array([0, 1, 1], [1, 0, 1], [1, 1, 0]), "sigma": 2}
+    expected_rv_op_params = {"W": np.array([0, 1, 1], [1, 0, 1], [1, 1, 0]), "sigma": 2}
+    sizes_to_check = [None, (1), (2), (2, 2)]
+    sizes_expected = [(3,), (3), (6), (12)]
+    checks_to_run = [
+        "check_pymc_params_match_rv_op",
+        "check_rv_size",
+    ]
 
-    assert pt.isclose(
-        pm.logp(RV, np.array([0.01, -0.03, 0.02, 0.00])).eval(), np.array(4.60022238)
-    ).eval(), "logp inaccuracy"
+    @pytest.mark.parametrize(
+        "W,msg",
+        [
+            (np.array([0, 1, 0, 0]), "W must be matrix with ndim=2"),
+            (np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1]]), "W must be a symmetric matrix"),
+            (
+                np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1], [0, 1, 1, 0]]),
+                "W must be a symmetric matrix",
+            ),
+            (
+                np.array([[0, 1, 1, 0], [1, 0, 0, 0.5], [1, 0, 0, 1], [0, 0.5, 1, 0]]),
+                "W must be composed of only 1s and 0s",
+            ),
+        ],
+    )
+    def test_icar_matrix_checks(W, msg):
+        with pytest.raises(ValueError, match=msg):
+            pm.ICAR.dist(W=W)
 
+    def test_icar_logp():
+        W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
 
-# test draw
+        with pm.Model() as m:
+            RV = pm.ICAR("phi", W=W)
 
+        assert pt.isclose(
+            pm.logp(RV, np.array([0.01, -0.03, 0.02, 0.00])).eval(), np.array(4.60022238)
+        ).eval(), "logp inaccuracy"
 
-def test_icar_rng_fn():
-    W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
+    def test_icar_rng_fn():
+        W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
 
-    RV = pm.ICAR.dist(W=W)
+        RV = pm.ICAR.dist(W=W)
 
-    with pytest.raises(NotImplementedError, match="Cannot sample from ICAR prior"):
-        pm.draw(RV)
+        with pytest.raises(NotImplementedError, match="Cannot sample from ICAR prior"):
+            pm.draw(RV)
 
 
 @pytest.mark.parametrize("sparse", [True, False])

From d496689aea72111a673debcac94706ce8750e426 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Mon, 24 Jul 2023 14:47:00 -0700
Subject: [PATCH 04/11] add ICAR to lists in distributions/__init__.py

---
 pymc/distributions/__init__.py     | 2 ++
 pymc/distributions/multivariate.py | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/pymc/distributions/__init__.py b/pymc/distributions/__init__.py
index 8af6789a12..bc3d9c7863 100644
--- a/pymc/distributions/__init__.py
+++ b/pymc/distributions/__init__.py
@@ -86,6 +86,7 @@
 )
 from pymc.distributions.multivariate import (
     CAR,
+    ICAR,
     Dirichlet,
     DirichletMultinomial,
     KroneckerNormal,
@@ -198,6 +199,7 @@
     "Truncated",
     "Censored",
     "CAR",
+    "ICAR",
     "PolyaGamma",
     "HurdleGamma",
     "HurdleLogNormal",
diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index 636a720887..29056666d9 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -2238,7 +2238,7 @@ def rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_strength):
 
 
 class ICAR(Continuous):
-    """
+    r"""
     The intrinsic conditional autoregressive prior. It is primarily used to model
     covariance between neighboring areas on large datasets. It is a special case
     of the :class:`~pymc.CAR` distribution where alpha is set to 1.

From 8114a7debe6334fe02a7e8c04ebe19eabef0cb7c Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Mon, 24 Jul 2023 14:56:44 -0700
Subject: [PATCH 05/11] fix check_pymc_params_match_rv_op test

---
 tests/distributions/test_multivariate.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/distributions/test_multivariate.py b/tests/distributions/test_multivariate.py
index c9bff28b66..1b0fb97666 100644
--- a/tests/distributions/test_multivariate.py
+++ b/tests/distributions/test_multivariate.py
@@ -2098,8 +2098,8 @@ def check_draws_match_expected(self):
 )
 class TestICAR(BaseTestDistributionRandom):
     pymc_dist = pm.ICAR
-    pymc_dist_params = {"W": np.array([0, 1, 1], [1, 0, 1], [1, 1, 0]), "sigma": 2}
-    expected_rv_op_params = {"W": np.array([0, 1, 1], [1, 0, 1], [1, 1, 0]), "sigma": 2}
+    pymc_dist_params = {"W": np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]), "sigma": 2}
+    expected_rv_op_params = {"W": np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]), "sigma": 2}
     sizes_to_check = [None, (1), (2), (2, 2)]
     sizes_expected = [(3,), (3), (6), (12)]
     checks_to_run = [

From bf8be008a5ff61d1d0c2008278113d5fb1fd5b78 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Mon, 24 Jul 2023 17:26:21 -0700
Subject: [PATCH 06/11] deleted unnecessary pytest.mark.parametize

---
 tests/distributions/test_multivariate.py | 17 ++---------------
 1 file changed, 2 insertions(+), 15 deletions(-)

diff --git a/tests/distributions/test_multivariate.py b/tests/distributions/test_multivariate.py
index 1b0fb97666..ed7b80f0bb 100644
--- a/tests/distributions/test_multivariate.py
+++ b/tests/distributions/test_multivariate.py
@@ -2082,20 +2082,6 @@ def check_draws_match_expected(self):
         assert np.all(np.abs(draw(x, random_seed=rng) - np.array([0.5, 0, 2.0])) < 0.01)
 
 
-# checks on adjacency matrix
-
-
-@pytest.mark.parametrize(
-    "W,msg",
-    [
-        (np.array([0, 1, 0, 0]), "W must be matrix with ndim=2"),
-        (np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1]]), "W must be a square matrix"),
-        (
-            np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1], [0, 1, 1, 0]]),
-            "W must be a symmetric matrix",
-        ),
-    ],
-)
 class TestICAR(BaseTestDistributionRandom):
     pymc_dist = pm.ICAR
     pymc_dist_params = {"W": np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]), "sigma": 2}
@@ -2124,7 +2110,8 @@ class TestICAR(BaseTestDistributionRandom):
     )
     def test_icar_matrix_checks(W, msg):
         with pytest.raises(ValueError, match=msg):
-            pm.ICAR.dist(W=W)
+            with pm.Model():
+                pm.ICAR("phi", W=W)
 
     def test_icar_logp():
         W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])

From 299ea68d48f4024d4c861e3b597b33bcd2735976 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Tue, 25 Jul 2023 13:42:23 -0700
Subject: [PATCH 07/11] fix check for square matrix + check_pymc_match_rv_op

---
 pymc/distributions/multivariate.py       |  4 ++
 tests/distributions/test_multivariate.py | 56 ++++++++++++------------
 2 files changed, 33 insertions(+), 27 deletions(-)

diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index 29056666d9..311fbe1ddb 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -2334,12 +2334,16 @@ class ICAR(Continuous):
     @classmethod
     def dist(cls, W, sigma=1, zero_sum_strength=0.001, **kwargs):
         # check that adjacency matrix is two dimensional,
+        # square,
         # symmetrical
         # and composed of 1s or 0s.
 
         if not W.ndim == 2:
             raise ValueError("W must be matrix with ndim=2")
 
+        if not W.shape[0] == W.shape[1]:
+            raise ValueError("W must be a square matrix")
+
         if not np.allclose(W.T, W):
             raise ValueError("W must be a symmetric matrix")
 
diff --git a/tests/distributions/test_multivariate.py b/tests/distributions/test_multivariate.py
index ed7b80f0bb..6158bbd31c 100644
--- a/tests/distributions/test_multivariate.py
+++ b/tests/distributions/test_multivariate.py
@@ -2085,19 +2085,39 @@ def check_draws_match_expected(self):
 class TestICAR(BaseTestDistributionRandom):
     pymc_dist = pm.ICAR
     pymc_dist_params = {"W": np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]), "sigma": 2}
-    expected_rv_op_params = {"W": np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]), "sigma": 2}
-    sizes_to_check = [None, (1), (2), (2, 2)]
-    sizes_expected = [(3,), (3), (6), (12)]
-    checks_to_run = [
-        "check_pymc_params_match_rv_op",
-        "check_rv_size",
-    ]
+    expected_rv_op_params = {
+        "W": np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]),
+        "node1": np.array([1, 2, 2]),
+        "node2": np.array([0, 0, 1]),
+        "N": 3,
+        "sigma": 2,
+        "zero_sum_strength": 0.001,
+    }
+    checks_to_run = ["check_pymc_params_match_rv_op"]
+
+    def test_icar_logp(self):
+        W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
+
+        with pm.Model() as m:
+            RV = pm.ICAR("phi", W=W)
+
+        assert pt.isclose(
+            pm.logp(RV, np.array([0.01, -0.03, 0.02, 0.00])).eval(), np.array(4.60022238)
+        ).eval(), "logp inaccuracy"
+
+    def test_icar_rng_fn(self):
+        W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
+
+        RV = pm.ICAR.dist(W=W)
+
+        with pytest.raises(NotImplementedError, match="Cannot sample from ICAR prior"):
+            pm.draw(RV)
 
     @pytest.mark.parametrize(
         "W,msg",
         [
             (np.array([0, 1, 0, 0]), "W must be matrix with ndim=2"),
-            (np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1]]), "W must be a symmetric matrix"),
+            (np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1]]), "W must be a square matrix"),
             (
                 np.array([[0, 1, 0, 0], [1, 0, 0, 1], [1, 0, 0, 1], [0, 1, 1, 0]]),
                 "W must be a symmetric matrix",
@@ -2108,29 +2128,11 @@ class TestICAR(BaseTestDistributionRandom):
             ),
         ],
     )
-    def test_icar_matrix_checks(W, msg):
+    def test_icar_matrix_checks(self, W, msg):
         with pytest.raises(ValueError, match=msg):
             with pm.Model():
                 pm.ICAR("phi", W=W)
 
-    def test_icar_logp():
-        W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
-
-        with pm.Model() as m:
-            RV = pm.ICAR("phi", W=W)
-
-        assert pt.isclose(
-            pm.logp(RV, np.array([0.01, -0.03, 0.02, 0.00])).eval(), np.array(4.60022238)
-        ).eval(), "logp inaccuracy"
-
-    def test_icar_rng_fn():
-        W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])
-
-        RV = pm.ICAR.dist(W=W)
-
-        with pytest.raises(NotImplementedError, match="Cannot sample from ICAR prior"):
-            pm.draw(RV)
-
 
 @pytest.mark.parametrize("sparse", [True, False])
 def test_car_rng_fn(sparse):

From 6f0579a3475fb1e1b6afbbf59975ae06be548fc8 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Tue, 25 Jul 2023 15:59:07 -0700
Subject: [PATCH 08/11] fix moment test

---
 tests/distributions/test_multivariate.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/distributions/test_multivariate.py b/tests/distributions/test_multivariate.py
index 6158bbd31c..6f6159ce3b 100644
--- a/tests/distributions/test_multivariate.py
+++ b/tests/distributions/test_multivariate.py
@@ -1069,7 +1069,7 @@ def test_car_moment(self, mu, size, expected):
     )
     def test_icar_moment(self, W, expected):
         with pm.Model() as model:
-            RV = pm.ICAR("phi", W=W)
+            RV = pm.ICAR("x", W=W)
         assert_moment_is_expected(model, expected)
 
     @pytest.mark.parametrize(

From e6a7acd476e5c1b4f519d2d4d5a014fcbddc198a Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Wed, 26 Jul 2023 14:42:14 -0700
Subject: [PATCH 09/11] removed asserts, added tests on sizes

---
 pymc/distributions/multivariate.py       | 10 ----------
 tests/distributions/test_multivariate.py | 10 +++++++++-
 2 files changed, 9 insertions(+), 11 deletions(-)

diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index 311fbe1ddb..2de9ab2af9 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -2259,12 +2259,6 @@ class ICAR(Continuous):
     constraint by finding the sum of the vector $\\phi$ and penalizing based on its
     distance from zero.
 
-    ========  ==========================================
-    Support   :math:`x \\in \\mathbb{R}^k`
-    Mean      :math:`0`
-    Variance  :math:`T^{-1}` ?
-    ========  ==========================================
-
     Parameters
     ----------
     W : ndarray of int
@@ -2365,14 +2359,10 @@ def dist(cls, W, sigma=1, zero_sum_strength=0.001, **kwargs):
         # check on sigma
 
         sigma = pt.as_tensor_variable(floatX(sigma))
-        sigma = Assert("sigma > 0")(sigma, pt.gt(sigma, 0))
 
         # check on centering_strength
 
         zero_sum_strength = pt.as_tensor_variable(floatX(zero_sum_strength))
-        zero_sum_strength = Assert("centering_strength > 0")(
-            zero_sum_strength, pt.gt(zero_sum_strength, 0)
-        )
 
         return super().dist([W, node1, node2, N, sigma, zero_sum_strength], **kwargs)
 
diff --git a/tests/distributions/test_multivariate.py b/tests/distributions/test_multivariate.py
index 6f6159ce3b..a43afa1e7a 100644
--- a/tests/distributions/test_multivariate.py
+++ b/tests/distributions/test_multivariate.py
@@ -2093,7 +2093,15 @@ class TestICAR(BaseTestDistributionRandom):
         "sigma": 2,
         "zero_sum_strength": 0.001,
     }
-    checks_to_run = ["check_pymc_params_match_rv_op"]
+    checks_to_run = ["check_pymc_params_match_rv_op", "check_rv_inferred_size"]
+
+    def check_rv_inferred_size(self):
+        sizes_to_check = [None, (), 1, (1,), 5, (4, 5), (2, 4, 2)]
+        sizes_expected = [(3,), (3,), (1, 3), (1, 3), (5, 3), (4, 5, 3), (2, 4, 2, 3)]
+        for size, expected in zip(sizes_to_check, sizes_expected):
+            pymc_rv = self.pymc_dist.dist(**self.pymc_dist_params, size=size)
+            expected_symbolic = tuple(pymc_rv.shape.eval())
+            assert expected_symbolic == expected
 
     def test_icar_logp(self):
         W = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]])

From 6b466044ec5d36f9bc7c5c85c9c138d12ebb1363 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Mon, 31 Jul 2023 10:20:28 -0700
Subject: [PATCH 10/11] Implementing Bill's comments

---
 pymc/distributions/multivariate.py | 87 ++++++++++++++----------------
 1 file changed, 41 insertions(+), 46 deletions(-)

diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index 2de9ab2af9..1f0bd79832 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -2226,11 +2226,11 @@ class ICARRV(RandomVariable):
     dtype = "floatX"
     _print_name = ("ICAR", "\\operatorname{ICAR}")
 
-    def __call__(self, W, node1, node2, N, sigma, zero_sum_strength, size=None, **kwargs):
-        return super().__call__(W, node1, node2, N, sigma, zero_sum_strength, size=size, **kwargs)
+    def __call__(self, W, node1, node2, N, sigma, zero_sum_stdev, size=None, **kwargs):
+        return super().__call__(W, node1, node2, N, sigma, zero_sum_stdev, size=size, **kwargs)
 
     @classmethod
-    def rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_strength):
+    def rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_stdev):
         raise NotImplementedError("Cannot sample from ICAR prior")
 
 
@@ -2240,30 +2240,29 @@ def rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_strength):
 class ICAR(Continuous):
     r"""
     The intrinsic conditional autoregressive prior. It is primarily used to model
-    covariance between neighboring areas on large datasets. It is a special case
+    covariance between neighboring areas. It is a special case
     of the :class:`~pymc.CAR` distribution where alpha is set to 1.
 
     The log probability density function is
 
     .. math::
-        f(\\phi| W,\\sigma) =
-          -\frac{1}{2\\sigma^{2}} \\sum_{i\\sim j} (\\phi_{i} - \\phi_{j})^2 -
-          \frac{1}{2}*\frac{\\sum_{i}{\\phi_{i}}}{0.001N}^{2} - \\ln{\\sqrt{2\\pi}} -
-          \\ln{0.001N}
+        f(\phi| W,\sigma) =
+          -\frac{1}{2\sigma^{2}} \sum_{i\sim j} (\phi_{i} - \phi_{j})^2 -
+          \frac{1}{2}*\frac{\sum_{i}{\phi_{i}}}{0.001N}^{2} - \ln{\sqrt{2\\pi}} -
+          \ln{0.001N}
 
     The first term represents the spatial covariance component. Each $\\phi_{i}$ is penalized
     based on the square distance from each of its neighbors. The notation $i\\sim j$
     indicates a sum over all the neighbors of $\\phi_{i}$. The last three terms are the
     Normal log density function where the mean is zero and the standard deviation is
-    $N * 0.001$ (where N is the length of the vector $\\phi$). This component imposed the zero-sum
-    constraint by finding the sum of the vector $\\phi$ and penalizing based on its
-    distance from zero.
+    $N * 0.001$ (where N is the length of the vector $\\phi$). This component imposes
+    a zero-sum constraint by finding the sum of the vector $\\phi$ and penalizing based
+    on its distance from zero.
 
     Parameters
     ----------
     W : ndarray of int
         Symmetric adjacency matrix of 1s and 0s indicating adjacency between elements.
-        Must pass either W or both node1 and node2.
 
     sigma : scalar, default 1
         Standard deviation of the vector of phi's. Putting a prior on sigma
@@ -2271,9 +2270,11 @@ class ICAR(Continuous):
         preferable to use a non-centered parameterization by using the default
         value and multiplying the resulting phi's by sigma. See the example below.
 
-    zero_sum_strength : scalar, default 0.001
-        Controls how strongly to enforce the zero-sum constraint. It sets the
-        standard deviation of a normal density function with mean zero.
+    zero_sum_stdev : scalar, default 0.001
+        Controls how strongly to enforce the zero-sum constraint. The sum of
+        phi is normally distributed with a mean of zero and small standard deviation.
+        This parameter sets the standard deviation of a normal density function with
+        mean zero.
 
 
     Examples
@@ -2289,25 +2290,23 @@ class ICAR(Continuous):
         # 4x4 adjacency matrix
         # arranged in a square lattice
 
-        W = np.array([[0,1,0,1],
-                      [1,0,1,0],
-                      [0,1,0,1],
-                      [1,0,1,0]])
+        W = np.array([
+            [0,1,0,1],
+            [1,0,1,0],
+            [0,1,0,1],
+            [1,0,1,0]
+        ])
 
         # centered parameterization
-
         with pm.Model():
-            sigma = pm.Exponential('sigma',1)
-            phi = pm.ICAR('phi',W=W,sigma=sigma)
-
+            sigma = pm.Exponential('sigma', 1)
+            phi = pm.ICAR('phi', W=W, sigma=sigma)
             mu = phi
 
         # non-centered parameterization
-
         with pm.Model():
-            sigma = pm.Exponential('sigma',1)
-            phi = pm.ICAR('phi',W=W)
-
+            sigma = pm.Exponential('sigma', 1)
+            phi = pm.ICAR('phi', W=W)
             mu = sigma * phi
 
     References
@@ -2326,12 +2325,7 @@ class ICAR(Continuous):
     rv_op = icar
 
     @classmethod
-    def dist(cls, W, sigma=1, zero_sum_strength=0.001, **kwargs):
-        # check that adjacency matrix is two dimensional,
-        # square,
-        # symmetrical
-        # and composed of 1s or 0s.
-
+    def dist(cls, W, sigma=1, zero_sum_stdev=0.001, **kwargs):
         if not W.ndim == 2:
             raise ValueError("W must be matrix with ndim=2")
 
@@ -2345,6 +2339,12 @@ def dist(cls, W, sigma=1, zero_sum_strength=0.001, **kwargs):
             raise ValueError("W must be composed of only 1s and 0s")
 
         # convert adjacency matrix to edgelist representation
+        # An edgelist is a pair of lists.
+        # If node i and node j are connected then one list
+        # will contain i and the other will contain j at the same
+        # index value.
+        # We only use the lower triangle here because adjacency
+        # is a undirected connection.
 
         node1, node2 = np.where(np.tril(W) == 1)
 
@@ -2356,30 +2356,25 @@ def dist(cls, W, sigma=1, zero_sum_strength=0.001, **kwargs):
         N = pt.shape(W)[0]
         N = pt.as_tensor_variable(N)
 
-        # check on sigma
-
         sigma = pt.as_tensor_variable(floatX(sigma))
+        zero_sum_stdev = pt.as_tensor_variable(floatX(zero_sum_stdev))
 
-        # check on centering_strength
-
-        zero_sum_strength = pt.as_tensor_variable(floatX(zero_sum_strength))
-
-        return super().dist([W, node1, node2, N, sigma, zero_sum_strength], **kwargs)
+        return super().dist([W, node1, node2, N, sigma, zero_sum_stdev], **kwargs)
 
-    def moment(rv, size, W, node1, node2, N, sigma, zero_sum_strength):
+    def moment(rv, size, W, node1, node2, N, sigma, zero_sum_stdev):
         return pt.zeros(N)
 
-    def logp(value, W, node1, node2, N, sigma, zero_sum_strength):
+    def logp(value, W, node1, node2, N, sigma, zero_sum_stdev):
         pairwise_difference = (-1 / (2 * sigma**2)) * pt.sum(
             pt.square(value[node1] - value[node2])
         )
-        soft_center = (
-            -0.5 * pt.pow(pt.sum(value) / (zero_sum_strength * N), 2)
+        zero_sum = (
+            -0.5 * pt.pow(pt.sum(value) / (zero_sum_stdev * N), 2)
             - pt.log(pt.sqrt(2.0 * np.pi))
-            - pt.log(zero_sum_strength * N)
+            - pt.log(zero_sum_stdev * N)
         )
 
-        return check_parameters(pairwise_difference + soft_center, sigma > 0, msg="sigma > 0")
+        return check_parameters(pairwise_difference + zero_sum, sigma > 0, msg="sigma > 0")
 
 
 class StickBreakingWeightsRV(RandomVariable):

From e66e1d0877c35d46ddedf1db6ae727d0b7f44d78 Mon Sep 17 00:00:00 2001
From: Daniel Saunders <dsaunders406@gmail.com>
Date: Mon, 31 Jul 2023 13:48:38 -0700
Subject: [PATCH 11/11] add _supp_shape_from_params()

---
 pymc/distributions/multivariate.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index 1f0bd79832..ef74026840 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -33,7 +33,10 @@
 from pytensor.tensor.nlinalg import det, eigh, matrix_inverse, trace
 from pytensor.tensor.random.basic import dirichlet, multinomial, multivariate_normal
 from pytensor.tensor.random.op import RandomVariable, default_supp_shape_from_params
-from pytensor.tensor.random.utils import broadcast_params
+from pytensor.tensor.random.utils import (
+    broadcast_params,
+    supp_shape_from_ref_param_shape,
+)
 from pytensor.tensor.slinalg import Cholesky, SolveTriangular
 from pytensor.tensor.type import TensorType
 from scipy import linalg, stats
@@ -2229,6 +2232,14 @@ class ICARRV(RandomVariable):
     def __call__(self, W, node1, node2, N, sigma, zero_sum_stdev, size=None, **kwargs):
         return super().__call__(W, node1, node2, N, sigma, zero_sum_stdev, size=size, **kwargs)
 
+    def _supp_shape_from_params(self, dist_params, param_shapes=None):
+        return supp_shape_from_ref_param_shape(
+            ndim_supp=self.ndim_supp,
+            dist_params=dist_params,
+            param_shapes=param_shapes,
+            ref_param_idx=0,
+        )
+
     @classmethod
     def rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_stdev):
         raise NotImplementedError("Cannot sample from ICAR prior")