diff --git a/pyest/gm/gm.py b/pyest/gm/gm.py
index e485502..ad97d36 100755
--- a/pyest/gm/gm.py
+++ b/pyest/gm/gm.py
@@ -31,6 +31,7 @@
     'v_eval_mvnpdfchol',
     'integral_gauss_product',
     'integral_gauss_product_chol',
+    'integral_squared_gm',
     'marginal_2d',
     'marginal_nd',
     'comp_bounds',
@@ -463,6 +464,24 @@ def integral_gauss_product(m1, P1, m2, P2, allow_singular=False):
     return multivariate_normal.pdf(m1, m2, P1 + P2, allow_singular=allow_singular)
 
 
+def integral_squared_gm(p):
+    ''' compute integral of squared Gaussian mixture
+
+    Parameters
+    ----------
+    p : GaussianMixture
+        Gaussian mixture
+
+    Returns
+    -------
+    integral : float
+        integral of the squared Gaussian mixture
+    '''
+    return np.sum([
+        wi*wj*eval_mvnpdf(mi, mj, Pi + Pj) for (wi, mi, Pi) in p for (wj, mj, Pj) in p
+    ])
+
+
 def marginal_2d(m, P, dimensions=[0, 1]):
     """ compute 2D marginal of GM"""
     return marginal_nd(m, P, dimensions)
diff --git a/pyest/metrics.py b/pyest/metrics.py
index 6b2e1a5..f2837ae 100644
--- a/pyest/metrics.py
+++ b/pyest/metrics.py
@@ -1,6 +1,7 @@
 import numpy as np
 import pyest.gm as pygm
 from scipy.linalg import solve_triangular
+from scipy.integrate import dblquad
 
 
 def l2_dist(p1, p2):
@@ -87,4 +88,106 @@ def madem(m, S, m_ref):
     """
     return np.linalg.norm(
         solve_triangular(S, m - m_ref, lower=True)
-    )
\ No newline at end of file
+    )
+
+
+def integral_squared_error_2d(p1, p2, a, b, c, d, epsabs=1.49e-2, epsrel=1.49e-2):
+    ''' compute integral squared error between two 2D densities
+
+    Parameters
+    ----------
+    p1: callable
+        first density, p1([x,y])
+    p2: callable
+        second density, p2([x,y])
+    a, b : float
+        The limits of integration in x: a<b
+    c, d : float
+        The limits of integration in y: c<d
+    epsabs : float, optional
+        absolute error tolerance for numerical integration
+    epsrel : float, optional
+        relative error tolerance for numerical integration
+
+
+    Returns
+    -------
+    ise: foat
+        kld
+    int_error:
+        numerical integration estimated error
+
+    See Also
+    --------
+    normalized_integral_squared_error_2d : compute normalized integral squared
+        error between two 2D densities
+    l2_dist : compute L2 distance (ISE) between two Gaussian mixtures
+
+    Notes
+    -----
+    This function is intended for use with generic callable densities and
+    makes no assumptions about the form of the densities. If both p1 and p2
+    are Gaussian mixtures, use l2_dist instead, which is exact and more
+    efficient.
+    '''
+    integrand_fun = lambda y,x: (p1([x,y])-p2([x,y]))**2
+    return dblquad(integrand_fun, a, b, c, d, epsabs=epsabs, epsrel=epsrel)
+
+
+def normalized_integral_squared_error_2d(p1, p2, a, b, c, d, epsabs=1.49e-2, epsrel=1.49e-2):
+    ''' compute normalized integral squared error between two 2D, densities
+
+    Parameters
+    ----------
+    p1: callable
+        first density
+    p2: callable
+        second density
+    a, b : float
+        The limits of integration in x: a<b
+    c, d : float
+        The limits of integration in y: c<d
+    epsabs : float, optional
+        absolute error tolerance for numerical integration
+    epsrel : float, optional
+        relative error tolerance for numerical integration
+
+    Returns
+    -------
+    nise: float
+        normalized integral squared error
+    ise:
+        integral squared error
+    err:
+        numerical integration estimated error in ISE computation
+
+    See Also
+    --------
+    integral_squared_error_2d : compute integral squared error between two
+        2D densities
+    l2_dist : compute L2 distance (ISE) between two Gaussian mixtures
+
+    Notes
+    -----
+    This function is intended for use with generic callable densities and
+    makes no assumptions about the form of the densities. If both p1 and p2
+    are Gaussian mixtures, use metrics.l2_dist and gm.integral_squared_gm
+    instead for the numerator and denominator terms separately, which is
+    exact and more efficient.
+
+    '''
+    ise, err = integral_squared_error_2d(p1, p2, a, b, c, d, epsabs=epsabs, epsrel=epsrel)
+    # if p1 is a GaussianMixtureRv, use l2_dist
+    if isinstance(p1, pygm.GaussianMixture):
+        int_p1_sq = pygm.integral_squared_gm(p1)
+    else:
+        p1_sq_integrand_fun = lambda y,x: p1([x,y])**2
+        int_p1_sq = dblquad(p1_sq_integrand_fun, a, b, c, d, epsabs=epsabs, epsrel=epsrel)[0]
+    if isinstance(p2, pygm.GaussianMixture):
+        int_p2_sq = pygm.integral_squared_gm(p2)
+    else:
+        p2_sq_integrand_fun = lambda y,x: p2([x,y])**2
+        int_p2_sq = dblquad(p2_sq_integrand_fun, a, b, c, d, epsabs=epsabs, epsrel=epsrel)[0]
+    nise = ise/(int_p1_sq + int_p2_sq)
+
+    return nise, ise, err
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
index 8d2cd64..1296aca 100644
--- a/tests/test_metrics.py
+++ b/tests/test_metrics.py
@@ -22,6 +22,48 @@ def test_l2_dist():
     npt.assert_approx_equal(gm.l2_dist(gm1, gm3), des_l2, significant=9)
 
 
+def test_integral_squared_error_2d():
+    """ test computation of integral squared error between 2D densities """
+
+    # test identical mixtures
+    p1 = gm.defaults.default_gm().marginal_2d([0, 1])
+    p2 = gm.defaults.default_gm().marginal_2d([0, 1])
+    lb, ub = gm.bounds(p1.m, p1.P, sigma_mult=3)
+
+    ise, int_err = metrics.integral_squared_error_2d(p1, p2, lb[0], ub[0], lb[1], ub[1])
+    assert abs(ise) < 1e-10
+
+    # now, create two mixtures with disjoint supports. The ISE in this case
+    # should trend toward the sum of the individual integrals of the squared densities.
+    shift = ub - lb
+    p2 = gm.defaults.default_gm(mean_shift=np.array([shift[0], shift[1], 0, 0])).marginal_2d([0, 1])
+    ise, int_err = metrics.integral_squared_error_2d(p1, p2, lb[0], ub[0] + shift[0], lb[1], ub[1] + shift[1])
+    int_p1_sq = gm.integral_squared_gm(p1)
+    int_p2_sq = gm.integral_squared_gm(p2)
+    des_ise = int_p1_sq + int_p2_sq
+    assert abs(ise - des_ise) < 1e-5
+
+
+def test_normalized_integral_squared_error_2d():
+    """ test computation of NISE between 2D densities """
+
+    # test identical mixtures
+    p1 = gm.defaults.default_gm().marginal_2d([0, 1])
+    p2 = gm.defaults.default_gm().marginal_2d([0, 1])
+    lb, ub = gm.bounds(p1.m, p1.P, sigma_mult=3)
+
+    nise, ise, int_err = metrics.normalized_integral_squared_error_2d(p1, p2, lb[0], ub[0], lb[1], ub[1])
+    assert abs(nise) < 1e-10
+
+    # now, create two mixtures with disjoint supports. The NISE in this case
+    # should trend toward 1
+    shift = ub - lb
+    p2 = gm.defaults.default_gm(mean_shift=np.array([shift[0], shift[1], 0, 0])).marginal_2d([0, 1])
+    nise, ise, int_err = metrics.normalized_integral_squared_error_2d(p1, p2, lb[0], ub[0] + shift[0], lb[1], ub[1] + shift[1])
+    des_nise = 1
+    assert abs(nise - des_nise) < 1e-4
+
+
 
 if __name__ == '__main__':
     pytest.main([__file__])