diff --git a/docs/api.rst b/docs/api.rst
index 57e09b70d2..56f65a211a 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -44,6 +44,8 @@ Making Models from PDFs
 
    ~pdf.Model
    ~pdf._ModelConfig
+   ~pdf._MainModel
+   ~pdf._ConstraintModel
    ~mixins._ChannelSummaryMixin
    ~workspace.Workspace
    ~patchset.PatchSet
diff --git a/docs/examples/notebooks/binderexample/StatisticalAnalysis.ipynb b/docs/examples/notebooks/binderexample/StatisticalAnalysis.ipynb
index 8385164ba9..b2c0ca5063 100644
--- a/docs/examples/notebooks/binderexample/StatisticalAnalysis.ipynb
+++ b/docs/examples/notebooks/binderexample/StatisticalAnalysis.ipynb
@@ -176,7 +176,7 @@
    "outputs": [],
    "source": [
     "def get_mc_counts(pars):\n",
-    "    deltas, factors = pdf._modifications(pars)\n",
+    "    deltas, factors = pdf.modifications(pars)\n",
     "    allsum = pyhf.tensorlib.concatenate(\n",
     "        deltas + [pyhf.tensorlib.astensor(pdf.nominal_rates)]\n",
     "    )\n",
diff --git a/src/pyhf/pdf.py b/src/pyhf/pdf.py
index 2b03da9020..ca051d1652 100644
--- a/src/pyhf/pdf.py
+++ b/src/pyhf/pdf.py
@@ -650,7 +650,21 @@ def logpdf(self, maindata, pars):
         """
         return self.make_pdf(pars).log_prob(maindata)
 
-    def _modifications(self, pars):
+    def modifications(self, pars):
+        """
+        Obtain the additive and multiplicative modifications to the expected
+        event rates for the model parameters.
+
+        Args:
+            pars (:obj:`tensor`): The model parameters
+
+        Returns:
+            Tuple of additive and multiplicative modifications:
+                - deltas (:obj:`list`) is the result of an ``apply(pars)`` of combined modifiers
+                  with ``"addition"`` ``op_code``
+                - factors (:obj:`list`) is the result of ``apply(pars)`` of combined modifiers
+                  with ``"multiplication"`` ``op_code``
+        """
         deltas = list(
             filter(
                 lambda x: x is not None,
@@ -695,7 +709,7 @@ def expected_data(self, pars, return_by_sample=False):
         """
         tensorlib, _ = get_backend()
         pars = tensorlib.astensor(pars)
-        deltas, factors = self._modifications(pars)
+        deltas, factors = self.modifications(pars)
 
         allsum = tensorlib.concatenate(deltas + [self.nominal_rates])
 
@@ -826,8 +840,12 @@ def expected_auxdata(self, pars):
         pars = tensorlib.astensor(pars)
         return self.make_pdf(pars)[1].expected_data()
 
-    def _modifications(self, pars):
-        return self.main_model._modifications(pars)
+    def modifications(self, pars):
+        """
+        The modifications applied to the :class:`~pyhf.pdf._MainModel`. See
+        :func:`pyhf.pdf._MainModel.modifications` for details.
+        """
+        return self.main_model.modifications(pars)
 
     @property
     def nominal_rates(self):
diff --git a/tests/test_modifiers.py b/tests/test_modifiers.py
index 6432e75e3b..52b453d6a8 100644
--- a/tests/test_modifiers.py
+++ b/tests/test_modifiers.py
@@ -77,7 +77,7 @@ def test_staterror_holes():
 
     model = pyhf.Model(spec, poi_name="")
     assert model.config.npars == 9
-    _, factors = model._modifications(
+    _, factors = model.modifications(
         pyhf.tensorlib.astensor([2, 2.0, 1.0, 1.0, 3.0, 4.0, 1.0, 5.0, 6.0])
     )
     assert model.config.param_set("staterror_1").suggested_fixed == [
@@ -151,7 +151,7 @@ def test_shapesys_holes():
     }
 
     model = pyhf.Model(spec, poi_name="mu")
-    _, factors = model._modifications(
+    _, factors = model.modifications(
         pyhf.tensorlib.astensor([1.0, 2.0, 1.0, 1.0, 3.0, 4.0, 1.0, 1.0, 5.0])
     )
     assert (factors[1][0, 0, 0, :] == [2.0, 1.0, 1.0, 3.0, 1.0, 1.0, 1.0, 1.0]).all()
diff --git a/tests/test_pdf.py b/tests/test_pdf.py
index e680dad3f6..f12e1448ab 100644
--- a/tests/test_pdf.py
+++ b/tests/test_pdf.py
@@ -754,7 +754,7 @@ def test_lumi_np_scaling():
         [[[1.0, 1.0]]],
         [[[1.0, 1.0]]],
     ]
-    assert pdf._modifications(np.array(pars))[1][0].tolist() == [mods]
+    assert pdf.modifications(np.array(pars))[1][0].tolist() == [mods]
     assert pdf.expected_data(pars).tolist() == [120.0, 110.0, 1.0]
 
     pars[poi_slice], pars[lumi_slice] = [[1.0], [alpha_lumi]]
@@ -763,7 +763,7 @@ def test_lumi_np_scaling():
         [[[1.0, 1.0]]],
         [[[alpha_lumi, alpha_lumi]]],
     ]
-    assert pdf._modifications(np.array(pars))[1][0].tolist() == [mods]
+    assert pdf.modifications(np.array(pars))[1][0].tolist() == [mods]
     assert pytest.approx(pdf.expected_data(pars).tolist()) == [
         100 + 20.0 * alpha_lumi,
         110.0 * alpha_lumi,