diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 0e6b4ae61d..4f2fbf66c2 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -764,11 +764,14 @@ class Muir(Univariate):
         cv.check_greater_than('Muir kt', kt, 0.0)
         self._kt = kt
 
+    @property
+    def std_dev(self):
+        return np.sqrt(4.*self.e0*self.kt/self.m_rat)
+
     def sample(self, n_samples=1, seed=None):
         # https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
         np.random.seed(seed)
-        sigma = np.sqrt(4.*self.e0*self.kt/self.m_rat)
-        return np.random.normal(self.e0, sigma, n_samples)
+        return np.random.normal(self.e0, self.std_dev, n_samples)
 
     def to_xml_element(self, element_name):
         """Return XML representation of the Watt distribution
diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index 7bf67c1fa2..532936e422 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -292,3 +292,14 @@ def test_muir():
     assert d.m_rat == pytest.approx(mass)
     assert d.kt == pytest.approx(temp)
     assert len(d) == 3
+
+    # sample normal distribution
+    n_samples = 10000
+    samples = d.sample(n_samples, seed=100)
+    samples = np.abs(samples - mean)
+    within_1_sigma = np.count_nonzero(samples < d.std_dev)
+    assert within_1_sigma / n_samples >= 0.68
+    within_2_sigma = np.count_nonzero(samples < 2*d.std_dev)
+    assert within_2_sigma / n_samples >= 0.95
+    within_3_sigma = np.count_nonzero(samples < 3*d.std_dev)
+    assert within_3_sigma / n_samples >= 0.99