diff --git a/cuqi/implicitprior/_restorator.py b/cuqi/implicitprior/_restorator.py
index 16569017a3..a7649e001c 100644
--- a/cuqi/implicitprior/_restorator.py
+++ b/cuqi/implicitprior/_restorator.py
@@ -98,7 +98,7 @@ def get_conditioning_variables(self):
 
 
 class MoreauYoshidaPrior(Distribution):
-    """    
+    r"""    
     This class defines (implicit) smoothed priors for which we can apply 
     gradient-based algorithms. The smoothing is performed using
     the Moreau-Yoshida envelope of the target prior potential.
@@ -106,42 +106,47 @@ class MoreauYoshidaPrior(Distribution):
     In the following we give a detailed explanation of the
     Moreau-Yoshida smoothing.
     
-    We consider a density such that - \log\pi(x) = -g(x) with g convex, lsc,
+    We consider a density such that :math:`- \log\pi(x) = -g(x)` with g convex, lsc,
     proper but not differentiable. Consequently, we cannot apply any
     algorithm requiring the gradient of g.
     Idea:
     We consider the Moreau envelope of g defined as
     
-    g_{smoothing_strength} (x) = inf_z 0.5*\| x-z \|_2^2/smoothing_strength + g(z).
+    .. math::
+
+        g_{\mathrm{smoothing\_strength}} (x) = \inf_z 0.5 \| x-z \|_2^2/\mathrm{smoothing\_strength} + g(z).
     
-    g_{smoothing_strength} has some nice properties
-        - g_{smoothing_strength}(x)-->g(x) as smoothing_strength-->0 for all x
-        - \nabla g_{smoothing_strength} is 1/smoothing_strength-Lipschitz
-        - \nabla g_{smoothing_strength}(x) = (x - prox_g^{smoothing_strength}(x))/smoothing_strength for all x with 
+    :math:`g_{\mathrm{smoothing\_strength}}` has some nice properties
+
+    - :math:`g_{\mathrm{smoothing\_strength}}(x) \to g(x)` as :math:`\mathrm{smoothing\_strength} \to 0` for all x
+    - :math:`\nabla g_{\mathrm{smoothing\_strength}}` is :math:`1/\mathrm{smoothing\_strength}`-Lipschitz
+    - :math:`\nabla g_{\mathrm{smoothing\_strength}}(x) = (x - \mathrm{prox}_g^{\mathrm{smoothing\_strength}}(x))/\mathrm{smoothing\_strength}` for all x with 
         
-        prox_g^{smoothing_strength}(x) = argmin_z 0.5*\| x-z \|_2^2/smoothing_strength + g(z) .
+        .. math::
+
+            \mathrm{prox}_g^{\mathrm{smoothing\_strength}}(x) = \operatorname{argmin}_z 0.5 \| x-z \|_2^2/\mathrm{smoothing\_strength} + g(z) .
 
     Consequently, we can apply any gradient-based algorithm with
-    g_{smoothing_strength} in lieu of g. These algorithms do not require the
-    full knowledge of g_{smoothing_strength} but only its gradient. The gradient
-    of g_{smoothing_strength} is fully determined by prox_g^{smoothing_strength}
-    and smoothing_strength.
+    :math:`g_{\mathrm{smoothing\_strength}}` in lieu of g. These algorithms do not require the
+    full knowledge of :math:`g_{\mathrm{smoothing\_strength}}` but only its gradient. The gradient
+    of :math:`g_{\mathrm{smoothing\_strength}}` is fully determined by :math:`\mathrm{prox}_g^{\mathrm{smoothing\_strength}}`
+    and :math:`\mathrm{smoothing\_strength}`.
     It is important as, although there exists an explicit formula for
-    g_{smoothing_strength}, it is rarely used in practice, as it would require
+    :math:`g_{\mathrm{smoothing\_strength}}`, it is rarely used in practice, as it would require
     us to solve an optimization problem each time we want to 
-    estimate g_{smoothing_strength}. Furthermore, there exist cases where we dont't
-    the regularization g with which the mapping prox_g^{smoothing_strength} is
+    estimate :math:`g_{\mathrm{smoothing\_strength}}`. Furthermore, there exist cases where we dont't
+    the regularization g with which the mapping :math:`\mathrm{prox}_g^{\mathrm{smoothing\_strength}}` is
     associated.
 
     Remark (Proximal operators are denoisers):
         We consider the denoising inverse problem x = u + n, with
-        n \sim \mathcal{N}(0, smoothing_strength I).
+        :math:`n \sim \mathcal{N}(0, \mathrm{smoothing\_strength} \, I)`.
         A mapping solving a denoising inverse problem is called denoiser. It takes
         the noisy observation x as an input and returns a less noisy version of x
         which is an estimate of u.
-        We assume a prior density \pi(u) \propto exp(- g(u)).
+        We assume a prior density :math:`\pi(u) \propto \exp(- g(u))`.
         Then the MAP estimate is given by 
-            x_MAP = \argmin_z 0.5 \| x - z \|_2^2/smoothing_strength + g(z) = prox_g^smoothing_strength(x)
+            :math:`x_{\mathrm{MAP}} = \operatorname{argmin}_z 0.5 \| x - z \|_2^2/\mathrm{smoothing\_strength} + g(z) = \mathrm{prox}_g^{\mathrm{smoothing\_strength}}(x)`
         Then proximal operators are denoisers. 
     
     Remark (Denoisers are not necessarily proximal operators): Data-driven
diff --git a/demos/howtos/defining_posterior.py b/demos/howtos/defining_posterior.py
index 39154634d5..c1b3636ecf 100644
--- a/demos/howtos/defining_posterior.py
+++ b/demos/howtos/defining_posterior.py
@@ -29,10 +29,8 @@
 # Then consider the following Bayesian model
 #
 # .. math::
-#    \begin{align*}
-#    \mathbf{x} &\sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
-#    \mathbf{y} &\sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\,\mathbf{I})
-#    \end{align*}
+#    \mathbf{x} \sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
+#    \mathbf{y} \sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\,\mathbf{I})
 #
 # which can be written in CUQIpy as
 
@@ -64,10 +62,8 @@
 # Suppose we had two forward models :math:`\mathbf{A}` and :math:`\mathbf{B}`:
 #
 # .. math::
-#    \begin{align*}
-#    \mathbf{y} &= \mathbf{A}\mathbf{x}\\
-#    \mathbf{d} &= \mathbf{B}\mathbf{x}\\
-#    \end{align*}
+#    \mathbf{y} = \mathbf{A}\mathbf{x}\\
+#    \mathbf{d} = \mathbf{B}\mathbf{x}
 
 # Both observations come from the same unknown x
 A, y_obs, _ = cuqi.testproblem.Deconvolution1D().get_components()
@@ -77,11 +73,9 @@
 # Then consider the following Bayesian model
 #
 # .. math::
-#    \begin{align*}
-#    \mathbf{x} &\sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
-#    \mathbf{y} &\sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\mathbf{I})\\
-#    \mathbf{d} &\sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01^2\mathbf{I})
-#    \end{align*}
+#    \mathbf{x} \sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
+#    \mathbf{y} \sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\mathbf{I})\\
+#    \mathbf{d} \sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01^2\mathbf{I})
 
 x = cuqi.distribution.Gaussian(np.zeros(A.domain_dim), 0.1)
 y = cuqi.distribution.Gaussian(A(x), 0.05**2)
@@ -115,11 +109,9 @@
 # we have the following 3 forward models
 #
 # .. math::
-#    \begin{align*}
-#    \mathbf{y} &= \mathbf{A}\mathbf{x}\\
-#    \mathbf{d} &= \mathbf{B}\mathbf{x}\\
-#    \mathbf{b} &= C(\mathbf{x})
-#    \end{align*}
+#    \mathbf{y} = \mathbf{A}\mathbf{x}\\
+#    \mathbf{d} = \mathbf{B}\mathbf{x}\\
+#    \mathbf{b} = C(\mathbf{x})
 #
 # where :math:`C` is a nonlinear function.
 
@@ -133,15 +125,13 @@
 # Then consider the following Bayesian model
 #
 # .. math::
-#    \begin{align*}
-#    q          &\sim \mathcal{U}(0.1, 10)\\
-#    l          &\sim \mathrm{Gamma}(1, 1)\\
-#    s          &\sim \mathrm{Gamma}(1, 10^{-2})\\
-#    \mathbf{x} &\sim \mathcal{N}(\mathbf{0}, l^{-1}\mathbf{I})\\
-#    \mathbf{y} &\sim \mathcal{N}(\mathbf{A}\mathbf{x}, s^{-1}\mathbf{I})\\
-#    \mathbf{d} &\sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01\mathbf{I})\\
-#    \mathbf{b} &\sim \mathcal{L}(\mathbf{C}(\mathbf{x}), q)
-#    \end{align*}
+#    q          \sim \mathcal{U}(0.1, 10)\\
+#    l          \sim \mathrm{Gamma}(1, 1)\\
+#    s          \sim \mathrm{Gamma}(1, 10^{-2})\\
+#    \mathbf{x} \sim \mathcal{N}(\mathbf{0}, l^{-1}\mathbf{I})\\
+#    \mathbf{y} \sim \mathcal{N}(\mathbf{A}\mathbf{x}, s^{-1}\mathbf{I})\\
+#    \mathbf{d} \sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01\mathbf{I})\\
+#    \mathbf{b} \sim \mathcal{L}(\mathbf{C}(\mathbf{x}), q)
 
 q = cuqi.distribution.Uniform(0.1, 10)
 l = cuqi.distribution.Gamma(1, 1)
diff --git a/demos/howtos/gibbs.py b/demos/howtos/gibbs.py
index 21c65dd469..a6ac36cf41 100644
--- a/demos/howtos/gibbs.py
+++ b/demos/howtos/gibbs.py
@@ -2,28 +2,28 @@
 Gibbs sampling
 ==============
 
-    This tutorial shows how to use CUQIpy to perform Gibbs sampling.
-    Gibbs sampling is a Markov chain Monte Carlo (MCMC) method for
-    sampling a joint probability distribution.
-
-    Opposed to jointly sampling the distribution simultaneously, Gibbs
-    sampling samples the variables of the distribution sequentially,
-    one variable at a time. When a variable represents a random vector, the
-    whole vector is sampled simultaneously.
-    
-    The sampling of each variable is done by sampling from the conditional
-    distribution of that variable given (fixed, previously sampled) values
-    of the other variables.
-
-    This is often a very efficient way of sampling from a joint
-    distribution if the conditional distributions are easy to sample
-    from. This is one way to exploit the structure of the joint
-    distribution. On the other hand, if the conditional distributions
-    are highly correlated and/or are difficult to sample from, then
-    Gibbs sampling can be very inefficient.
-
-    For these reasons, Gibbs sampling is often a double-edged sword,
-    that needs to be used in the right context.
+This tutorial shows how to use CUQIpy to perform Gibbs sampling.
+Gibbs sampling is a Markov chain Monte Carlo (MCMC) method for
+sampling a joint probability distribution.
+
+Opposed to jointly sampling the distribution simultaneously, Gibbs
+sampling samples the variables of the distribution sequentially,
+one variable at a time. When a variable represents a random vector, the
+whole vector is sampled simultaneously.
+
+The sampling of each variable is done by sampling from the conditional
+distribution of that variable given (fixed, previously sampled) values
+of the other variables.
+
+This is often a very efficient way of sampling from a joint
+distribution if the conditional distributions are easy to sample
+from. This is one way to exploit the structure of the joint
+distribution. On the other hand, if the conditional distributions
+are highly correlated and/or are difficult to sample from, then
+Gibbs sampling can be very inefficient.
+
+For these reasons, Gibbs sampling is often a double-edged sword,
+that needs to be used in the right context.
 
 """
 # %%
@@ -77,12 +77,10 @@
 #
 # .. math::
 #
-#     \begin{align}
-#         d &\sim \mathrm{Gamma}(1, 10^{-4}) \\
-#         l &\sim \mathrm{Gamma}(1, 10^{-4}) \\
-#         \mathbf{x} &\sim \mathrm{GMRF}(\mathbf{0}, d) \\
-#         \mathbf{y} &\sim \mathcal{N}(\mathbf{A} \mathbf{x}, l^{-1} \mathbf{I}_m)
-#     \end{align}
+#     d \sim \mathrm{Gamma}(1, 10^{-4}) \\
+#     l \sim \mathrm{Gamma}(1, 10^{-4}) \\
+#     \mathbf{x} \sim \mathrm{GMRF}(\mathbf{0}, d) \\
+#     \mathbf{y} \sim \mathcal{N}(\mathbf{A} \mathbf{x}, l^{-1} \mathbf{I}_m)
 #
 # where :math:`\mathbf{y}` is the observed data, and :math:`\mathbf{x}`
 # is the unknown signal. The hyperparameters :math:`d` and :math:`l` are
diff --git a/docs/api/index.rst b/docs/api/index.rst
index edc414633a..d4aed56474 100644
--- a/docs/api/index.rst
+++ b/docs/api/index.rst
@@ -7,7 +7,7 @@ The main modules are:
    - :doc:`cuqi.distribution <_autosummary/cuqi.distribution>` for defining probability distributions.
    - :doc:`cuqi.model <_autosummary/cuqi.model>` for defining deterministic models.
    - :doc:`cuqi.sampler <_autosummary/cuqi.sampler>` for sampling from probability distributions.
-   - :doc:`cuqi.inference <_autosummary/cuqi.implicitprior>` for defining implicit priors.
+   - :doc:`cuqi.implicitprior <_autosummary/cuqi.implicitprior>` for defining implicit priors.
 
 The following modules provide higher-level interfaces:
    - :doc:`cuqi.testproblem <_autosummary/cuqi.testproblem>` for defining specific test problems.