diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py
index f76a98546..eab0128f5 100644
--- a/pymc/distributions/multivariate.py
+++ b/pymc/distributions/multivariate.py
@@ -67,6 +67,7 @@
 )
 from pymc.distributions.shape_utils import (
     _change_dist_size,
+    build_icar_covariance,
     change_dist_size,
     get_support_shape,
     implicit_size_from_params,
@@ -2451,7 +2452,17 @@ class ICAR(Continuous):
 
     rv_op = icar
 
-    @classmethod
+    
+    def ICAR(name, W, tau=1.0, mu=None, method="eig", **kwargs):
+        W = pt.as_tensor_variable(W)
+
+        if mu is None:
+            mu = pt.zeros(W.shape[0])
+
+        cov = build_icar_covariance(W, tau)   # private python helper that gets eig pseudo inverse
+
+        return pm.MvNormal.dist(mu=mu, cov=cov, method=method, name=name, **kwargs)
+
     def dist(cls, W, sigma=1, zero_sum_stdev=0.001, **kwargs):
         # Note: These checks are forcing W to be non-symbolic
         if not W.ndim == 2:
diff --git a/pymc/distributions/shape_utils.py b/pymc/distributions/shape_utils.py
index cdab3046b..12b8d2c9c 100644
--- a/pymc/distributions/shape_utils.py
+++ b/pymc/distributions/shape_utils.py
@@ -46,6 +46,16 @@
 from pymc.pytensorf import PotentialShapeType
 
 
+def build_icar_covariance(W, tau):
+    D = np.diag(W.sum(axis=1))
+    Q = tau*(D-W)
+    vals, vecs = np.linalg.eigh(Q)
+    tol = np.max(vals)*1e-12
+    inv_vals = np.where(vals > tol, 1/vals, 0.0)
+    cov = (vecs*inv_vals) @ vecs.T
+    cov = 0.5*(cov+cov.T)
+    return cov
+
 def to_tuple(shape):
     """Convert ints, arrays, and Nones to tuples.