diff --git a/docs/source/uml/classes_mmm.png b/docs/source/uml/classes_mmm.png
index 5bf93583..c2cb8683 100644
Binary files a/docs/source/uml/classes_mmm.png and b/docs/source/uml/classes_mmm.png differ
diff --git a/pymc_marketing/mmm/__init__.py b/pymc_marketing/mmm/__init__.py
index 1c452e70..c636272a 100644
--- a/pymc_marketing/mmm/__init__.py
+++ b/pymc_marketing/mmm/__init__.py
@@ -25,6 +25,7 @@
 )
 from pymc_marketing.mmm.components.saturation import (
     HillSaturation,
+    HillSaturationSigmoid,
     InverseScaledLogisticSaturation,
     LogisticSaturation,
     MichaelisMentenSaturation,
@@ -50,6 +51,7 @@
     "DelayedSaturatedMMM",
     "GeometricAdstock",
     "HillSaturation",
+    "HillSaturationSigmoid",
     "LogisticSaturation",
     "InverseScaledLogisticSaturation",
     "MMM",
diff --git a/pymc_marketing/mmm/components/saturation.py b/pymc_marketing/mmm/components/saturation.py
index b510dddb..7cba45aa 100644
--- a/pymc_marketing/mmm/components/saturation.py
+++ b/pymc_marketing/mmm/components/saturation.py
@@ -77,7 +77,8 @@ def function(self, x, b):
 
 from pymc_marketing.mmm.components.base import Transformation
 from pymc_marketing.mmm.transformers import (
-    hill_saturation,
+    hill_function,
+    hill_saturation_sigmoid,
     inverse_scaled_logistic_saturation,
     logistic_saturation,
     michaelis_menten,
@@ -343,7 +344,7 @@ class MichaelisMentenSaturation(SaturationTransformation):
 class HillSaturation(SaturationTransformation):
     """Wrapper around Hill saturation function.
 
-    For more information, see :func:`pymc_marketing.mmm.transformers.hill_saturation`.
+    For more information, see :func:`pymc_marketing.mmm.transformers.hill_function`.
 
     .. plot::
         :context: close-figs
@@ -364,7 +365,41 @@ class HillSaturation(SaturationTransformation):
 
     lookup_name = "hill"
 
-    function = hill_saturation
+    def function(self, x, slope, kappa, beta):
+        return beta * hill_function(x, slope, kappa)
+
+    default_priors = {
+        "slope": Prior("HalfNormal", sigma=1.5),
+        "kappa": Prior("HalfNormal", sigma=1.5),
+        "beta": Prior("HalfNormal", sigma=1.5),
+    }
+
+
+class HillSaturationSigmoid(SaturationTransformation):
+    """Wrapper around Hill saturation sigmoid function.
+
+    For more information, see :func:`pymc_marketing.mmm.transformers.hill_saturation_sigmoid`.
+
+    .. plot::
+        :context: close-figs
+
+        import matplotlib.pyplot as plt
+        import numpy as np
+        from pymc_marketing.mmm import HillSaturationSigmoid
+
+        rng = np.random.default_rng(0)
+
+        adstock = HillSaturationSigmoid()
+        prior = adstock.sample_prior(random_seed=rng)
+        curve = adstock.sample_curve(prior)
+        adstock.plot_curve(curve, sample_kwargs={"rng": rng})
+        plt.show()
+
+    """
+
+    lookup_name = "hill_sigmoid"
+
+    function = hill_saturation_sigmoid
 
     default_priors = {
         "sigma": Prior("HalfNormal", sigma=1.5),
@@ -415,6 +450,7 @@ def function(self, x, alpha, beta):
         TanhSaturationBaselined,
         MichaelisMentenSaturation,
         HillSaturation,
+        HillSaturationSigmoid,
         RootSaturation,
     ]
 }
diff --git a/pymc_marketing/mmm/transformers.py b/pymc_marketing/mmm/transformers.py
index 83f4e6d2..2ae28d86 100644
--- a/pymc_marketing/mmm/transformers.py
+++ b/pymc_marketing/mmm/transformers.py
@@ -899,13 +899,86 @@ def michaelis_menten(
     return alpha * x / (lam + x)
 
 
-def hill_saturation(
+def hill_function(
+    x: pt.TensorLike, slope: pt.TensorLike, kappa: pt.TensorLike
+) -> pt.TensorVariable:
+    r"""Hill Function
+
+    .. math::
+        f(x) = 1 - \frac{\kappa^s}{\kappa^s + x^s}
+
+    where:
+     - :math:`s` is the slope of the hill.
+     - :math:`\kappa` is the half-saturation point as :math:`f(\kappa) = 0.5` for any value of :math:`s` and :math:`\kappa`.
+     - :math:`x` is the independent variable and must be non-negative.
+
+    Hill function from Equation (5) in the paper [1]_.
+
+    .. plot::
+        :context: close-figs
+
+        import numpy as np
+        import matplotlib.pyplot as plt
+        from pymc_marketing.mmm.transformers import hill_function
+        x = np.linspace(0, 10, 100)
+        # Varying slope
+        slopes = [0.3, 0.7, 1.2]
+        fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
+        for i, slope in enumerate(slopes):
+            plt.subplot(1, 3, i+1)
+            y = hill_function(x, slope, 2).eval()
+            plt.plot(x, y)
+            plt.xlabel('x')
+            plt.title(f'Slope = {slope}')
+        plt.subplot(1,3,1)
+        plt.ylabel('Hill Saturation Sigmoid')
+        plt.tight_layout()
+        plt.show()
+        # Varying kappa
+        kappas = [1, 5, 10]
+        fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
+        for i, kappa in enumerate(kappas):
+            plt.subplot(1, 3, i+1)
+            y = hill_function(x, 1, kappa).eval()
+            plt.plot(x, y)
+            plt.xlabel('x')
+            plt.title(f'Kappa = {kappa}')
+        plt.subplot(1,3,1)
+        plt.ylabel('Hill Saturation Sigmoid')
+        plt.tight_layout()
+        plt.show()
+
+    Parameters
+    ----------
+    x : float or array-like
+        The independent variable, typically representing the concentration of a
+        substrate or the intensity of a stimulus.
+    slope : float
+        The slope of the hill. Must pe non-positive.
+    kappa : float
+        The half-saturation point as :math:`f(\kappa) = 0.5` for any value of :math:`s` and :math:`\kappa`.
+
+    Returns
+    -------
+    float
+        The value of the Hill function given the parameters.
+
+    References
+    ----------
+    .. [1] Jin, Yuxue, et al. “Bayesian methods for media mix modeling with carryover and shape effects.” (2017).
+    """  # noqa: E501
+    return pt.as_tensor_variable(
+        1 - pt.power(kappa, slope) / (pt.power(kappa, slope) + pt.power(x, slope))
+    )
+
+
+def hill_saturation_sigmoid(
     x: pt.TensorLike,
     sigma: pt.TensorLike,
     beta: pt.TensorLike,
     lam: pt.TensorLike,
 ) -> pt.TensorVariable:
-    r"""Hill Saturation Function
+    r"""Hill Saturation Sigmoid Function
 
     .. math::
         f(x) = \frac{\sigma}{1 + e^{-\beta(x - \lambda)}} - \frac{\sigma}{1 + e^{\beta\lambda}}
@@ -929,19 +1002,19 @@ def hill_saturation(
 
         import numpy as np
         import matplotlib.pyplot as plt
-        from pymc_marketing.mmm.transformers import hill_saturation
+        from pymc_marketing.mmm.transformers import hill_saturation_sigmoid
         x = np.linspace(0, 10, 100)
         # Varying sigma
         sigmas = [0.5, 1, 1.5]
         fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
         for i, sigma in enumerate(sigmas):
             plt.subplot(1, 3, i+1)
-            y = hill_saturation(x, sigma, 2, 5).eval()
+            y = hill_saturation_sigmoid(x, sigma, 2, 5).eval()
             plt.plot(x, y)
             plt.xlabel('x')
             plt.title(f'Sigma = {sigma}')
         plt.subplot(1,3,1)
-        plt.ylabel('Hill Saturation')
+        plt.ylabel('Hill Saturation Sigmoid')
         plt.tight_layout()
         plt.show()
         # Varying beta
@@ -949,12 +1022,12 @@ def hill_saturation(
         fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
         for i, beta in enumerate(betas):
             plt.subplot(1, 3, i+1)
-            y = hill_saturation(x, 1, beta, 5).eval()
+            y = hill_saturation_sigmoid(x, 1, beta, 5).eval()
             plt.plot(x, y)
             plt.xlabel('x')
             plt.title(f'Beta = {beta}')
         plt.subplot(1,3,1)
-        plt.ylabel('Hill Saturation')
+        plt.ylabel('Hill Saturation Sigmoid')
         plt.tight_layout()
         plt.show()
         # Varying lam
@@ -962,12 +1035,12 @@ def hill_saturation(
         fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
         for i, lam in enumerate(lams):
             plt.subplot(1, 3, i+1)
-            y = hill_saturation(x, 1, 2, lam).eval()
+            y = hill_saturation_sigmoid(x, 1, 2, lam).eval()
             plt.plot(x, y)
             plt.xlabel('x')
             plt.title(f'Lambda = {lam}')
         plt.subplot(1,3,1)
-        plt.ylabel('Hill Saturation')
+        plt.ylabel('Hill Saturation Sigmoid')
         plt.tight_layout()
         plt.show()
 
@@ -977,8 +1050,8 @@ def hill_saturation(
         The independent variable, typically representing the concentration of a
         substrate or the intensity of a stimulus.
     sigma : float
-        The upper asymptote of the curve, representing the maximum value the
-        function will approach as x grows large.
+        The upper asymptote of the curve, representing the approximate maximum value the
+        function will approach as x grows large. The true maximum value is at `sigma * (1 - 1 / (1 + exp(beta * lam)))`
     beta : float
         The slope parameter, determining the steepness of the curve.
     lam : float
@@ -988,7 +1061,7 @@ def hill_saturation(
     Returns
     -------
     float or array-like
-        The value of the Hill function for each input value of x.
+        The value of the Hill saturation sigmoid function for each input value of x.
     """
     return sigma / (1 + pt.exp(-beta * (x - lam))) - sigma / (1 + pt.exp(beta * lam))
 
diff --git a/tests/mmm/components/test_saturation.py b/tests/mmm/components/test_saturation.py
index 43245f9b..cebd9cd6 100644
--- a/tests/mmm/components/test_saturation.py
+++ b/tests/mmm/components/test_saturation.py
@@ -22,6 +22,7 @@
 
 from pymc_marketing.mmm import (
     HillSaturation,
+    HillSaturationSigmoid,
     InverseScaledLogisticSaturation,
     LogisticSaturation,
     MichaelisMentenSaturation,
@@ -48,6 +49,7 @@ def saturation_functions():
         TanhSaturationBaselined(),
         MichaelisMentenSaturation(),
         HillSaturation(),
+        HillSaturationSigmoid(),
         RootSaturation(),
     ]
 
@@ -104,6 +106,7 @@ def test_support_for_lift_test_integrations(saturation) -> None:
         ("tanh_baselined", TanhSaturationBaselined),
         ("michaelis_menten", MichaelisMentenSaturation),
         ("hill", HillSaturation),
+        ("hill_sigmoid", HillSaturationSigmoid),
         ("root", RootSaturation),
     ],
 )
diff --git a/tests/mmm/test_transformers.py b/tests/mmm/test_transformers.py
index 98179c39..c1454c32 100644
--- a/tests/mmm/test_transformers.py
+++ b/tests/mmm/test_transformers.py
@@ -27,7 +27,8 @@
     batched_convolution,
     delayed_adstock,
     geometric_adstock,
-    hill_saturation,
+    hill_function,
+    hill_saturation_sigmoid,
     inverse_scaled_logistic_saturation,
     logistic_saturation,
     michaelis_menten,
@@ -466,9 +467,9 @@ def test_michaelis_menten(self, x, alpha, lam, expected):
             (3, 2, -1),
         ],
     )
-    def test_hill_monotonicity(self, sigma, beta, lam):
+    def test_hill_sigmoid_monotonicity(self, sigma, beta, lam):
         x = np.linspace(-10, 10, 100)
-        y = hill_saturation(x, sigma, beta, lam).eval()
+        y = hill_saturation_sigmoid(x, sigma, beta, lam).eval()
         assert np.all(np.diff(y) >= 0), "The function is not monotonic."
 
     @pytest.mark.parametrize(
@@ -479,8 +480,8 @@ def test_hill_monotonicity(self, sigma, beta, lam):
             (3, 2, -1),
         ],
     )
-    def test_hill_zero(self, sigma, beta, lam):
-        y = hill_saturation(0, sigma, beta, lam).eval()
+    def test_hill_sigmoid_zero(self, sigma, beta, lam):
+        y = hill_saturation_sigmoid(0, sigma, beta, lam).eval()
         assert y == pytest.approx(0.0)
 
     @pytest.mark.parametrize(
@@ -491,8 +492,8 @@ def test_hill_zero(self, sigma, beta, lam):
             (-3, 3, 2, -1),
         ],
     )
-    def test_hill_sigma_upper_bound(self, x, sigma, beta, lam):
-        y = hill_saturation(x, sigma, beta, lam).eval()
+    def test_hill_sigmoid_sigma_upper_bound(self, x, sigma, beta, lam):
+        y = hill_saturation_sigmoid(x, sigma, beta, lam).eval()
         assert y <= sigma, f"The output {y} exceeds the upper bound sigma {sigma}."
 
     @pytest.mark.parametrize(
@@ -503,8 +504,8 @@ def test_hill_sigma_upper_bound(self, x, sigma, beta, lam):
             (-1, 3, 2, -1, 1.5),
         ],
     )
-    def test_hill_behavior_at_lambda(self, x, sigma, beta, lam, expected):
-        y = hill_saturation(x, sigma, beta, lam).eval()
+    def test_hill_sigmoid_behavior_at_lambda(self, x, sigma, beta, lam, expected):
+        y = hill_saturation_sigmoid(x, sigma, beta, lam).eval()
         offset = sigma / (1 + np.exp(beta * lam))
         expected_with_offset = expected - offset
         np.testing.assert_almost_equal(
@@ -522,8 +523,8 @@ def test_hill_behavior_at_lambda(self, x, sigma, beta, lam, expected):
             (np.array([1, 2, 3]), 3, 2, 2),
         ],
     )
-    def test_hill_vectorized_input(self, x, sigma, beta, lam):
-        y = hill_saturation(x, sigma, beta, lam).eval()
+    def test_hill_sigmoid_vectorized_input(self, x, sigma, beta, lam):
+        y = hill_saturation_sigmoid(x, sigma, beta, lam).eval()
         assert (
             y.shape == x.shape
         ), "The function did not return the correct shape for vectorized input."
@@ -536,9 +537,9 @@ def test_hill_vectorized_input(self, x, sigma, beta, lam):
             (3, 2, -1),
         ],
     )
-    def test_hill_asymptotic_behavior(self, sigma, beta, lam):
+    def test_hill_sigmoid_asymptotic_behavior(self, sigma, beta, lam):
         x = 1e6  # A very large value to approximate infinity
-        y = hill_saturation(x, sigma, beta, lam).eval()
+        y = hill_saturation_sigmoid(x, sigma, beta, lam).eval()
         offset = sigma / (1 + np.exp(beta * lam))
         expected = sigma - offset
         np.testing.assert_almost_equal(
@@ -548,6 +549,102 @@ def test_hill_asymptotic_behavior(self, sigma, beta, lam):
             err_msg="The function does not approach sigma as x approaches infinity.",
         )
 
+    @pytest.mark.parametrize(
+        argnames=["slope", "kappa"],
+        argvalues=[
+            (1, 1),
+            (2, 0.5),
+            (3, 2),
+        ],
+        ids=["slope=1, kappa=1", "slope=2, kappa=0.5", "slope=3, kappa=2"],
+    )
+    def test_hill_monotonicity(self, slope, kappa):
+        x = np.linspace(0, 10, 100)
+        y = hill_function(x, slope, kappa).eval()
+        assert np.all(np.diff(y) >= 0), "The function is not monotonic."
+
+    @pytest.mark.parametrize(
+        argnames=["slope", "kappa"],
+        argvalues=[
+            (1, 1),
+            (2, 0.5),
+            (3, 2),
+        ],
+        ids=["slope=1, kappa=1", "slope=2, kappa=0.5", "slope=3, kappa=2"],
+    )
+    def test_hill_zero(self, slope, kappa):
+        y = hill_function(0, slope, kappa).eval()
+        assert y == pytest.approx(0.0)
+
+    @pytest.mark.parametrize(
+        argnames=["x", "slope", "kappa"],
+        argvalues=[
+            (1, 1, 1),
+            (2, 0.5, 0.5),
+            (3, 2, 2),
+        ],
+        ids=[
+            "x=1, slope=1, kappa=1",
+            "x=2, slope=0.5, kappa=0.5",
+            "x=3, slope=2, kappa=2",
+        ],
+    )
+    def test_hill_upper_bound(self, x, slope, kappa):
+        y = hill_function(x, slope, kappa).eval()
+        assert y <= 1, f"The output {y} exceeds the upper bound 1."
+
+    @pytest.mark.parametrize(
+        argnames=["slope", "kappa"],
+        argvalues=[
+            (1, 1),
+            (2, 0.5),
+            (3, 2),
+        ],
+        ids=["slope=1, kappa=1", "slope=2, kappa=0.5", "slope=3, kappa=2"],
+    )
+    def test_hill_behavior_at_midpoint(self, slope, kappa):
+        y = hill_function(kappa, slope, kappa).eval()
+        expected = 0.5
+        np.testing.assert_almost_equal(
+            y,
+            expected,
+            decimal=5,
+            err_msg="The function does not behave as expected at the midpoint.",
+        )
+
+    @pytest.mark.parametrize(
+        "x, slope, kappa",
+        [
+            (np.array([0, 1, 2]), 1, 1),
+            (np.array([-1, 0, 1]), 2, 0.5),
+            (np.array([1, 2, 3]), 3, 2),
+        ],
+    )
+    def test_hill_vectorized_input(self, x, slope, kappa):
+        y = hill_function(x, slope, kappa).eval()
+        assert (
+            y.shape == x.shape
+        ), "The function did not return the correct shape for vectorized input."
+
+    @pytest.mark.parametrize(
+        "slope, kappa",
+        [
+            (1, 1),
+            (2, 0.5),
+            (3, 2),
+        ],
+    )
+    def test_hill_asymptotic_behavior(self, slope, kappa):
+        x = 1e6  # A very large value to approximate infinity
+        y = hill_function(x, slope, kappa).eval()
+        expected = 1
+        np.testing.assert_almost_equal(
+            y,
+            expected,
+            decimal=5,
+            err_msg="The function does not approach sigma as x approaches infinity.",
+        )
+
 
 class TestTransformersComposition:
     @pytest.mark.parametrize(