diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst
index c74769be6..454409f1f 100644
--- a/docs/source/pythonapi/data.rst
+++ b/docs/source/pythonapi/data.rst
@@ -61,6 +61,7 @@ Core Functions
 
     atomic_mass
     atomic_weight
+    combine_distributions
     decay_constant
     dose_coefficients
     gnd_name
diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 51a594bc8..2c2505bf5 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -9,7 +9,9 @@ from uncertainties import ufloat, UFloat
 
 import openmc.checkvalue as cv
 from openmc.mixin import EqualityMixin
+from openmc.stats import Discrete, Tabular, combine_distributions
 from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
+from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
 
 
@@ -314,6 +316,12 @@ class Decay(EqualityMixin):
         'excited_state', 'mass', 'stable', 'spin', and 'parity'.
     spectra : dict
         Resulting radiation spectra for each radiation type.
+    sources : dict
+        Radioactive decay source distributions represented as a dictionary
+        mapping particle types (e.g., 'photon') to instances of
+        :class:`openmc.stats.Univariate`.
+
+        .. versionadded:: 0.13.1
 
     """
     def __init__(self, ev_or_filename):
@@ -329,6 +337,7 @@ class Decay(EqualityMixin):
         self.modes = []
         self.spectra = {}
         self.average_energies = {}
+        self._sources = None
 
         # Get head record
         items = get_head_record(file_obj)
@@ -495,3 +504,69 @@ class Decay(EqualityMixin):
 
         """
         return cls(ev_or_filename)
+
+    @property
+    def sources(self):
+        """Radioactive decay source distributions"""
+        # If property has been computed already, return it
+        # TODO: Replace with functools.cached_property when support is Python 3.9+
+        if self._sources is not None:
+            return self._sources
+
+        sources = {}
+        name = self.nuclide['name']
+        decay_constant = self.decay_constant.n
+        for particle, spectra in self.spectra.items():
+            # Set particle type based on 'particle' above
+            particle_type = {
+                'gamma': 'photon',
+                'beta-': 'electron',
+                'ec/beta+': 'positron',
+                'alpha': 'alpha',
+                'n': 'neutron',
+                'sf': 'fragment',
+                'p': 'proton',
+                'e-': 'electron',
+                'xray': 'photon',
+                'anti-neutrino': 'anti-neutrino',
+                'neutrino': 'neutrino',
+            }[particle]
+
+            if particle_type not in sources:
+                sources[particle_type] = []
+
+            # Create distribution for discrete
+            if spectra['continuous_flag'] in ('discrete', 'both'):
+                energies = []
+                intensities = []
+                for discrete_data in spectra['discrete']:
+                    energies.append(discrete_data['energy'].n)
+                    intensities.append(discrete_data['intensity'].n)
+                energies = np.array(energies)
+                intensity = spectra['discrete_normalization'].n
+                rates = decay_constant * intensity * np.array(intensities)
+                dist_discrete = Discrete(energies, rates)
+                sources[particle_type].append(dist_discrete)
+
+            # Create distribution for continuous
+            if spectra['continuous_flag'] in ('continuous', 'both'):
+                f = spectra['continuous']['probability']
+                if len(f.interpolation) > 1:
+                    raise NotImplementedError("Multiple interpolation regions: {name}, {particle}")
+                interpolation = INTERPOLATION_SCHEME[f.interpolation[0]]
+                if interpolation not in ('histogram', 'linear-linear'):
+                    raise NotImplementedError("Continuous spectra with {interpolation} interpolation ({name}, {particle}) not supported")
+
+                intensity = spectra['continuous_normalization'].n
+                rates = decay_constant * intensity * f.y
+                dist_continuous = Tabular(f.x, rates, interpolation)
+                sources[particle_type].append(dist_continuous)
+
+        # Combine discrete distributions
+        merged_sources = {}
+        for particle_type, dist_list in sources.items():
+            merged_sources[particle_type] = combine_distributions(
+                dist_list, [1.0]*len(dist_list))
+
+        self._sources = merged_sources
+        return self._sources
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 5528bdb30..58b7c6172 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -1,6 +1,7 @@
 from abc import ABC, abstractmethod
 from collections import defaultdict
 from collections.abc import Iterable
+from copy import deepcopy
 from numbers import Real
 from xml.etree import ElementTree as ET
 
@@ -95,9 +96,9 @@ class Discrete(Univariate):
 
     Attributes
     ----------
-    x : Iterable of float
+    x : numpy.ndarray
         Values of the random variable
-    p : Iterable of float
+    p : numpy.ndarray
         Discrete probability for each value
 
     """
@@ -122,7 +123,7 @@ class Discrete(Univariate):
         if isinstance(x, Real):
             x = [x]
         cv.check_type('discrete values', x, Iterable, Real)
-        self._x = x
+        self._x = np.array(x, dtype=float)
 
     @p.setter
     def p(self, p):
@@ -131,14 +132,15 @@ class Discrete(Univariate):
         cv.check_type('discrete probabilities', p, Iterable, Real)
         for pk in p:
             cv.check_greater_than('discrete probability', pk, 0.0, True)
-        self._p = p
+        self._p = np.array(p, dtype=float)
 
     def cdf(self):
         return np.insert(np.cumsum(self.p), 0, 0.0)
 
     def sample(self, n_samples=1, seed=None):
         np.random.seed(seed)
-        return np.random.choice(self.x, n_samples, p=self.p)
+        p = self.p / self.p.sum()
+        return np.random.choice(self.x, n_samples, p=p)
 
     def normalize(self):
         """Normalize the probabilities stored on the distribution"""
@@ -220,6 +222,18 @@ class Discrete(Univariate):
         p_arr = np.array([p_merged[x] for x in x_arr])
         return cls(x_arr, p_arr)
 
+    def integral(self):
+        """Return integral of distribution
+
+        .. versionadded:: 0.13.1
+
+        Returns
+        -------
+        float
+            Integral of discrete distribution
+        """
+        return np.sum(self.p)
+
 class Uniform(Univariate):
     """Distribution with constant probability over a finite interval [a,b]
 
@@ -825,9 +839,9 @@ class Tabular(Univariate):
 
     Attributes
     ----------
-    x : Iterable of float
+    x : numpy.ndarray
         Tabulated values of the random variable
-    p : Iterable of float
+    p : numpy.ndarray
         Tabulated probabilities
     interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'}, optional
         Indicate whether the density function is constant between tabulated
@@ -860,7 +874,7 @@ class Tabular(Univariate):
     @x.setter
     def x(self, x):
         cv.check_type('tabulated values', x, Iterable, Real)
-        self._x = x
+        self._x = np.array(x, dtype=float)
 
     @p.setter
     def p(self, p):
@@ -868,7 +882,7 @@ class Tabular(Univariate):
         if not self._ignore_negative:
             for pk in p:
                 cv.check_greater_than('tabulated probability', pk, 0.0, True)
-        self._p = p
+        self._p = np.array(p, dtype=float)
 
     @interpolation.setter
     def interpolation(self, interpolation):
@@ -881,8 +895,8 @@ class Tabular(Univariate):
                                       'distributions using histogram or '
                                       'linear-linear interpolation')
         c = np.zeros_like(self.x)
-        x = np.asarray(self.x)
-        p = np.asarray(self.p)
+        x = self.x
+        p = self.p
 
         if self.interpolation == 'histogram':
             c[1:] = p[:-1] * np.diff(x)
@@ -922,7 +936,7 @@ class Tabular(Univariate):
 
     def normalize(self):
         """Normalize the probabilities stored on the distribution"""
-        self.p = np.asarray(self.p) / self.cdf().max()
+        self.p /= self.cdf().max()
 
     def sample(self, n_samples=1, seed=None):
         if not self.interpolation in ('histogram', 'linear-linear'):
@@ -931,10 +945,11 @@ class Tabular(Univariate):
                                       'linear-linear interpolation')
         np.random.seed(seed)
         xi = np.random.rand(n_samples)
-        cdf = self.cdf()
-        cdf /= cdf.max()
+
         # always use normalized probabilities when sampling
+        cdf = self.cdf()
         p = self.p / cdf.max()
+        cdf /= cdf.max()
 
         # get CDF bins that are above the
         # sampled values
@@ -949,7 +964,7 @@ class Tabular(Univariate):
         # the random number is less than the next cdf
         # entry
         x_i = self.x[cdf_idx]
-        p_i = self.p[cdf_idx]
+        p_i = p[cdf_idx]
 
         if self.interpolation == 'histogram':
             # mask where probability is greater than zero
@@ -967,7 +982,7 @@ class Tabular(Univariate):
             # get variable and probability values for the
             # next entry
             x_i1 = self.x[cdf_idx + 1]
-            p_i1 = self.p[cdf_idx + 1]
+            p_i1 = p[cdf_idx + 1]
             # compute slope between entries
             m = (p_i1 - p_i) / (x_i1 - x_i)
             # set values for zero slope
@@ -1027,6 +1042,24 @@ class Tabular(Univariate):
         p = params[len(params)//2:]
         return cls(x, p, interpolation)
 
+    def integral(self):
+        """Return integral of distribution
+
+        .. versionadded: 0.13.1
+
+        Returns
+        -------
+        float
+            Integral of tabular distrbution
+        """
+        if self.interpolation == 'histogram':
+            return np.sum(np.diff(self.x) * self.p[:-1])
+        elif self.interpolation == 'linear-linear':
+            return np.trapz(self.p, self.x)
+        else:
+            raise NotImplementedError(
+                f'integral() not supported for {self.inteprolation} interpolation')
+
 
 class Legendre(Univariate):
     r"""Probability density given by a Legendre polynomial expansion
@@ -1135,9 +1168,10 @@ class Mixture(Univariate):
 
     def sample(self, n_samples=1, seed=None):
         np.random.seed(seed)
-        idx = np.random.choice(self.distribution, n_samples, p=self.probability)
+        idx = np.random.choice(range(len(self.distribution)),
+                               n_samples, p=self.probability)
 
-        out = np.zeros_like(idx)
+        out = np.empty_like(idx, dtype=float)
         for i in np.unique(idx):
             n_dist_samples = np.count_nonzero(idx == i)
             samples = self.distribution[i].sample(n_dist_samples)
@@ -1199,3 +1233,68 @@ class Mixture(Univariate):
             distribution.append(Univariate.from_xml_element(pair.find("dist")))
 
         return cls(probability, distribution)
+
+    def integral(self):
+        """Return integral of the distribution
+
+        .. versionadded:: 0.13.1
+
+        Returns
+        -------
+        float
+            Integral of the distribution
+        """
+        return sum([
+            p*dist.integral()
+            for p, dist in zip(self.probability, self.distribution)
+        ])
+
+
+def combine_distributions(dists, probs):
+    """Combine distributions with specified probabilities
+
+    This function can be used to combine multiple instances of
+    :class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular`. Multiple
+    discrete distributions are merged into a single distribution and the
+    remainder of the distributions are put into a :class:`~openmc.stats.Mixture`
+    distribution.
+
+    .. versionadded:: 0.13.1
+
+    Parameters
+    ----------
+    dists : iterable of openmc.stats.Univariate
+        Distributions to combine
+    probs : iterable of float
+        Probability (or intensity) of each distribution
+
+    """
+    # Get copy of distribution list so as not to modify the argument
+    dist_list = deepcopy(dists)
+
+    # Get list of discrete/continuous distribution indices
+    discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
+    cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)]
+
+    # Apply probabilites to continuous distributions
+    for i in cont_index:
+        dist = dist_list[i]
+        dist.p *= probs[i]
+
+    if discrete_index:
+        # Create combined discrete distribution
+        dist_discrete = [dist_list[i] for i in discrete_index]
+        discrete_probs = [probs[i] for i in discrete_index]
+        combined_dist = Discrete.merge(dist_discrete, discrete_probs)
+
+        # Replace multiple discrete distributions with merged
+        for idx in reversed(discrete_index):
+            dist_list.pop(idx)
+        dist_list.append(combined_dist)
+
+    # Combine discrete and continuous if present
+    if len(dist_list) > 1:
+        probs = [d.integral() for d in dist_list]
+        dist_list[:] = [Mixture(probs, dist_list.copy())]
+
+    return dist_list[0]
diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py
index 06b8e6bed..f0bda1bd9 100644
--- a/tests/unit_tests/test_data_decay.py
+++ b/tests/unit_tests/test_data_decay.py
@@ -23,6 +23,14 @@ def nb90():
     return openmc.data.Decay.from_endf(filename)
 
 
+@pytest.fixture(scope='module')
+def ba137m():
+    """Ba137_m1 decay data."""
+    endf_data = os.environ['OPENMC_ENDF_DATA']
+    filename = os.path.join(endf_data, 'decay', 'dec-056_Ba_137m1.endf')
+    return openmc.data.Decay.from_endf(filename)
+
+
 @pytest.fixture(scope='module')
 def u235_yields():
     """U235 fission product yield data."""
@@ -48,6 +56,7 @@ def test_nb90_halflife(nb90):
     ufloat_close(nb90.decay_constant, log(2.)/nb90.half_life)
     ufloat_close(nb90.decay_energy, ufloat(2265527.5, 25159.400474401213))
 
+
 def test_nb90_nuclide(nb90):
     assert nb90.nuclide['atomic_number'] == 41
     assert nb90.nuclide['mass_number'] == 90
@@ -91,3 +100,29 @@ def test_fpy(u235_yields):
     assert len(u235_yields.independent) == 3
     thermal = u235_yields.independent[0]
     ufloat_close(thermal['I135'], ufloat(0.0292737, 0.000819663))
+
+
+def test_sources(ba137m, nb90):
+    # Running .sources twice should give same objects
+    sources = ba137m.sources
+    sources2 = ba137m.sources
+    for key in sources:
+        assert sources[key] is sources2[key]
+
+    # Each source should be a univariate distribution
+    for dist in sources.values():
+        assert isinstance(dist, openmc.stats.Univariate)
+
+    # Check for presence of 662 keV gamma ray in decay of Ba137m
+    gamma_source = ba137m.sources['photon']
+    assert isinstance(gamma_source, openmc.stats.Discrete)
+    b = np.isclose(gamma_source.x, 661657.)
+    assert np.count_nonzero(b) == 1
+
+    # Check value of decay/s/atom
+    idx = np.flatnonzero(b)[0]
+    assert gamma_source.p[idx] == pytest.approx(0.004069614)
+
+    # Nb90 decays by β+ and should emit positrons, electrons, and photons
+    sources = nb90.sources
+    assert len(set(sources.keys()) ^ {'positron', 'electron', 'photon'}) == 0
diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index b57f9578a..1cfcf00c3 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -6,6 +6,12 @@ import openmc
 import openmc.stats
 
 
+def assert_sample_mean(samples, expected_mean):
+    std_dev = samples.std() / np.sqrt(samples.size)
+    assert np.abs(expected_mean - samples.mean()) < 3*std_dev
+
+
+
 def test_discrete():
     x = [0.0, 1.0, 10.0]
     p = [0.3, 0.2, 0.5]
@@ -13,8 +19,8 @@ def test_discrete():
     elem = d.to_xml_element('distribution')
 
     d = openmc.stats.Discrete.from_xml_element(elem)
-    assert d.x == x
-    assert d.p == p
+    np.testing.assert_array_equal(d.x, x)
+    np.testing.assert_array_equal(d.p, p)
     assert len(d) == len(x)
 
     d = openmc.stats.Univariate.from_xml_element(elem)
@@ -33,13 +39,11 @@ def test_discrete():
 
     d3 = openmc.stats.Discrete(vals, probs)
 
-    # sample discrete distribution
+    # sample discrete distribution and check that the mean of the samples is
+    # within 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d3.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_merge_discrete():
@@ -76,17 +80,16 @@ def test_uniform():
     assert len(d) == 2
 
     t = d.to_tabular()
-    assert t.x == [a, b]
-    assert t.p == [1/(b-a), 1/(b-a)]
+    np.testing.assert_array_equal(t.x, [a, b])
+    np.testing.assert_array_equal(t.p, [1/(b-a), 1/(b-a)])
     assert t.interpolation == 'histogram'
 
+    # Sample distribution and check that the mean of the samples is within 3
+    # std. dev. of the expected mean
     exp_mean = 0.5 * (a + b)
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_powerlaw():
@@ -102,13 +105,11 @@ def test_powerlaw():
 
     exp_mean = 100.0 * (n+1) / (n+2)
 
-    # sample power law distribution
+    # sample power law distribution and check that the mean of the samples is
+    # within 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_maxwell():
@@ -122,21 +123,15 @@ def test_maxwell():
 
     exp_mean = 3/2 * theta
 
-    # sample maxwell distribution
+    # sample maxwell distribution and check that the mean of the samples is
+    # within 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
     # A second sample with a different seed
     samples_2 = d.sample(n_samples, seed=200)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples_2.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples_2.mean()) < 3*std_dev
-
+    assert_sample_mean(samples_2, exp_mean)
     assert samples_2.mean() != samples.mean()
 
 
@@ -156,13 +151,11 @@ def test_watt():
     # https://doi.org/10.1016/j.physletb.2003.09.048
     exp_mean = 3/2 * a + a**2 * b / 4
 
-    # sample Watt distribution
+    # sample Watt distribution and check that the mean of the samples is within
+    # 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_tabular():
@@ -228,6 +221,11 @@ def test_mixture():
     assert mix.distribution == [d1, d2]
     assert len(mix) == 4
 
+    # Sample and make sure sample mean is close to expected mean
+    n_samples = 1_000_000
+    samples = mix.sample(n_samples)
+    assert_sample_mean(samples, (2.5 + 5.0)/2)
+
     elem = mix.to_xml_element('distribution')
 
     d = openmc.stats.Mixture.from_xml_element(elem)
@@ -396,3 +394,47 @@ def test_muir():
     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
+
+
+def test_combine_distributions():
+    # Combine two discrete (same data as in test_merge_discrete)
+    x1 = [0.0, 1.0, 10.0]
+    p1 = [0.3, 0.2, 0.5]
+    d1 = openmc.stats.Discrete(x1, p1)
+    x2 = [0.5, 1.0, 5.0]
+    p2 = [0.4, 0.5, 0.1]
+    d2 = openmc.stats.Discrete(x2, p2)
+
+    # Merged distribution should have x values sorted and probabilities
+    # appropriately combined. Duplicate x values should appear once.
+    merged = openmc.stats.combine_distributions([d1, d2], [0.6, 0.4])
+    assert isinstance(merged, openmc.stats.Discrete)
+    assert merged.x == pytest.approx([0.0, 0.5, 1.0, 5.0, 10.0])
+    assert merged.p == pytest.approx(
+        [0.6*0.3, 0.4*0.4, 0.6*0.2 + 0.4*0.5, 0.4*0.1, 0.6*0.5])
+
+    # Probabilities add up but are not normalized
+    d1 = openmc.stats.Discrete([3.0], [1.0])
+    triple = openmc.stats.combine_distributions([d1, d1, d1], [1.0, 2.0, 3.0])
+    assert triple.x == pytest.approx([3.0])
+    assert triple.p == pytest.approx([6.0])
+
+    # Combine discrete and tabular
+    t1 = openmc.stats.Tabular(x2, p2)
+    mixed = openmc.stats.combine_distributions([d1, t1], [0.5, 0.5])
+    assert isinstance(mixed, openmc.stats.Mixture)
+    assert len(mixed.distribution) == 2
+    assert len(mixed.probability) == 2
+
+    # Combine 1 discrete and 2 tabular -- the tabular distributions should
+    # combine to produce a uniform distribution with mean 0.5. The combined
+    # distribution should have a mean of 0.25.
+    t1 = openmc.stats.Tabular([0., 1.], [2.0, 0.0])
+    t2 = openmc.stats.Tabular([0., 1.], [0.0, 2.0])
+    d1 = openmc.stats.Discrete([0.0], [1.0])
+    combined = openmc.stats.combine_distributions([t1, t2, d1], [0.25, 0.25, 0.5])
+
+    # Sample the combined distribution and make sure the sample mean is within
+    # uncertainty of the expected value
+    samples = combined.sample(10)
+    assert_sample_mean(samples, 0.25)