Cannot sample from ICAR prior

Questions v5
1

I’m trying to build a spatially correlated Bayesian model using PyMC. By querying the API, I found that there is a pymc.ICAR function in version 5.10.4. However, when I tried to run the sample code provided in the API, it didn’t work. I want to ask if there are any issues with the function? Here is the sample code I used:

W = np.array([
    [0,1,0,1],
    [1,0,1,0],
    [0,1,0,1],
    [1,0,1,0]
])
with pm.Model():
    sigma = pm.Exponential('sigma', 1)
    phi = pm.ICAR('phi', W=W, sigma=sigma)
    mu = phi

    print(pm.draw(phi))

I encountered the following error:

NotImplementedError: Cannot sample from ICAR prior

During handling of the above exception, another exception occurred:

NotImplementedError                       Traceback (most recent call last)
/usr/local/lib/python3.10/dist-packages/pymc/distributions/multivariate.py in rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_stdev)
   2323     @classmethod
   2324     def rng_fn(cls, rng, size, W, node1, node2, N, sigma, zero_sum_stdev):
-> 2325         raise NotImplementedError("Cannot sample from ICAR prior")
   2326 
   2327 

NotImplementedError: Cannot sample from ICAR prior
Apply node that caused the error: icar_rv{1, (2, 1, 1, 0, 0, 0), floatX, True}(RandomGeneratorSharedVariable(<Generator(PCG64) at 0x7F9DDC07F140>), [], 11, [[0 1 0 1] ... [1 0 1 0]], [1 2 3 3], [0 1 0 2], 4, sigma, 0.001)
Toposort index: 1
Inputs types: [RandomGeneratorType, TensorType(int64, shape=(0,)), TensorType(int64, shape=()), TensorType(int64, shape=(4, 4)), TensorType(int64, shape=(4,)), TensorType(int64, shape=(4,)), TensorType(int64, shape=()), TensorType(float64, shape=()), TensorType(float64, shape=())]
Inputs shapes: ['No shapes', (0,), (), (4, 4), (4,), (4,), (), (), ()]
Inputs strides: ['No strides', (0,), (), (32, 8), (8,), (8,), (), (), ()]
Inputs values: [Generator(PCG64) at 0x7F9DDC07F140, array([], dtype=int64), array(11), 'not shown', array([1, 2, 3, 3]), array([0, 1, 0, 2]), array(4), array(2.22426998), array(0.001)]
Outputs clients: [['output'], ['output']]

HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.
2

Yeah that’s the expected behaviour for the ICAR. It doesn’t have a random draw function. The natural way to take random draws from ICAR would be to generate a precision matrix from the spatial structure, invert it into a covariance matrix and pass it over to the multivariate normal. Unfortunately, the precision matrix has no inverse so this is not an option.

There are two work arounds:

  1. take random draws from the CAR with a very high alpha. As alpha approaches 1, it becomes an ICAR. If alpha is just below 1, you don’t get the matrix inverse problem and you can still sample from the multivariate normal. So if you just want to get a feel for how this distribution behaves, this is a nice trick.

     N = 4
 adj_matrix = np.array([[0,1,0,0],
                       [1,0,1,0],
                       [0,1,0,1],
                       [0,0,1,0]])
 
 phi = pm.CAR.dist(mu=np.zeros(N), tau=1, alpha=0.99, W=adj_matrix)

 pm.draw(phi)

  2. You can run NUTS or any other MCMC to sample from the distribution. It has a well-defined log-likelihood so MCMC strategies can sample from it, even if there is no natural random draw function.

     W = np.array([
     [0,1,0,1],
     [1,0,1,0],
     [0,1,0,1],
     [1,0,1,0]
 ])
 with pm.Model():
     sigma = pm.Exponential('sigma', 1)
     phi = pm.ICAR('phi', W=W, sigma=sigma)
     trace = pm.sample()
1 Like