diff --git a/species/data/database.py b/species/data/database.py
index 7fc5ccd..c091ef6 100644
--- a/species/data/database.py
+++ b/species/data/database.py
@@ -2163,7 +2163,7 @@ def add_samples(
         samples : np.ndarray
             Samples of the posterior.
         ln_prob : np.ndarray
-            Log posterior for each sample.
+            Log-likelihood for each sample.
         modelpar : list(str)
             List with the model parameter names.
         bounds : dict
@@ -2275,7 +2275,7 @@ def get_probable_sample(
     ) -> Dict[str, float]:
         """
         Function for extracting the sample parameters
-        with the highest posterior probability.
+        with the maximum likelihood.
 
         Parameters
         ----------
@@ -2293,11 +2293,11 @@ def get_probable_sample(
         -------
         dict
             Parameters and values for the sample with the
-            maximum posterior probability.
+            maximum likelihood.
         """
 
         if verbose:
-            print_section("Get sample with highest probability")
+            print_section("Get sample with the maximum likelihood")
             print(f"Database tag: {tag}")
 
         if burnin is None:
@@ -2327,10 +2327,8 @@ def get_probable_sample(
                 samples = np.reshape(samples, (-1, n_param))
                 ln_prob = np.reshape(ln_prob, -1)
 
-            index_max = np.unravel_index(ln_prob.argmax(), ln_prob.shape)
-
-            # max_prob = ln_prob[index_max]
-            max_sample = samples[index_max]
+            index_max = np.argmax(ln_prob)
+            max_sample = samples[index_max, ]
 
             prob_sample = {}
 
@@ -3310,7 +3308,7 @@ def get_samples(
                                 f"({norm_prior[0]}, {norm_prior[1]})"
                             )
 
-                            normal_priors[f"{bound_item}/{filter_item}"] = (
+                            normal_priors[f"{prior_item}/{filter_item}"] = (
                                 norm_prior[0],
                                 norm_prior[1],
                             )
@@ -3320,7 +3318,7 @@ def get_samples(
 
                         print(f"   - {prior_item} = ({norm_prior[0]}, {norm_prior[1]})")
 
-                        normal_priors[bound_item] = (norm_prior[0], norm_prior[1])
+                        normal_priors[prior_item] = (norm_prior[0], norm_prior[1])
 
         median_sample = self.get_median_sample(tag, burnin, verbose=False)
         prob_sample = self.get_probable_sample(tag, burnin, verbose=False)
diff --git a/species/plot/plot_mcmc.py b/species/plot/plot_mcmc.py
index 2f25847..8735433 100644
--- a/species/plot/plot_mcmc.py
+++ b/species/plot/plot_mcmc.py
@@ -247,8 +247,10 @@ def plot_posterior(
         Axis limits of all parameters. Automatically set if the
         argument is set to ``None``.
     max_prob : bool
-        Plot the position of the sample with the
-        maximum posterior probability.
+        Plot the position of the sample with the maximum likelihood
+        probability. The sample may need to match with the maximum
+        posterior probability, depending on the use of any normal
+        priors.
     vmr : bool
         Plot the volume mixing ratios (i.e. number fractions)
         instead of the mass fractions of the retrieved species with
@@ -998,7 +1000,7 @@ def plot_posterior(
 
     if max_prob:
         max_idx = np.argmax(box.ln_prob)
-        max_sample = samples[max_idx, :]
+        max_sample = samples[max_idx, ]
 
     if isinstance(title_fmt, list) and len(title_fmt) != ndim:
         raise ValueError(