diff --git a/docs/mmm/lifttest/index.ipy b/docs/mmm/lifttest/index.ipy
new file mode 100644
index 0000000..56eff42
--- /dev/null
+++ b/docs/mmm/lifttest/index.ipy
@@ -0,0 +1,278 @@
+# %% [markdown]
+# # Lift test
+
+import matplotlib.pyplot as plt
+import numpy as np
+# %%
+import pandas as pd
+
+index = pd.period_range("2000-01-01", "2005-01-01", freq="D")
+rng = np.random.default_rng(0)
+
+X = pd.DataFrame(
+    {
+        "investment1": np.cumsum(rng.normal(0, 1, size=len(index))),
+        "investment2": np.cumsum(rng.normal(0, 1, size=len(index))),
+    },
+    index=index,
+)
+X -= X.min()
+X /= X.max()
+X += 0.05
+
+X["investment2"] = X["investment1"] + 0.1 + rng.normal(0, 0.01, size=len(index))
+X.plot.line(alpha=0.9)
+
+import numpyro.distributions as dist
+
+# %%
+from prophetverse.effects import (HillEffect, LinearEffect,
+                                  LinearFourierSeasonality)
+from prophetverse.effects.trend import PiecewiseLinearTrend
+from prophetverse.engine import MAPInferenceEngine
+from prophetverse.engine.optimizer import LBFGSSolver
+from prophetverse.sktime import Prophetverse
+from prophetverse.utils.regex import exact, no_input_columns
+
+model = Prophetverse(
+    trend=PiecewiseLinearTrend(
+        changepoint_interval=100,
+        changepoint_prior_scale=0.001,
+        changepoint_range=-100,
+    ),
+    exogenous_effects=[
+        (
+            "seasonality",
+            LinearFourierSeasonality(
+                freq="D",
+                sp_list=[365.25],
+                fourier_terms_list=[3],
+                prior_scale=0.1,
+                effect_mode="multiplicative",
+            ),
+            no_input_columns,
+        ),
+        (
+            "investment1",
+            HillEffect(
+                half_max_prior=dist.HalfNormal(0.2),
+                slope_prior=dist.Gamma(2,1),
+                max_effect_prior=dist.HalfNormal(1.5),
+                effect_mode="additive",
+            ),
+            exact("investment1"),
+        ),
+        (
+            "investment2",
+            LinearEffect(
+                prior=dist.HalfNormal(0.5),
+                effect_mode="additive",
+            ),
+            exact("investment2"),
+        ),
+    ],
+    inference_engine=MAPInferenceEngine(num_steps=1),
+)
+
+
+from prophetverse.experimental.simulate import simulate
+
+samples = simulate(
+    model=model,
+    fh=X.index,
+    X=X,
+)
+
+# %%
+
+y = pd.DataFrame(data={"sales" : samples["obs"][0].flatten()}, index=index)
+true_effect = pd.DataFrame(
+    data={
+        "investment1": samples["investment1"][0].flatten(),
+        "investment2": samples["investment2"][0].flatten(),
+    },
+    index=index,
+    
+)
+
+fig, ax = plt.subplots(figsize=(10, 5))
+y.plot.line(ax=ax)
+fig.show()
+
+fig, ax = plt.subplots(figsize=(10, 5))
+true_effect.plot.line(ax=ax)
+fig.show()
+
+
+# %%
+
+import logging
+
+from prophetverse.engine.optimizer import LBFGSSolver
+from prophetverse.logger import logger
+
+logging.basicConfig()
+logging.root.setLevel(logging.DEBUG)
+
+model = model.set_params(
+    inference_engine=MAPInferenceEngine(
+                        num_steps=1000,
+                        optimizer=LBFGSSolver(memory_size=100,
+                                              max_linesearch_steps=100)
+))
+model.fit(y=y, X=X)
+
+components = model.predict_components(fh=index, X=X)
+
+
+# %%
+
+fig, ax = plt.subplots(figsize=(10, 5))
+components["obs"].plot.line(ax=ax)
+y.plot.line(ax=ax, color="black")
+# %%
+
+plt.figure()
+plt.scatter(X["investment1"], components["investment1"])
+plt.scatter(X["investment1"], true_effect["investment1"], color="black")
+
+
+plt.figure()
+plt.scatter(X["investment2"], components["investment2"])
+plt.scatter(X["investment2"], true_effect["investment2"], color="black")
+
+
+# %%
+def get_simulated_lift_test(X, n=10):
+    X_b = X.copy()
+
+    for col in ["investment1", "investment2"]:
+        
+        
+        X_b[col] = X_b[col]*rng.uniform(0.1, 0.9, size=X.shape[0])
+
+    samples_b = simulate(
+        model=model.clone().set_params(inference_engine__num_steps=1),
+        fh=X.index,
+        X=X_b,
+        do={k: v[0] for k, v in samples.items()},
+    )
+
+    true_effect_b = pd.DataFrame(
+        index=X_b.index,
+        data={
+            "investment1": samples_b["investment1"][0].flatten(),
+            "investment2": samples_b["investment2"][0].flatten(),
+        },
+    )
+
+    lift = np.abs(true_effect_b  - true_effect)
+
+    outs = []
+
+    for col in ["investment1", "investment2"]:
+        lift_test_dataframe = pd.DataFrame(
+            index=X.index,
+            data={
+                "lift": lift[col],
+                "x_start": X.loc[:, col],
+                "x_end": X_b.loc[:, col],
+                "y_start": true_effect.loc[:, col],
+                "y_end": true_effect_b.loc[:, col],
+            },
+        )
+        outs.append(lift_test_dataframe.sample(n=n, replace=False))
+
+    return tuple(outs)
+
+
+lift_test_dataframe1, lift_test_dataframe2 = get_simulated_lift_test(X,  n=7)
+# %%
+
+from prophetverse.effects.lift_likelihood import LiftExperimentLikelihood
+
+lift_experiment_effect1 = LiftExperimentLikelihood(
+       effect=model.get_params()["investment1"],
+        lift_test_results=lift_test_dataframe1,
+        prior_scale=1e-2,
+        likelihood_scale=1,
+)
+
+lift_experiment_effect2 = LiftExperimentLikelihood(
+    effect=model.get_params()["investment2"],
+    lift_test_results=lift_test_dataframe2,
+    prior_scale=1e-2,
+    likelihood_scale=1,
+)
+
+# %%
+
+from numpyro.infer.initialization import init_to_feasible, init_to_sample
+
+from prophetverse.engine.optimizer import LBFGSSolver
+
+new_model = model.clone()
+new_model.set_params(
+    investment1=lift_experiment_effect1,
+    investment2=lift_experiment_effect2,
+)
+new_model.fit(y=y, X=X)
+
+# %%
+new_components = new_model.predict_components(fh=index, X=X)
+
+# %%
+
+fig, ax = plt.subplots(figsize=(10, 5))
+components["obs"].plot.line(ax=ax)
+y.plot.line(ax=ax, color="black")
+new_components["obs"].plot.line(ax=ax)
+
+# %%
+
+plt.figure()
+plt.scatter(X["investment1"], components["investment1"])
+plt.scatter(X["investment1"], true_effect["investment1"], color="black")
+plt.scatter(X["investment1"], new_components["investment1"], color="tab:orange")
+
+for date in lift_test_dataframe1.index:
+    baseline = true_effect.loc[date, "investment1"]
+    plt.plot(
+        [
+            lift_test_dataframe1.loc[date, "x_start"],
+            lift_test_dataframe1.loc[date, "x_end"],
+        ],
+        [
+            lift_test_dataframe1.loc[date, "y_start"],
+            lift_test_dataframe1.loc[date, "y_end"],
+        ],
+        color="red",
+        alpha=0.8,
+    )
+plt.show()
+
+
+plt.figure()
+plt.scatter(X["investment2"], components["investment2"])
+plt.scatter(X["investment2"], true_effect["investment2"], color="black")
+plt.scatter(X["investment2"], new_components["investment2"], color="tab:orange")
+
+
+for date in lift_test_dataframe2.index:
+    baseline = true_effect.loc[date, "investment1"]
+    plt.plot(
+        [
+            lift_test_dataframe2.loc[date, "x_start"],
+            lift_test_dataframe2.loc[date, "x_end"],
+        ],
+        [ 
+            lift_test_dataframe2.loc[date, "y_start"],
+            lift_test_dataframe2.loc[date, "y_end"],
+        ],
+        color="red",
+        alpha=0.9,
+    )
+plt.show()
+
+
+# %%
diff --git a/extension_templates/effect.py b/extension_templates/effect.py
index f08f8e9..02319fe 100644
--- a/extension_templates/effect.py
+++ b/extension_templates/effect.py
@@ -9,6 +9,66 @@
 from prophetverse.utils.frame_to_array import series_to_tensor_or_array
 
 
+class MySimpleEffectName(BaseEffect):
+    """Base class for effects."""
+
+    _tags = {
+        # Supports multivariate data? Can this
+        # Effect be used with Multiariate prophet?
+        "supports_multivariate": False,
+        # If no columns are found, should
+        # _predict be skipped?
+        "skip_predict_if_no_match": True,
+        # Should only the indexes related to the forecasting horizon be passed to
+        # _transform?
+        "filter_indexes_with_forecating_horizon_at_transform": True,
+    }
+
+    def __init__(self, param1: Any, param2: Any):
+        self.param1 = param1
+        self.param2 = param2
+
+    def _sample_params(self, data, predicted_effects):
+        # call numpyro.sample to sample the parameters of the effect
+        # return a dictionary with the sampled parameters, where
+        # key is the name of the parameter and value is the sampled value
+        return {}
+
+    def _predict(
+        self,
+        data: Any,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
+    ) -> jnp.ndarray:
+        """Apply and return the effect values.
+
+        Parameters
+        ----------
+        data : Any
+            Data obtained from the transformed method.
+
+        predicted_effects : Dict[str, jnp.ndarray], optional
+            A dictionary containing the predicted effects, by default None.
+
+        params : Dict[str, jnp.ndarray]
+            A dictionary containing the sampled parameters of the effect.
+
+        Returns
+        -------
+        jnp.ndarray
+            An array with shape (T,1) for univariate timeseries, or (N, T, 1) for
+            multivariate timeseries, where T is the number of timepoints and N is the
+            number of series.
+        """
+        # predicted effects come with the following shapes:
+        # (T, 1) shaped array for univariate timeseries
+        # (N, T, 1) shaped array for multivariate timeseries, where N is the number of
+        # series
+
+        # Here you use the params to compute the effect.
+        raise NotImplementedError("Subclasses must implement _predict()")
+
+
 class MyEffectName(BaseEffect):
     """Base class for effects."""
 
@@ -76,10 +136,17 @@ def _transform(self, X: pd.DataFrame, fh: pd.Index) -> Any:
         array = series_to_tensor_or_array(X)
         return array
 
-    def predict(
+    def _sample_params(self, data, predicted_effects):
+        # call numpyro.sample to sample the parameters of the effect
+        # return a dictionary with the sampled parameters, where
+        # key is the name of the parameter and value is the sampled value
+        return {}
+
+    def _predict(
         self,
-        data: Dict,
+        data: Any,
         predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -91,6 +158,9 @@ def predict(
         predicted_effects : Dict[str, jnp.ndarray], optional
             A dictionary containing the predicted effects, by default None.
 
+        params : Dict[str, jnp.ndarray]
+            A dictionary containing the sampled parameters of the effect.
+
         Returns
         -------
         jnp.ndarray
@@ -98,11 +168,10 @@ def predict(
             multivariate timeseries, where T is the number of timepoints and N is the
             number of series.
         """
-        # Get the trend
+        # predicted effects come with the following shapes:
         # (T, 1) shaped array for univariate timeseries
         # (N, T, 1) shaped array for multivariate timeseries, where N is the number of
         # series
-        # trend: jnp.ndarray = predicted_effects["trend"]
-        # Or user predicted_effects.get("trend") to return None if the trend is
-        # not found
+
+        # Here you use the params to compute the effect.
         raise NotImplementedError("Subclasses must implement _predict()")
diff --git a/src/prophetverse/effects/__init__.py b/src/prophetverse/effects/__init__.py
index 1841178..901192a 100644
--- a/src/prophetverse/effects/__init__.py
+++ b/src/prophetverse/effects/__init__.py
@@ -1,9 +1,10 @@
 """Effects that define relationships between variables and the target."""
 
 from .base import BaseEffect
+from .exact_likelihood import ExactLikelihood
 from .fourier import LinearFourierSeasonality
 from .hill import HillEffect
-from .lift_experiment import LiftExperimentLikelihood
+from .lift_likelihood import LiftExperimentLikelihood
 from .linear import LinearEffect
 from .log import LogEffect
 
@@ -12,6 +13,7 @@
     "HillEffect",
     "LinearEffect",
     "LogEffect",
+    "ExactLikelihood",
     "LiftExperimentLikelihood",
     "LinearFourierSeasonality",
 ]
diff --git a/src/prophetverse/effects/base.py b/src/prophetverse/effects/base.py
index d2b1b2e..bfcc0c6 100644
--- a/src/prophetverse/effects/base.py
+++ b/src/prophetverse/effects/base.py
@@ -223,6 +223,7 @@ def predict(
         self,
         data: Dict,
         predicted_effects: Optional[Dict[str, jnp.ndarray]] = None,
+        params: Optional[Dict[str, jnp.ndarray]] = None,
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -244,14 +245,65 @@ def predict(
         if predicted_effects is None:
             predicted_effects = {}
 
-        x = self._predict(data, predicted_effects)
+        if params is None:
+            params = self.sample_params(data, predicted_effects)
+
+        x = self._predict(data, predicted_effects, params)
 
         return x
 
-    def _predict(
+    def sample_params(
         self,
         data: Dict,
+        predicted_effects: Optional[Dict[str, jnp.ndarray]] = None,
+    ):
+        """Sample parameters from the prior distribution.
+
+        Parameters
+        ----------
+        data : Dict
+            The data to be used for sampling the parameters, obtained from
+            `transform` method.
+
+        predicted_effects : Optional[Dict[str, jnp.ndarray]]
+            A dictionary containing the predicted effects, by default None.
+
+        Returns
+        -------
+        Dict
+            A dictionary containing the sampled parameters.
+        """
+        if predicted_effects is None:
+            predicted_effects = {}
+
+        return self._sample_params(data, predicted_effects)
+
+    def _sample_params(
+        self,
+        data: Any,
         predicted_effects: Dict[str, jnp.ndarray],
+    ):
+        """Sample parameters from the prior distribution.
+
+        Should be implemented by subclasses to provide the actual sampling logic.
+
+        Parameters
+        ----------
+        data : Any
+            The data to be used for sampling the parameters, obtained from
+            `transform` method.
+        predicted_effects : Dict[str, jnp.ndarray]
+            A dictionary containing the predicted effects, by default None.
+
+        Returns
+        -------
+        Dict
+            A dictionary containing the sampled parameters.
+        """
+        return {}
+
+    def _predict(
+        self, data: Dict, predicted_effects: Dict[str, jnp.ndarray], params: Dict
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -273,7 +325,10 @@ def _predict(
         raise NotImplementedError("Subclasses must implement _predict()")
 
     def __call__(
-        self, data: Dict, predicted_effects: Dict[str, jnp.ndarray]
+        self,
+        data: Dict,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: Optional[Dict[str, jnp.ndarray]] = None,
     ) -> jnp.ndarray:
         """Run the processes to calculate effect as a function."""
         return self.predict(data=data, predicted_effects=predicted_effects)
@@ -301,9 +356,10 @@ class BaseAdditiveOrMultiplicativeEffect(BaseEffect):
         or "multiplicative".
     """
 
-    def __init__(self, effect_mode="additive"):
+    def __init__(self, effect_mode="additive", base_effect_name: str = "trend"):
 
         self.effect_mode = effect_mode
+        self.base_effect_name = base_effect_name
 
         if effect_mode not in ["additive", "multiplicative"]:
             raise ValueError(
@@ -317,6 +373,7 @@ def predict(
         self,
         data: Any,
         predicted_effects: Optional[Dict[str, jnp.ndarray]] = None,
+        params: Optional[Dict[str, jnp.ndarray]] = None,
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -336,19 +393,29 @@ def predict(
             number of series.
         """
         if predicted_effects is None:
+            predicted_effects = {}
+
+        if params is None:
+            params = self.sample_params(data, predicted_effects)
+
+        if (
+            self.base_effect_name not in predicted_effects
+            and self.effect_mode == "multiplicative"
+        ):
             raise ValueError(
                 "BaseAdditiveOrMultiplicativeEffect requires trend in"
                 + " predicted_effects"
             )
 
-        trend = predicted_effects["trend"]
-        if trend.ndim == 1:
-            trend = trend.reshape((-1, 1))
-
-        x = super().predict(data=data, predicted_effects=predicted_effects)
-        x = x.reshape(trend.shape)
+        x = super().predict(
+            data=data, predicted_effects=predicted_effects, params=params
+        )
 
         if self.effect_mode == "additive":
             return x
 
-        return trend * x
+        base_effect = predicted_effects[self.base_effect_name]
+        if base_effect.ndim == 1:
+            base_effect = base_effect.reshape((-1, 1))
+        x = x.reshape(base_effect.shape)
+        return base_effect * x
diff --git a/src/prophetverse/effects/lift_experiment.py b/src/prophetverse/effects/exact_likelihood.py
similarity index 73%
rename from src/prophetverse/effects/lift_experiment.py
rename to src/prophetverse/effects/exact_likelihood.py
index c5416c2..ce2e3f5 100644
--- a/src/prophetverse/effects/lift_experiment.py
+++ b/src/prophetverse/effects/exact_likelihood.py
@@ -11,10 +11,10 @@
 
 from .base import BaseEffect
 
-__all__ = ["LiftExperimentLikelihood"]
+__all__ = ["ExactLikelihood"]
 
 
-class LiftExperimentLikelihood(BaseEffect):
+class ExactLikelihood(BaseEffect):
     """Wrap an effect and applies a normal likelihood to its output.
 
     This class uses an input as a reference for the effect, and applies a normal
@@ -22,10 +22,10 @@ class LiftExperimentLikelihood(BaseEffect):
 
     Parameters
     ----------
-    effect : BaseEffect
-        The effect to wrap.
-    lift_test_results : pd.DataFrame
-        A dataframe with the lift test results. Should be in sktime format, and must
+    effect_name : str
+        The effect to use in the likelihood.
+    reference_df : pd.DataFrame
+        A dataframe with the reference values. Should be in sktime format, and must
         have the same index as the input data.
     prior_scale : float
         The scale of the prior distribution for the likelihood.
@@ -35,20 +35,20 @@ class LiftExperimentLikelihood(BaseEffect):
 
     def __init__(
         self,
-        effect: BaseEffect,
-        lift_test_results: pd.DataFrame,
+        effect_name: str,
+        reference_df: pd.DataFrame,
         prior_scale: float,
     ):
 
-        self.effect = effect
-        self.lift_test_results = lift_test_results
+        self.effect_name = effect_name
+        self.reference_df = reference_df
         self.prior_scale = prior_scale
 
         assert self.prior_scale > 0, "prior_scale must be greater than 0"
 
         super().__init__()
 
-    def fit(self, y: pd.DataFrame, X: pd.DataFrame, scale: float = 1):
+    def _fit(self, y: pd.DataFrame, X: pd.DataFrame, scale: float = 1):
         """Initialize the effect.
 
         This method is called during `fit()` of the forecasting model.
@@ -74,9 +74,7 @@ def fit(self, y: pd.DataFrame, X: pd.DataFrame, scale: float = 1):
         -------
         None
         """
-        self.effect.fit(X=X, y=y, scale=scale)
         self.timeseries_scale = scale
-        super().fit(X=X, y=y, scale=scale)
 
     def _transform(self, X: pd.DataFrame, fh: pd.Index) -> Dict[str, Any]:
         """Prepare input data to be passed to numpyro model.
@@ -99,17 +97,19 @@ def _transform(self, X: pd.DataFrame, fh: pd.Index) -> Dict[str, Any]:
             Dictionary with data for the lift and for the inner effect
         """
         data_dict = {}
-        data_dict["inner_effect_data"] = self.effect._transform(X, fh=fh)
 
-        X_lift = self.lift_test_results.reindex(fh, fill_value=jnp.nan)
+        X_lift = self.reference_df.reindex(fh, fill_value=jnp.nan)
         lift_array = series_to_tensor_or_array(X_lift)
-        data_dict["observed_lift"] = lift_array / self.timeseries_scale
-        data_dict["obs_mask"] = ~jnp.isnan(data_dict["observed_lift"])
+        data_dict["observed_reference_value"] = lift_array / self.timeseries_scale
+        data_dict["obs_mask"] = ~jnp.isnan(data_dict["observed_reference_value"])
 
         return data_dict
 
     def _predict(
-        self, data: Dict, predicted_effects: Dict[str, jnp.ndarray]
+        self,
+        data: Dict,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -126,18 +126,16 @@ def _predict(
         jnp.ndarray
             An array with shape (T,1) for univariate timeseries.
         """
-        observed_lift = data["observed_lift"]
+        observed_reference_value = data["observed_reference_value"]
         obs_mask = data["obs_mask"]
 
-        x = self.effect.predict(
-            data=data["inner_effect_data"], predicted_effects=predicted_effects
-        )
+        x = predicted_effects[self.effect_name]
 
         with numpyro.handlers.mask(mask=obs_mask):
             numpyro.sample(
-                "lift_experiment",
+                "exact_likelihood:ignore",
                 dist.Normal(x, self.prior_scale),
-                obs=observed_lift,
+                obs=observed_reference_value,
             )
 
         return x
diff --git a/src/prophetverse/effects/fourier.py b/src/prophetverse/effects/fourier.py
index 71fd37b..3e167b0 100644
--- a/src/prophetverse/effects/fourier.py
+++ b/src/prophetverse/effects/fourier.py
@@ -124,10 +124,26 @@ def _transform(self, X: pd.DataFrame, fh: pd.Index) -> jnp.ndarray:
 
         return array
 
+    def _sample_params(self, data, predicted_effects=None):
+        """Sample parameters from the prior distribution.
+
+        Parameters
+        ----------
+        data : jnp.ndarray
+            The data to be used for sampling the parameters.
+
+        Returns
+        -------
+        dict
+            A dictionary containing the sampled parameters.
+        """
+        return self.linear_effect_.sample_params(data, predicted_effects)
+
     def _predict(
         self,
         data: Dict,
         predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -147,5 +163,7 @@ def _predict(
             number of series.
         """
         return self.linear_effect_.predict(
-            data=data, predicted_effects=predicted_effects
+            data=data,
+            predicted_effects=predicted_effects,
+            params=params,
         )
diff --git a/src/prophetverse/effects/hill.py b/src/prophetverse/effects/hill.py
index 227cc46..278df88 100644
--- a/src/prophetverse/effects/hill.py
+++ b/src/prophetverse/effects/hill.py
@@ -44,10 +44,35 @@ def __init__(
 
         super().__init__(effect_mode=effect_mode)
 
+    def _sample_params(
+        self, data, predicted_effects: Dict[str, jnp.ndarray]
+    ) -> Dict[str, jnp.ndarray]:
+        """
+        Sample the parameters of the effect.
+
+        Parameters
+        ----------
+        data : Any
+            Data obtained from the transformed method.
+        predicted_effects : Dict[str, jnp.ndarray]
+            A dictionary containing the predicted effects
+
+        Returns
+        -------
+        Dict[str, jnp.ndarray]
+            A dictionary containing the sampled parameters of the effect.
+        """
+        return {
+            "half_max": numpyro.sample("half_max", self.half_max_prior),
+            "slope": numpyro.sample("slope", self.slope_prior),
+            "max_effect": numpyro.sample("max_effect", self.max_effect_prior),
+        }
+
     def _predict(
         self,
         data: Dict[str, jnp.ndarray],
         predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -64,10 +89,11 @@ def _predict(
         jnp.ndarray
             An array with shape (T,1) for univariate timeseries.
         """
-        half_max = numpyro.sample("half_max", self.half_max_prior)
-        slope = numpyro.sample("slope", self.slope_prior)
-        max_effect = numpyro.sample("max_effect", self.max_effect_prior)
+        half_max = params["half_max"]
+        slope = params["slope"]
+        max_effect = params["max_effect"]
 
+        data = jnp.clip(data, 1e-9, None)
         x = _exponent_safe(data / half_max, -slope)
         effect = max_effect / (1 + x)
         return effect
diff --git a/src/prophetverse/effects/lift_likelihood.py b/src/prophetverse/effects/lift_likelihood.py
new file mode 100644
index 0000000..d86f5dd
--- /dev/null
+++ b/src/prophetverse/effects/lift_likelihood.py
@@ -0,0 +1,211 @@
+"""Composition of effects (Effects that wrap other effects)."""
+
+from typing import Any, Dict
+
+import jax.numpy as jnp
+import numpyro
+import pandas as pd
+
+from prophetverse.distributions import GammaReparametrized
+from prophetverse.utils.frame_to_array import series_to_tensor_or_array
+
+from .base import BaseEffect
+
+__all__ = ["LiftExperimentLikelihood"]
+
+
+class LiftExperimentLikelihood(BaseEffect):
+    """Wrap an effect and applies a normal likelihood to its output.
+
+    This class uses an input as a reference for the effect, and applies a normal
+    likelihood to the output of the effect.
+
+    Parameters
+    ----------
+    effect : BaseEffect
+        The effect to wrap.
+    lift_test_results : pd.DataFrame
+        A dataframe with the lift test results. Should be in sktime format, and must
+        have the same index as the input data.
+    prior_scale : float
+        The scale of the prior distribution for the likelihood.
+    """
+
+    _tags = {"skip_predict_if_no_match": False, "supports_multivariate": False}
+
+    def __init__(
+        self,
+        effect: BaseEffect,
+        lift_test_results: pd.DataFrame,
+        prior_scale: float,
+        likelihood_scale: float = 1,
+    ):
+
+        self.effect = effect
+        self.lift_test_results = lift_test_results
+        self.prior_scale = prior_scale
+        self.likelihood_scale = likelihood_scale
+
+        super().__init__()
+
+        assert self.prior_scale > 0, "prior_scale must be greater than 0"
+
+        mandatory_columns = ["x_start", "x_end", "lift"]
+        assert all(
+            column in self.lift_test_results.columns for column in mandatory_columns
+        ), f"lift_test_results must have the following columns: {mandatory_columns}"
+
+    def fit(self, y: pd.DataFrame, X: pd.DataFrame, scale: float = 1):
+        """Initialize the effect.
+
+        This method is called during `fit()` of the forecasting model.
+        It receives the Exogenous variables DataFrame and should be used to initialize
+        any necessary parameters or data structures, such as detecting the columns that
+        match the regex pattern.
+
+        This method MUST set _input_feature_columns_names to a list of column names
+
+        Parameters
+        ----------
+        y : pd.DataFrame
+            The timeseries dataframe
+
+        X : pd.DataFrame
+            The DataFrame to initialize the effect.
+
+        scale : float, optional
+            The scale of the timeseries. For multivariate timeseries, this is
+            a dataframe. For univariate, it is a simple float.
+
+        Returns
+        -------
+        None
+        """
+        self.effect_ = self.effect.clone()
+        self.effect_.fit(X=X, y=y, scale=scale)
+        self.timeseries_scale = scale
+        super().fit(X=X, y=y, scale=scale)
+
+    def _transform(self, X: pd.DataFrame, fh: pd.Index) -> Dict[str, Any]:
+        """Prepare input data to be passed to numpyro model.
+
+        Returns a dictionary with the data for the lift and for the inner effect.
+
+        Parameters
+        ----------
+        X : pd.DataFrame
+            The input DataFrame containing the exogenous variables for the training
+            time indexes, if passed during fit, or for the forecasting time indexes, if
+            passed during predict.
+
+        fh : pd.Index
+            The forecasting horizon as a pandas Index.
+
+        Returns
+        -------
+        Dict[str, Any]
+            Dictionary with data for the lift and for the inner effect
+        """
+        data_dict = {}
+        data_dict["inner_effect_data"] = self.effect_._transform(X, fh=fh)
+
+        # Check if fh and self.lift_test_results have same index type
+        if not isinstance(fh, self.lift_test_results.index.__class__):
+            raise TypeError(
+                "fh and self.lift_test_results must have the same index type"
+            )
+        X_lift = self.lift_test_results.reindex(fh, fill_value=jnp.nan)
+
+        data_dict["observed_lift"] = (
+            series_to_tensor_or_array(X_lift["lift"].dropna()) / self.timeseries_scale
+        )
+        data_dict["x_start"] = series_to_tensor_or_array(X_lift["x_start"].dropna())
+        data_dict["x_end"] = series_to_tensor_or_array(X_lift["x_end"].dropna())
+        data_dict["obs_mask"] = ~jnp.isnan(series_to_tensor_or_array(X_lift["lift"]))
+
+        return data_dict
+
+    def _sample_params(self, data, predicted_effects):
+        """
+        Sample the parameters of the effect.
+
+        Calls the sample_params method of the inner effect.
+
+        Parameters
+        ----------
+        data : Any
+            Data obtained from the transformed method.
+        predicted_effects : Dict[str, jnp.ndarray]
+            A dictionary containing the predicted effects
+
+        Returns
+        -------
+        Dict[str, jnp.ndarray]
+            A dictionary containing the sampled parameters of the effect.
+        """
+        return self.effect_.sample_params(
+            data=data["inner_effect_data"], predicted_effects=predicted_effects
+        )
+
+    def _predict(
+        self,
+        data: Dict,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
+    ) -> jnp.ndarray:
+        """Apply and return the effect values.
+
+        Parameters
+        ----------
+        data : Any
+            Data obtained from the transformed method.
+
+        predicted_effects : Dict[str, jnp.ndarray], optional
+            A dictionary containing the predicted effects, by default None.
+
+        Returns
+        -------
+        jnp.ndarray
+            An array with shape (T,1) for univariate timeseries.
+        """
+        observed_lift = data["observed_lift"].reshape((-1, 1))
+        x_start = data["x_start"].reshape((-1, 1))
+        x_end = data["x_end"].reshape((-1, 1))
+        obs_mask = data["obs_mask"]
+
+        predicted_effects_masked = {
+            k: v[obs_mask] for k, v in predicted_effects.items()
+        }
+
+        # Call the effect a first time
+        x = self.effect_.predict(
+            data=data["inner_effect_data"],
+            predicted_effects=predicted_effects,
+            params=params,
+        )
+
+        # Get the start and end values
+        y_start = self.effect_.predict(
+            data=x_start,
+            predicted_effects=predicted_effects_masked,
+            params=params,
+        )
+        y_end = self.effect_.predict(
+            data=x_end, predicted_effects=predicted_effects_masked, params=params
+        )
+
+        # Calculate the delta_y
+        delta_y = jnp.abs(y_end - y_start)
+
+        with numpyro.handlers.scale(scale=self.likelihood_scale):
+            distribution = GammaReparametrized(delta_y, self.prior_scale)
+
+            # Add :ignore so that the model removes this
+            # sample when organizing the output dataframe
+            numpyro.sample(
+                "lift_experiment:ignore",
+                distribution,
+                obs=observed_lift,
+            )
+
+        return x
diff --git a/src/prophetverse/effects/linear.py b/src/prophetverse/effects/linear.py
index f7077f1..f6c26fb 100644
--- a/src/prophetverse/effects/linear.py
+++ b/src/prophetverse/effects/linear.py
@@ -40,10 +40,22 @@ def __init__(
 
         super().__init__(effect_mode=effect_mode)
 
+    def _sample_params(self, data, predicted_effects):
+
+        n_features = data.shape[-1]
+
+        with numpyro.plate("features_plate", n_features, dim=-1):
+            coefficients = numpyro.sample("coefs", self.prior)
+
+        return {
+            "coefficients": coefficients,
+        }
+
     def _predict(
         self,
         data: Any,
-        predicted_effects: Optional[Dict[str, jnp.ndarray]] = None,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -52,8 +64,8 @@ def _predict(
         data : Any
             Data obtained from the transformed method.
 
-        predicted_effects : Dict[str, jnp.ndarray], optional
-            A dictionary containing the predicted effects, by default None.
+        predicted_effects : Dict[str, jnp.ndarray]
+            A dictionary containing the predicted effects
 
         Returns
         -------
@@ -62,10 +74,7 @@ def _predict(
             multivariate timeseries, where T is the number of timepoints and N is the
             number of series.
         """
-        n_features = data.shape[-1]
-
-        with numpyro.plate("features_plate", n_features, dim=-1):
-            coefficients = numpyro.sample("coefs", self.prior)
+        coefficients = params["coefficients"]
 
         if coefficients.ndim == 1:
             coefficients = jnp.expand_dims(coefficients, axis=-1)
diff --git a/src/prophetverse/effects/log.py b/src/prophetverse/effects/log.py
index dee6fba..c6235e7 100644
--- a/src/prophetverse/effects/log.py
+++ b/src/prophetverse/effects/log.py
@@ -38,10 +38,19 @@ def __init__(
         self.rate_prior = rate_prior or dist.Gamma(1, 1)
         super().__init__(effect_mode=effect_mode)
 
+    def _sample_params(self, data, predicted_effects):
+        scale = numpyro.sample("log_scale", self.scale_prior)
+        rate = numpyro.sample("log_rate", self.rate_prior)
+        return {
+            "scale": scale,
+            "rate": rate,
+        }
+
     def _predict(  # type: ignore[override]
         self,
         data: jnp.ndarray,
-        predicted_effects: Optional[Dict[str, jnp.ndarray]] = None,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
     ) -> jnp.ndarray:
         """Apply and return the effect values.
 
@@ -60,8 +69,9 @@ def _predict(  # type: ignore[override]
             multivariate timeseries, where T is the number of timepoints and N is the
             number of series.
         """
-        scale = numpyro.sample("log_scale", self.scale_prior)
-        rate = numpyro.sample("log_rate", self.rate_prior)
+        scale = params["scale"]
+        rate = params["rate"]
+
         effect = scale * jnp.log(jnp.clip(rate * data + 1, 1e-8, None))
 
         return effect
diff --git a/src/prophetverse/effects/trend/flat.py b/src/prophetverse/effects/trend/flat.py
index ce07c16..d44fec1 100644
--- a/src/prophetverse/effects/trend/flat.py
+++ b/src/prophetverse/effects/trend/flat.py
@@ -1,5 +1,7 @@
 """Flat trend model."""
 
+from typing import Any, Dict
+
 import jax.numpy as jnp
 import numpyro
 import numpyro.distributions as dist
@@ -56,8 +58,34 @@ def _transform(self, X: pd.DataFrame, fh: pd.Index) -> dict:
         idx = X.index
         return jnp.ones((len(idx), 1))
 
+    def _sample_params(self, data: Any, predicted_effects: Dict[str, jnp.ndarray]):
+        """Sample parameters from the prior distribution.
+
+        Parameters
+        ----------
+        data : jnp.ndarray
+            A constant vector with the size of the series time indexes
+
+        Returns
+        -------
+        dict
+            dictionary containing the sampled parameters for the trend model
+        """
+        return {
+            "trend_flat_coefficient": numpyro.sample(
+                "trend_flat_coefficient",
+                dist.Gamma(
+                    rate=self.changepoint_prior_loc / self.changepoint_prior_scale**2,
+                    concentration=self.changepoint_prior_loc,
+                ),
+            ),
+        }
+
     def _predict(  # type: ignore[override]
-        self, data: jnp.ndarray, predicted_effects=None
+        self,
+        data: jnp.ndarray,
+        predicted_effects: dict,
+        params: dict,
     ) -> jnp.ndarray:
         """Apply the trend.
 
@@ -74,18 +102,6 @@ def _predict(  # type: ignore[override]
         # Alias for clarity
         constant_vector = data
 
-        mean = self.changepoint_prior_loc
-        var = self.changepoint_prior_scale**2
-
-        rate = mean / var
-        concentration = mean * rate
-
-        coefficient = numpyro.sample(
-            "trend_flat_coefficient",
-            dist.Gamma(
-                rate=rate,
-                concentration=concentration,
-            ),
-        )
+        coefficient = params["trend_flat_coefficient"]
 
         return constant_vector * coefficient
diff --git a/src/prophetverse/effects/trend/piecewise.py b/src/prophetverse/effects/trend/piecewise.py
index 308782d..42b9010 100644
--- a/src/prophetverse/effects/trend/piecewise.py
+++ b/src/prophetverse/effects/trend/piecewise.py
@@ -5,7 +5,7 @@
 """
 
 import itertools
-from typing import Dict, Tuple, Union
+from typing import Any, Dict, Tuple, Union
 
 import jax.numpy as jnp
 import numpy as np
@@ -157,8 +157,33 @@ def _transform(self, X: pd.DataFrame, fh: pd.Index) -> dict:
         idx = self._fh_to_index(fh)
         return self.get_changepoint_matrix(idx)
 
+    def _sample_params(self, data: Any, predicted_effects: Dict[str, jnp.ndarray]):
+
+        changepoint_matrix = data
+
+        offset = numpyro.sample(
+            "offset",
+            dist.Normal(self._offset_prior_loc, self._offset_prior_scale),
+        )
+        changepoint_coefficients = numpyro.sample(
+            "changepoint_coefficients",
+            dist.Laplace(self._changepoint_prior_loc, self._changepoint_prior_scale),
+        )
+
+        if changepoint_matrix.ndim == 3:
+            changepoint_coefficients = changepoint_coefficients.reshape((1, -1, 1))
+            offset = offset.reshape((-1, 1, 1))
+
+        return {
+            "changepoint_coefficients": changepoint_coefficients,
+            "offset": offset,
+        }
+
     def _predict(
-        self, data: jnp.ndarray, predicted_effects: Dict[str, jnp.ndarray]
+        self,
+        data: jnp.ndarray,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: dict,
     ) -> jnp.ndarray:
         """
         Compute the trend based on the given changepoint matrix.
@@ -178,20 +203,8 @@ def _predict(
         """
         # alias for clarity
         changepoint_matrix = data
-        offset = numpyro.sample(
-            "offset",
-            dist.Normal(self._offset_prior_loc, self._offset_prior_scale),
-        )
-
-        changepoint_coefficients = numpyro.sample(
-            "changepoint_coefficients",
-            dist.Laplace(self._changepoint_prior_loc, self._changepoint_prior_scale),
-        )
-
-        # If multivariate
-        if changepoint_matrix.ndim == 3:
-            changepoint_coefficients = changepoint_coefficients.reshape((1, -1, 1))
-            offset = offset.reshape((-1, 1, 1))
+        changepoint_coefficients = params["changepoint_coefficients"]
+        offset = params["offset"]
 
         trend = (changepoint_matrix) @ changepoint_coefficients + offset
 
@@ -533,8 +546,38 @@ def _suggest_global_trend_and_offset(
 
         return global_rates, offset
 
+    def _sample_params(self, data, predicted_effects):
+        """
+        Sample params for the effect.
+
+        Use super to sample the changepoint coefficients and offset, and then sample
+        the capacity using the capacity prior.
+
+        Parameters
+        ----------
+        data : Any
+            The input data.
+        predicted_effects : Dict[str, jnp.ndarray]
+            The predicted effects
+
+        Returns
+        -------
+        dict
+            The sampled parameters.
+        """
+        with numpyro.plate("series", self.n_series, dim=-3):
+            capacity = numpyro.sample("capacity", self.capacity_prior)
+
+        return {
+            "capacity": capacity,
+            **super()._sample_params(data=data, predicted_effects=predicted_effects),
+        }
+
     def _predict(  # type: ignore[override]
-        self, data: jnp.ndarray, predicted_effects=None
+        self,
+        data: Any,
+        predicted_effects: Dict[str, jnp.ndarray],
+        params: Dict[str, jnp.ndarray],
     ) -> jnp.ndarray:
         """
         Compute the trend for the given changepoint matrix.
@@ -549,10 +592,11 @@ def _predict(  # type: ignore[override]
         jnp.ndarray
             The computed trend.
         """
-        with numpyro.plate("series", self.n_series, dim=-3):
-            capacity = numpyro.sample("capacity", self.capacity_prior)
+        trend = super()._predict(
+            data=data, predicted_effects=predicted_effects, params=params
+        )
 
-        trend = super()._predict(data=data, predicted_effects=predicted_effects)
+        capacity = params["capacity"]
         if self.n_series == 1:
             capacity = capacity.squeeze()
 
diff --git a/src/prophetverse/sktime/base.py b/src/prophetverse/sktime/base.py
index 0c409e7..065172c 100644
--- a/src/prophetverse/sktime/base.py
+++ b/src/prophetverse/sktime/base.py
@@ -430,6 +430,13 @@ def _get_predictive_samples_dict(
         predict_data = self._get_predict_data(X=X, fh=fh)
 
         predictive_samples_ = self.inference_engine_.predict(**predict_data)
+
+        keys_to_delete = []
+        for key in predictive_samples_.keys():
+            if key.endswith(":ignore"):
+                keys_to_delete.append(key)
+        for key in keys_to_delete:
+            del predictive_samples_[key]
         return predictive_samples_
 
     def predict_all_sites_samples(self, fh, X=None):
diff --git a/tests/conftest.py b/tests/conftest.py
index ec879c8..288482b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,14 +1,8 @@
 """Configure tests and declare global fixtures."""
 
-import warnings  # noqa: F401
-
-import jax.numpy as jnp
-import numpyro  # noqa: F401
 import pandas as pd
 import pytest
 
-from prophetverse.effects.base import BaseAdditiveOrMultiplicativeEffect
-
 # warnings.filterwarnings("ignore")
 
 
@@ -28,25 +22,3 @@ def sample_data():
             "lin_x2": [0.2 * i for i in range(10, 20)],
         }
     )
-
-
-class ConcreteEffect(BaseAdditiveOrMultiplicativeEffect):
-    """Most simple class to test abstracteffect methods."""
-
-    _tags = {"skip_predict_if_no_match": False}
-
-    def _predict(self, data: jnp.ndarray, predicted_effects=None) -> jnp.ndarray:
-        """Calculate simple effect."""
-        return jnp.mean(data, axis=1)
-
-
-@pytest.fixture(name="effect_with_regex")
-def effect_with_regex():
-    """Most simple class of abstracteffect with optional regex."""
-    return ConcreteEffect()
-
-
-@pytest.fixture
-def effect_without_regex():
-    """Most simple class of abstracteffect without optional regex."""
-    return ConcreteEffect()
diff --git a/tests/effects/test_base.py b/tests/effects/test_base.py
index 97a71ff..29e2dd4 100644
--- a/tests/effects/test_base.py
+++ b/tests/effects/test_base.py
@@ -5,6 +5,28 @@
 from prophetverse.effects.base import BaseAdditiveOrMultiplicativeEffect, BaseEffect
 
 
+class ConcreteEffect(BaseAdditiveOrMultiplicativeEffect):
+    """Most simple class to test abstracteffect methods."""
+
+    _tags = {"skip_predict_if_no_match": False}
+
+    def _predict(self, data, predicted_effects, params) -> jnp.ndarray:
+        """Calculate simple effect."""
+        return jnp.mean(data, axis=1, keepdims=True)
+
+
+@pytest.fixture(name="effect_with_regex")
+def effect_with_regex():
+    """Most simple class of abstracteffect with optional regex."""
+    return ConcreteEffect()
+
+
+@pytest.fixture
+def effect_without_regex():
+    """Most simple class of abstracteffect without optional regex."""
+    return ConcreteEffect()
+
+
 @pytest.mark.smoke
 def test__predict(effect_with_regex):
     trend = jnp.array([1.0, 2.0, 3.0]).reshape((-1, 1))
@@ -13,14 +35,9 @@ def test__predict(effect_with_regex):
     expected_result = jnp.mean(data, axis=1).reshape((-1, 1))
     assert jnp.allclose(result, expected_result)
 
+    call_result = effect_with_regex(data=data, predicted_effects={"trend": trend})
 
-@pytest.mark.smoke
-def test_call(effect_with_regex):
-    trend = jnp.array([1.0, 2.0, 3.0]).reshape((-1, 1))
-    data = jnp.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]).reshape((-1, 2))
-    result = effect_with_regex(data=data, predicted_effects={"trend": trend})
-    expected_result = jnp.mean(data, axis=1).reshape((-1, 1))
-    assert jnp.allclose(result, expected_result)
+    assert jnp.all(call_result == result)
 
 
 def test_bad_effect_mode():
diff --git a/tests/effects/test_exact_likelihood.py b/tests/effects/test_exact_likelihood.py
new file mode 100644
index 0000000..23e6827
--- /dev/null
+++ b/tests/effects/test_exact_likelihood.py
@@ -0,0 +1,88 @@
+import jax.numpy as jnp
+import numpyro.distributions as dist
+import pandas as pd
+import pytest
+from numpyro import handlers
+
+from prophetverse.effects import ExactLikelihood, LinearEffect
+
+
+@pytest.fixture
+def exact_likelihood_results():
+    return pd.DataFrame(
+        data={"test_results": [1, 2, 3, 4, 5, 6]},
+        index=pd.date_range("2021-01-01", periods=6),
+    )
+
+
+@pytest.fixture
+def inner_effect():
+    return LinearEffect(prior=dist.Delta(2))
+
+
+@pytest.fixture
+def exact_likelihood_effect_instance(exact_likelihood_results):
+    return ExactLikelihood(
+        effect_name="exog",
+        reference_df=exact_likelihood_results,
+        prior_scale=1.0,
+    )
+
+
+@pytest.fixture
+def X():
+    return pd.DataFrame(
+        data={"exog": [10, 20, 30, 40, 50, 60]},
+        index=pd.date_range("2021-01-01", periods=6),
+    )
+
+
+@pytest.fixture
+def y(X):
+    return pd.DataFrame(index=X.index, data=[1] * len(X))
+
+
+def test_exact_likelihood_initialization(
+    exact_likelihood_effect_instance, exact_likelihood_results
+):
+    assert exact_likelihood_effect_instance.reference_df.equals(
+        exact_likelihood_results
+    )
+    assert exact_likelihood_effect_instance.prior_scale == 1.0
+
+
+def test_exact_likelihood_fit(X, exact_likelihood_effect_instance):
+
+    exact_likelihood_effect_instance.fit(y=y, X=X, scale=1)
+    assert exact_likelihood_effect_instance.timeseries_scale == 1
+    assert exact_likelihood_effect_instance._is_fitted
+
+
+def test_exact_likelihood_transform_train(X, y, exact_likelihood_effect_instance):
+    fh = y.index.get_level_values(-1).unique()
+    exact_likelihood_effect_instance.fit(X=X, y=y)
+    transformed = exact_likelihood_effect_instance.transform(
+        X,
+        fh=fh,
+    )
+    assert "observed_reference_value" in transformed
+    assert transformed["observed_reference_value"] is not None
+
+
+def test_exact_likelihood_predict(X, y, exact_likelihood_effect_instance):
+    fh = X.index.get_level_values(-1).unique()
+
+    exog = jnp.array([1, 2, 3, 4, 5, 6]).reshape((-1, 1))
+    exact_likelihood_effect_instance.fit(X=X, y=y)
+    data = exact_likelihood_effect_instance.transform(X=X, fh=fh)
+
+    exec_trace = handlers.trace(exact_likelihood_effect_instance.predict).get_trace(
+        data=data, predicted_effects={"exog": exog}
+    )
+
+    assert len(exec_trace) == 1
+
+    trace_likelihood = exec_trace["exact_likelihood:ignore"]
+    assert trace_likelihood["type"] == "sample"
+    assert jnp.all(trace_likelihood["value"] == exog)
+    assert trace_likelihood["is_observed"]
diff --git a/tests/effects/test_lift_experiment.py b/tests/effects/test_lift_experiment.py
deleted file mode 100644
index 9f755cd..0000000
--- a/tests/effects/test_lift_experiment.py
+++ /dev/null
@@ -1,83 +0,0 @@
-import jax.numpy as jnp
-import numpyro.distributions as dist
-import pandas as pd
-import pytest
-
-from prophetverse.effects import LiftExperimentLikelihood, LinearEffect
-
-
-@pytest.fixture
-def lift_test_results():
-    return pd.DataFrame(
-        data={"test_results": [1, 2, 3, 4, 5, 6]},
-        index=pd.date_range("2021-01-01", periods=6),
-    )
-
-
-@pytest.fixture
-def inner_effect():
-    return LinearEffect(prior=dist.Delta(2))
-
-
-@pytest.fixture
-def lift_experiment_effect_instance(lift_test_results, inner_effect):
-    return LiftExperimentLikelihood(
-        effect=inner_effect, lift_test_results=lift_test_results, prior_scale=1.0
-    )
-
-
-@pytest.fixture
-def X():
-    return pd.DataFrame(
-        data={"exog": [10, 20, 30, 40, 50, 60]},
-        index=pd.date_range("2021-01-01", periods=6),
-    )
-
-
-@pytest.fixture
-def y(X):
-    return pd.DataFrame(index=X.index, data=[1] * len(X))
-
-
-def test_liftexperimentlikelihood_initialization(
-    lift_experiment_effect_instance, inner_effect, lift_test_results
-):
-    assert lift_experiment_effect_instance.effect == inner_effect
-    assert lift_experiment_effect_instance.lift_test_results.equals(lift_test_results)
-    assert lift_experiment_effect_instance.prior_scale == 1.0
-
-
-def test_liftexperimentlikelihood_fit(X, lift_experiment_effect_instance):
-
-    lift_experiment_effect_instance.fit(y=y, X=X, scale=1)
-    assert lift_experiment_effect_instance.timeseries_scale == 1
-    assert lift_experiment_effect_instance.effect._is_fitted
-
-
-def test_liftexperimentlikelihood_transform_train(
-    X, y, lift_experiment_effect_instance
-):
-    fh = y.index.get_level_values(-1).unique()
-    lift_experiment_effect_instance.fit(X=X, y=y)
-    transformed = lift_experiment_effect_instance.transform(
-        X,
-        fh=fh,
-    )
-    assert "observed_lift" in transformed
-    assert transformed["observed_lift"] is not None
-
-
-def test_liftexperimentlikelihood_predict(X, y, lift_experiment_effect_instance):
-    fh = X.index.get_level_values(-1).unique()
-
-    trend = jnp.array([1, 2, 3, 4, 5, 6])
-    lift_experiment_effect_instance.fit(X=X, y=y)
-    data = lift_experiment_effect_instance.transform(X=X, fh=fh)
-    predicted = lift_experiment_effect_instance.predict(
-        data=data, predicted_effects={"trend": trend}
-    )
-    inner_effect_data = lift_experiment_effect_instance.effect.transform(X, fh=fh)
-    inner_effect_predict = lift_experiment_effect_instance.effect.predict(
-        data=inner_effect_data, predicted_effects={"trend": trend}
-    )
-    assert jnp.all(predicted == inner_effect_predict)
diff --git a/tests/effects/test_lift_test_likelihood.py b/tests/effects/test_lift_test_likelihood.py
new file mode 100644
index 0000000..2f833c5
--- /dev/null
+++ b/tests/effects/test_lift_test_likelihood.py
@@ -0,0 +1,90 @@
+import jax.numpy as jnp
+import numpyro.distributions as dist
+import pandas as pd
+import pytest
+from numpyro import handlers
+
+from prophetverse.effects import LiftExperimentLikelihood, LinearEffect
+
+
+@pytest.fixture
+def lift_test_results():
+    index = pd.date_range("2021-01-01", periods=2)
+    return pd.DataFrame(
+        index=index, data={"x_start": [1, 2], "x_end": [10, 20], "lift": [2, 4]}
+    )
+
+
+@pytest.fixture
+def inner_effect():
+    return LinearEffect(prior=dist.Delta(2), effect_mode="additive")
+
+
+@pytest.fixture
+def lift_experiment_likelihood_effect_instance(lift_test_results, inner_effect):
+    return LiftExperimentLikelihood(
+        effect=inner_effect,
+        lift_test_results=lift_test_results,
+        prior_scale=1.0,
+    )
+
+
+@pytest.fixture
+def X():
+    return pd.DataFrame(
+        data={"exog": [10, 20, 30, 40, 50, 60]},
+        index=pd.date_range("2021-01-01", periods=6),
+    )
+
+
+@pytest.fixture
+def y(X):
+    return pd.DataFrame(index=X.index, data=[1] * len(X))
+
+
+def test_lift_experiment_likelihood_initialization(
+    lift_experiment_likelihood_effect_instance, lift_test_results
+):
+    assert lift_experiment_likelihood_effect_instance.lift_test_results.equals(
+        lift_test_results
+    )
+    assert lift_experiment_likelihood_effect_instance.prior_scale == 1.0
+
+
+def test_lift_experiment_likelihood_fit(X, lift_experiment_likelihood_effect_instance):
+
+    lift_experiment_likelihood_effect_instance.fit(y=y, X=X, scale=1)
+    assert lift_experiment_likelihood_effect_instance.timeseries_scale == 1
+    assert lift_experiment_likelihood_effect_instance.effect_._is_fitted
+
+
+def test_lift_experiment_likelihood_transform_train(
+    X, y, lift_experiment_likelihood_effect_instance, lift_test_results
+):
+    fh = y.index.get_level_values(-1).unique()
+    lift_experiment_likelihood_effect_instance.fit(X=X, y=y)
+    transformed = lift_experiment_likelihood_effect_instance.transform(
+        X,
+        fh=fh,
+    )
+    assert "observed_lift" in transformed
+    assert len(transformed["observed_lift"]) == len(lift_test_results)
+
+
+def test_lift_experiment_likelihood_predict(
+    X, y, lift_experiment_likelihood_effect_instance
+):
+    fh = X.index.get_level_values(-1).unique()
+
+    exog = jnp.array([1, 2, 3, 4, 5, 6]).reshape((-1, 1))
+    lift_experiment_likelihood_effect_instance.fit(X=X, y=y)
+    data = lift_experiment_likelihood_effect_instance.transform(X=X, fh=fh)
+
+    exec_trace = handlers.trace(
+        lift_experiment_likelihood_effect_instance.predict
+    ).get_trace(data=data, predicted_effects={"exog": exog})
+
+    assert "lift_experiment:ignore" in exec_trace
+    trace_likelihood = exec_trace["lift_experiment:ignore"]
+    assert trace_likelihood["type"] == "sample"
+    assert trace_likelihood["is_observed"]