diff --git a/halotools/empirical_models/component_model_templates/scatter_models.py b/halotools/empirical_models/component_model_templates/scatter_models.py
index df0484721..39c332aec 100644
--- a/halotools/empirical_models/component_model_templates/scatter_models.py
+++ b/halotools/empirical_models/component_model_templates/scatter_models.py
@@ -101,7 +101,7 @@ def mean_scatter(self, **kwargs):
         elif 'prim_haloprop' in list(kwargs.keys()):
             mass = kwargs['prim_haloprop']
         else:
-            raise KeyError("Must pass one of the following keyword arguments to mean_occupation:\n"
+            raise KeyError("Must pass one of the following keyword arguments to mean_scatter:\n"
                 "``table`` or ``prim_haloprop``")
 
         self._update_interpol()
diff --git a/halotools/empirical_models/smhm_models/tests/test_zu_mandelbaum15.py b/halotools/empirical_models/smhm_models/tests/test_zu_mandelbaum15.py
new file mode 100644
index 000000000..be742c0a1
--- /dev/null
+++ b/halotools/empirical_models/smhm_models/tests/test_zu_mandelbaum15.py
@@ -0,0 +1,15 @@
+"""
+"""
+from __future__ import division, print_function, absolute_import, unicode_literals
+
+import numpy as np
+
+from ...smhm_models import ZuMandelbaum15SmHm
+
+
+__all__ = ('test_mc_scatter1', )
+
+
+def test_mc_scatter1():
+    model = ZuMandelbaum15SmHm()
+    sm = model.mc_stellar_mass(prim_haloprop=1e12)
diff --git a/halotools/empirical_models/smhm_models/zu_mandelbaum15.py b/halotools/empirical_models/smhm_models/zu_mandelbaum15.py
index 15f3c615d..c7fbf7494 100644
--- a/halotools/empirical_models/smhm_models/zu_mandelbaum15.py
+++ b/halotools/empirical_models/smhm_models/zu_mandelbaum15.py
@@ -5,6 +5,7 @@
 from __future__ import division, print_function, absolute_import, unicode_literals
 import numpy as np
 from warnings import warn
+from astropy.utils.misc import NumpyRNGContext
 
 from ..component_model_templates import PrimGalpropModel
 
@@ -185,3 +186,27 @@ def scatter_ln_mstar(self, halo_mass):
         eta = self.param_dict['smhm_sigma_slope']
 
         return np.where(halo_mass < m1, sigma, sigma + eta*np.log10(halo_mass/m1))
+
+    def mean_scatter(self, **kwargs):
+        if 'table' in kwargs.keys():
+            halo_mass = kwargs['table'][self.prim_haloprop_key]
+        else:
+            halo_mass = np.atleast_1d(kwargs['prim_haloprop'])
+        return np.log10(np.e)*self.scatter_ln_mstar(halo_mass)
+
+    def scatter_realization(self, **kwargs):
+        """ Monte Carlo realization of stellar mass stochasticity
+        """
+        seed = kwargs.get('seed', None)
+
+        scatter_scale = np.atleast_1d(self.mean_scatter(**kwargs))
+
+        # initialize result with zero scatter result
+        result = np.zeros(len(scatter_scale))
+
+        # only draw from a normal distribution for non-zero values of scatter
+        mask = (scatter_scale > 0.0)
+        with NumpyRNGContext(seed):
+            result[mask] = np.random.normal(loc=0, scale=scatter_scale[mask])
+
+        return result