def deviance_information_criterions(mean_posterior_lls, ll_per_sample):
    r"""Calculates the Deviance Information Criteria (DIC) using three methods.

    This returns a dictionary returning the ``DIC_2002``, the ``DIC_2004`` and the ``DIC_Ando_2011`` method.
    The first is based on Spiegelhalter et al (2002), the second based on Gelman et al. (2004) and the last on
    Ando (2011). All cases differ in how they calculate model complexity, i.e. the effective number of parameters
    in the model. In all cases the model with the smallest DIC is preferred.

    All these DIC methods measure fitness using the deviance, which is, for a likelihood :math:`p(y | \theta)`
    defined as:

    .. math::

        D(\theta) = -2\log p(y|\theta)

    From this, the posterior mean deviance,

    .. math::

        \bar{D} = \mathbb{E}_{\theta}[D(\theta)]

    is then used as a measure of how well the model fits the data.

    The complexity, or measure of effective number of parameters, can be measured in see ways, see
    Spiegelhalter et al. (2002), Gelman et al (2004) and Ando (2011). The first method calculated the parameter
    deviance as:

    .. math::
        :nowrap:

        \begin{align}
        p_{D} &= \mathbb{E}_{\theta}[D(\theta)] - D(\mathbb{E}[\theta)]) \\
              &= \bar{D} - D(\bar{\theta})
        \end{align}

    i.e. posterior mean deviance minus the deviance evaluated at the posterior mean of the parameters.

    The second method calculated :math:`p_{D}` as:

    .. math::

        p_{D} = p_{V} = \frac{1}{2}\hat{var}(D(\theta))

    i.e. half the variance of the deviance is used as an estimate of the number of free parameters in the model.

    The third method calculates the parameter deviance as:

    .. math::

        p_{D} = 2 \cdot (\bar{D} - D(\bar{\theta}))

    That is, twice the complexity of that of the first method.

    Finally, the DIC is (for all cases) defined as:

    .. math::

        DIC = \bar{D} + p_{D}

    Args:
        mean_posterior_lls (ndarray): a 1d matrix containing the log likelihood for the average posterior
            point estimate. That is, the single log likelihood of the average parameters.
        ll_per_sample (ndarray): a (d, n) array with for d problems the n log likelihoods.
            This is the log likelihood per sample.

    Returns:
        dict: a dictionary containing the ``DIC_2002``, the ``DIC_2004`` and the ``DIC_Ando_2011`` information
            criterion maps.
    """
    mean_deviance = -2 * np.mean(ll_per_sample, axis=1)
    deviance_at_mean = -2 * mean_posterior_lls

    pd_2002 = mean_deviance - deviance_at_mean
    pd_2004 = np.var(ll_per_sample, axis=1) / 2.0

    return {'DIC_2002': np.nan_to_num(mean_deviance + pd_2002),
            'DIC_2004': np.nan_to_num(mean_deviance + pd_2004),
            'DIC_Ando_2011': np.nan_to_num(mean_deviance + 2 * pd_2002)}