diff --git a/gammapy/datasets/core.py b/gammapy/datasets/core.py
index e58d3b2a69..78cd3662c6 100644
--- a/gammapy/datasets/core.py
+++ b/gammapy/datasets/core.py
@@ -10,6 +10,7 @@
 from gammapy.data import GTI
 from gammapy.modeling.models import DatasetModels, Models
 from gammapy.utils.scripts import make_name, make_path, read_yaml, to_yaml, write_yaml
+from gammapy.stats import prior_fit_statistic
 
 log = logging.getLogger(__name__)
 
@@ -68,15 +69,22 @@ def mask(self):
             return self.mask_safe
 
     def stat_sum(self):
-        """Total statistic given the current model parameters and priors."""
+        """Total statistic given the current model parameters."""
         stat = self.stat_array()
 
         if self.mask is not None:
             stat = stat[self.mask.data]
-        prior_stat_sum = 0.0
+
+        return np.sum(stat, dtype=np.float64)
+
+    def stat_sum_prior(self):
+        """Total statistic given the priors set on the models."""
         if self.models is not None:
-            prior_stat_sum = self.models.parameters.prior_stat_sum()
-        return np.sum(stat, dtype=np.float64) + prior_stat_sum
+            return prior_fit_statistic(self.models.priors)
+
+    def stat_sum_posterior(self):
+        """Total statistic given the current model parameters and priors set on the models."""
+        return self.stat_sum() + self.stat_sum_prior()
 
     @abc.abstractmethod
     def stat_array(self):
@@ -240,6 +248,18 @@ def stat_sum(self):
             stat_sum += dataset.stat_sum()
         return stat_sum
 
+    def stat_sum_prior(self):
+        """Compute joint statistic function value for priors."""
+        prior_stat_sum = 0
+        # Compute prior_fit_statistics from all datasets
+        for dataset in self:
+            prior_stat_sum += dataset.stat_sum_prior()
+        return prior_stat_sum
+
+    def stat_sum_posterior(self):
+        """Compute joint statistic function value for the posterior (parameter values and priors)."""
+        return self.stat_sum() + self.stat_sum_prior()
+
     def select_time(self, time_min, time_max, atol="1e-6 s"):
         """Select datasets in a given time interval.
 
diff --git a/gammapy/datasets/map.py b/gammapy/datasets/map.py
index adb6211231..8a55146214 100644
--- a/gammapy/datasets/map.py
+++ b/gammapy/datasets/map.py
@@ -17,6 +17,7 @@
     cash,
     cash_sum_cython,
     get_wstat_mu_bkg,
+    prior_fit_statistic,
     wstat,
 )
 from gammapy.utils.fits import HDULocation, LazyFitsData
@@ -407,7 +408,6 @@ def __init__(
     ):
         self._name = make_name(name)
         self._evaluators = {}
-
         self.counts = counts
         self.exposure = exposure
         self.background = background
@@ -564,7 +564,6 @@ def models(self, models):
                         use_cache=USE_NPRED_CACHE,
                     )
                     self._evaluators[model.name] = evaluator
-
         self._models = models
 
     @property
@@ -1335,20 +1334,22 @@ def plot_residuals(
         return ax_spatial, ax_spectral
 
     def stat_sum(self):
-        """Total statistic function value given the current model parameters and priors."""
-        prior_stat_sum = 0.0
-        if self.models is not None:
-            prior_stat_sum = self.models.parameters.prior_stat_sum()
-
+        """Total statistic function value given the current model parameters."""
         counts, npred = self.counts.data.astype(float), self.npred().data
 
         if self.mask is not None:
-            return (
-                cash_sum_cython(counts[self.mask.data], npred[self.mask.data])
-                + prior_stat_sum
-            )
+            return cash_sum_cython(counts[self.mask.data], npred[self.mask.data])
         else:
-            return cash_sum_cython(counts.ravel(), npred.ravel()) + prior_stat_sum
+            return cash_sum_cython(counts.ravel(), npred.ravel())
+
+    def stat_sum_prior(self):
+        """Total statistic function value given the priors set on the models."""
+        if self.models is not None:
+            return prior_fit_statistic(self.models.priors)
+
+    def stat_sum_posterior(self):
+        """Total statistic function value given the current model parameters and priors set on the models."""
+        return self.stat_sum() + self.stat_sum_prior()
 
     def fake(self, random_state="random-seed"):
         """Simulate fake counts for the current model and reduced IRFs.
diff --git a/gammapy/datasets/tests/test_datasets.py b/gammapy/datasets/tests/test_datasets.py
index d926349386..0d2d2af5dc 100644
--- a/gammapy/datasets/tests/test_datasets.py
+++ b/gammapy/datasets/tests/test_datasets.py
@@ -31,6 +31,9 @@ def test_datasets_types(datasets):
 
 
 def test_datasets_likelihood(datasets):
+
+    # rename the model of the second dataset for unique model names
+    datasets[1].models[0].name = datasets[1].name
     likelihood = datasets.stat_sum()
     assert_allclose(likelihood, 14472200.0002)
 
diff --git a/gammapy/datasets/tests/test_map.py b/gammapy/datasets/tests/test_map.py
index 40f07cec95..aaf8243bfb 100644
--- a/gammapy/datasets/tests/test_map.py
+++ b/gammapy/datasets/tests/test_map.py
@@ -747,16 +747,20 @@ def test_prior_stat_sum(sky_model, geom, geom_etrue):
     datasets.models = models
     dataset.counts = dataset.npred()
 
-    uniformprior = UniformPrior(min=-np.inf, max=0, weight=1)
-    datasets.models.parameters["amplitude"].prior = uniformprior
-    assert_allclose(datasets.stat_sum(), 12825.9370, rtol=1e-3)
+    UniformPrior(
+        min=-np.inf,
+        max=0,
+        weight=1,
+        modelparameters=datasets.models.parameters["amplitude"],
+    )
+    assert_allclose(datasets.stat_sum(), 12824.5062, rtol=1e-3)
 
     datasets.models.parameters["amplitude"].value = -1e-12
-    stat_sum_neg = datasets.stat_sum()
-    assert_allclose(stat_sum_neg, 470298.864993, rtol=1e-3)
+    stat_sum_post = datasets.stat_sum_posterior()
+    assert_allclose(stat_sum_post, 470332.1951, rtol=1e-3)
 
     datasets.models.parameters["amplitude"].prior.weight = 100
-    assert_allclose(datasets.stat_sum() - stat_sum_neg, 99, rtol=1e-3)
+    assert_allclose(datasets.stat_sum_posterior() - stat_sum_post, 99, rtol=1e-3)
 
 
 @requires_data()
diff --git a/gammapy/modeling/models/__init__.py b/gammapy/modeling/models/__init__.py
index 7aebf58f91..7b035df5bc 100644
--- a/gammapy/modeling/models/__init__.py
+++ b/gammapy/modeling/models/__init__.py
@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Built-in models in Gammapy."""
+
 from gammapy.utils.registry import Registry
 from .core import DatasetModels, Model, ModelBase, Models
 from .cube import (
@@ -8,7 +9,7 @@
     TemplateNPredModel,
     create_fermi_isotropic_diffuse_model,
 )
-from .prior import GaussianPrior, Prior, UniformPrior
+from .prior import MultiVariateGaussianPrior, GaussianPrior, Prior, UniformPrior
 from .spatial import (
     ConstantFluxSpatialModel,
     ConstantSpatialModel,
@@ -197,6 +198,7 @@
 
 PRIOR_REGISTRY = Registry(
     [
+        MultiVariateGaussianPrior,
         UniformPrior,
         GaussianPrior,
     ]
diff --git a/gammapy/modeling/models/core.py b/gammapy/modeling/models/core.py
index 821abf1be3..12ef350a40 100644
--- a/gammapy/modeling/models/core.py
+++ b/gammapy/modeling/models/core.py
@@ -37,7 +37,6 @@ def _get_model_class_from_dict(data):
     """Get a model class from a dictionary."""
     from . import (
         MODEL_REGISTRY,
-        PRIOR_REGISTRY,
         SPATIAL_MODEL_REGISTRY,
         SPECTRAL_MODEL_REGISTRY,
         TEMPORAL_MODEL_REGISTRY,
@@ -51,8 +50,6 @@ def _get_model_class_from_dict(data):
         cls = SPECTRAL_MODEL_REGISTRY.get_cls(data["spectral"]["type"])
     elif "temporal" in data:
         cls = TEMPORAL_MODEL_REGISTRY.get_cls(data["temporal"]["type"])
-    elif "prior" in data:
-        cls = PRIOR_REGISTRY.get_cls(data["prior"]["type"])
     return cls
 
 
@@ -419,6 +416,19 @@ def parameters(self):
         """Parameters as a `~gammapy.modeling.Parameters` object."""
         return Parameters.from_stack([_.parameters for _ in self._models])
 
+    @property
+    def priors(self):
+        """Priors (list).
+
+        Duplicate prior objects have been removed.
+        """
+        priors = {}
+
+        for parameter in self.parameters:
+            if parameter.prior is not None:
+                priors[parameter.prior] = parameter.prior
+        return list(priors)
+
     @property
     def parameters_unique_names(self):
         """List of unique parameter names. Return formatted as model_name.par_type.par_name."""
@@ -1219,10 +1229,6 @@ def insert(self, idx, model):
         _check_name_unique(model, self.names)
         self._models.insert(idx, model)
 
-    def set_prior(self, parameters, priors):
-        for parameter, prior in zip(parameters, priors):
-            parameter.prior = prior
-
 
 class restore_models_status:
     def __init__(self, models, restore_values=True):
diff --git a/gammapy/modeling/models/prior.py b/gammapy/modeling/models/prior.py
index 8897ae799a..10f1e95016 100644
--- a/gammapy/modeling/models/prior.py
+++ b/gammapy/modeling/models/prior.py
@@ -1,12 +1,13 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Priors for Gammapy."""
+
 import logging
 import numpy as np
 import astropy.units as u
-from gammapy.modeling import PriorParameter, PriorParameters
+from gammapy.modeling import Parameter, Parameters, PriorParameter, PriorParameters
 from .core import ModelBase
 
-__all__ = ["GaussianPrior", "UniformPrior"]
+__all__ = ["MultiVariateGaussianPrior", "GaussianPrior", "UniformPrior"]
 
 log = logging.getLogger(__name__)
 
@@ -37,7 +38,14 @@ class Prior(ModelBase):
 
     _unit = ""
 
-    def __init__(self, **kwargs):
+    def __init__(self, modelparameters, **kwargs):
+        if isinstance(modelparameters, Parameter):
+            self._modelparameters = Parameters([modelparameters])
+        elif isinstance(modelparameters, Parameters):
+            self._modelparameters = modelparameters
+        else:
+            raise ValueError(f"Invalid model type {modelparameters}")
+
         # Copy default parameters from the class to the instance
         default_parameters = self.default_parameters.copy()
 
@@ -57,6 +65,13 @@ def __init__(self, **kwargs):
         else:
             self._weight = 1
 
+        for par in self._modelparameters:
+            par.prior = self
+
+    @property
+    def modelparameters(self):
+        return self._modelparameters
+
     @property
     def parameters(self):
         """Prior parameters as a `~gammapy.modeling.PriorParameters` object."""
@@ -72,18 +87,18 @@ def __init_subclass__(cls, **kwargs):
 
     @property
     def weight(self):
-        """Weight mulitplied to the prior when evaluated."""
+        """Weight multiplied to the prior when evaluated."""
         return self._weight
 
     @weight.setter
     def weight(self, value):
         self._weight = value
 
-    def __call__(self, value):
+    def __call__(self):
         """Call evaluate method."""
-        # assuming the same unit as the PriorParamater here
+        # assuming the same unit as the PriorParameter here
         kwargs = {par.name: par.value for par in self.parameters}
-        return self.weight * self.evaluate(value.value, **kwargs)
+        return self.weight * self.evaluate(self._modelparameters.value, **kwargs)
 
     def to_dict(self, full_output=False):
         """Create dictionary for YAML serialisation."""
@@ -105,31 +120,83 @@ def to_dict(self, full_output=False):
                     ):
                         del par[item]
 
-        data = {"type": tag, "parameters": params, "weight": self.weight}
+        data = {
+            "type": tag,
+            "parameters": params,
+            "weight": self.weight,
+            "modelparameters": self._modelparameters,
+        }
 
-        if self.type is None:
-            return data
-        else:
-            return {self.type: data}
+        return data
 
     @classmethod
-    def from_dict(cls, data, **kwargs):
+    def from_dict(cls, data):
         """Get prior parameters from dictionary."""
+        from . import PRIOR_REGISTRY
+
+        prior_cls = PRIOR_REGISTRY.get_cls(data["type"])
         kwargs = {}
 
-        key0 = next(iter(data))
-        if key0 in ["prior"]:
-            data = data[key0]
-        if data["type"] not in cls.tag:
+        if data["type"] not in prior_cls.tag:
             raise ValueError(
                 f"Invalid model type {data['type']} for class {cls.__name__}"
             )
-
         priorparameters = _build_priorparameters_from_dict(
-            data["parameters"], cls.default_parameters
+            data["parameters"], prior_cls.default_parameters
         )
         kwargs["weight"] = data["weight"]
-        return cls.from_parameters(priorparameters, **kwargs)
+        kwargs["modelparameters"] = data["modelparameters"]
+
+        return prior_cls.from_parameters(priorparameters, **kwargs)
+
+
+class MultiVariateGaussianPrior(Prior):
+    """Multi-dimensional Gaussian Prior.
+
+    Parameters
+    ----------
+    modelparameters :
+        Meaning
+    covariance_matrix :
+        Meaning
+    """
+
+    tag = ["MultiVariateGaussianPrior"]
+    _type = "prior"
+
+    def __init__(self, modelparameters, covariance_matrix):
+        super().__init__(modelparameters)
+
+        self._covariance_matrix = covariance_matrix
+
+        value = np.asanyarray(self.covariance_matrix)
+        npars = len(self._modelparameters)
+        shape = (npars, npars)
+        if value.shape != shape:
+            raise ValueError(
+                f"Invalid covariance matrix shape: {value.shape}, expected {shape}"
+            )
+
+        # Do we want this?
+        self.dimension = value.shape[-1]
+
+        for par in self._modelparameters:
+            par.prior = self
+
+    # def from_subcovariance(self, parameters, subcovar):
+    #    idx = [parameters.index(par) for par in parameters]
+
+    def __call__(self):
+        """Call evaluate method"""
+        return self.evaluate(self._modelparameters.value)
+
+    @property
+    def covariance_matrix(self):
+        return self._covariance_matrix
+
+    def evaluate(self, values):
+        """Evaluate the MultiVariateGaussianPrior."""
+        return values.T @ self.covariance_matrix @ values
 
 
 class GaussianPrior(Prior):
diff --git a/gammapy/modeling/models/tests/test_prior.py b/gammapy/modeling/models/tests/test_prior.py
index 6f2fc4a978..f7ad3178f8 100644
--- a/gammapy/modeling/models/tests/test_prior.py
+++ b/gammapy/modeling/models/tests/test_prior.py
@@ -1,23 +1,26 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import pytest
 from numpy.testing import assert_allclose
-import astropy.units as u
-from gammapy.modeling.models import PRIOR_REGISTRY, GaussianPrior, Model, UniformPrior
+from gammapy.modeling import Parameter
+from gammapy.modeling.models import PRIOR_REGISTRY, GaussianPrior, Prior, UniformPrior
 from gammapy.utils.testing import assert_quantity_allclose
 
+TEST_PARAMETER = [dict(testparameter=Parameter(name="test", value=1))]
+
+
 TEST_PRIORS = [
     dict(
         name="gaussian",
-        model=GaussianPrior(mu=4.0, sigma=1.0),
-        val_at_0=16.0,
-        val_at_1=9.0,
-        val_with_weight_2=32.0,
+        model=GaussianPrior(mu=2.3, sigma=1.0, modelparameters=Parameter("index", 2.3)),
+        val_at_default=0.0,
+        val_at_0=5.29,
+        val_with_weight_2=10.58,
     ),
     dict(
         name="uni",
-        model=UniformPrior(min=0.0),
+        model=UniformPrior(min=0.0, modelparameters=Parameter("amplitude", 1e-12)),
+        val_at_default=1.0,
         val_at_0=0.0,
-        val_at_1=1.0,
         val_with_weight_2=0.0,
     ),
 ]
@@ -29,14 +32,15 @@ def test_priors(prior):
     for p in model.parameters:
         assert p.type == "prior"
 
-    value_0 = model(0.0 * u.Unit(""))
-    value_1 = model(1.0 * u.Unit(""))
+    value_default = model()
+    model.modelparameters[0].value = 0.0
+    value_0 = model()
+    assert_allclose(value_default, prior["val_at_default"], rtol=1e-7)
     assert_allclose(value_0, prior["val_at_0"], rtol=1e-7)
-    assert_allclose(value_1, prior["val_at_1"], rtol=1e-7)
 
     model.weight = 2.0
-    value_0_weight = model(0.0 * u.Unit(""))
-    assert_allclose(value_0_weight, prior["val_with_weight_2"], rtol=1e-7)
+    value_2_weight = model()
+    assert_allclose(value_2_weight, prior["val_with_weight_2"], rtol=1e-7)
 
 
 def test_to_from_dict():
@@ -45,9 +49,9 @@ def test_to_from_dict():
     model.weight = 2.0
     model_dict = model.to_dict()
     # Here we reverse the order of parameters list to ensure assignment is correct
-    model_dict["prior"]["parameters"].reverse()
+    model_dict["parameters"].reverse()
 
-    model_class = PRIOR_REGISTRY.get_cls(model_dict["prior"]["type"])
+    model_class = PRIOR_REGISTRY.get_cls(model_dict["type"])
     new_model = model_class.from_dict(model_dict)
 
     assert isinstance(new_model, UniformPrior)
@@ -57,7 +61,7 @@ def test_to_from_dict():
     assert_quantity_allclose(actual, desired)
     assert_allclose(model.weight, new_model.weight, rtol=1e-7)
 
-    new_model = Model.from_dict(model_dict)
+    new_model = Prior.from_dict(model_dict)
 
     assert isinstance(new_model, UniformPrior)
 
diff --git a/gammapy/modeling/parameter.py b/gammapy/modeling/parameter.py
index 6a23fe415b..5f6e27700f 100644
--- a/gammapy/modeling/parameter.py
+++ b/gammapy/modeling/parameter.py
@@ -153,7 +153,7 @@ def __init__(
         self.scan_n_sigma = scan_n_sigma
         self.interp = interp
         self.scale_method = scale_method
-        self.prior = prior
+        self._prior = prior
 
     def __get__(self, instance, owner):
         if instance is None:
@@ -184,20 +184,23 @@ def prior(self, value):
         if value is not None:
             from .models import Prior
 
+            # creating Prior out of dict
             if isinstance(value, dict):
-                from .models import Model
 
-                self._prior = Model.from_dict({"prior": value})
-            elif isinstance(value, Prior):
-                self._prior = value
+                value = Prior.from_dict(value)
+                if self in value.modelparameters:
+                    self._prior = value
+
+            # testing of self in prior.modelparameter
+            if isinstance(value, Prior):
+                if self in value.modelparameters:
+                    self._prior = value
+                else:
+                    raise TypeError(
+                        f"Invalid modelparameter: {value.modelparameters!r}"
+                    )
             else:
                 raise TypeError(f"Invalid type: {value!r}")
-        else:
-            self._prior = value
-
-    def prior_stat_sum(self):
-        if self.prior is not None:
-            return self.prior(self)
 
     @property
     def type(self):
@@ -465,7 +468,7 @@ def to_dict(self):
         if self._link_label_io is not None:
             output["link"] = self._link_label_io
         if self.prior is not None:
-            output["prior"] = self.prior.to_dict()["prior"]
+            output["prior"] = self.prior.to_dict()
         return output
 
     def autoscale(self):
diff --git a/gammapy/stats/__init__.py b/gammapy/stats/__init__.py
index b5216e5602..9f3f794f1d 100644
--- a/gammapy/stats/__init__.py
+++ b/gammapy/stats/__init__.py
@@ -2,7 +2,14 @@
 """Statistics."""
 
 from .counts_statistic import CashCountsStatistic, WStatCountsStatistic
-from .fit_statistics import cash, cstat, get_wstat_gof_terms, get_wstat_mu_bkg, wstat
+from .fit_statistics import (
+    cash,
+    cstat,
+    get_wstat_gof_terms,
+    get_wstat_mu_bkg,
+    prior_fit_statistic,
+    wstat,
+)
 from .fit_statistics_cython import (
     cash_sum_cython,
     f_cash_root_cython,
diff --git a/gammapy/stats/fit_statistics.py b/gammapy/stats/fit_statistics.py
index f67035242f..01df96655c 100644
--- a/gammapy/stats/fit_statistics.py
+++ b/gammapy/stats/fit_statistics.py
@@ -3,12 +3,41 @@
 
 see :ref:`fit-statistics`
 """
+
 import numpy as np
 from gammapy.stats.fit_statistics_cython import TRUNCATION_VALUE
 
 __all__ = ["cash", "cstat", "wstat", "get_wstat_mu_bkg", "get_wstat_gof_terms"]
 
 
+def prior_fit_statistic(priors):
+    r"""Prior.
+
+    Evaluating a list of priors.
+    Mulitdimensional priors are only counted once.
+
+    Parameters
+    ----------
+    priors : list
+        list of priors
+
+     Returns
+    -------
+    prior_stat_sum : float
+        sum of priors
+
+    """
+    prior_stat_sum = 0.0
+
+    if priors is not None:
+        _evaluated_priors = []
+        for prior in priors:
+            if prior not in _evaluated_priors:
+                prior_stat_sum += prior()
+                _evaluated_priors.append(prior)
+    return prior_stat_sum
+
+
 def cash(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
     r"""Cash statistic, for Poisson data.
 
@@ -170,9 +199,7 @@ def wstat(n_on, n_off, alpha, mu_sig, mu_bkg=None, extra_terms=True):
     with np.errstate(divide="ignore", invalid="ignore"):
         # This is a false positive error from pylint
         # See https://github.com/PyCQA/pylint/issues/2436
-        term2_ = -n_on * np.log(
-            mu_sig + alpha * mu_bkg
-        )  # pylint:disable=invalid-unary-operand-type
+        term2_ = -n_on * np.log(mu_sig + alpha * mu_bkg)  # pylint:disable=invalid-unary-operand-type
     # Handle n_on == 0
     condition = n_on == 0
     term2 = np.where(condition, 0, term2_)