def prior_GP_var_inv_gamma(y_invK_y, n_y, tau_range):
    """ Imposing an inverse-Gamma prior onto the variance (tau^2)
        parameter of a Gaussian Process, which is in turn a prior
        imposed over an unknown function y = f(x).
        The inverse-Gamma prior of tau^2, tau^2 ~ invgamma(shape, scale)
        is described by a shape parameter alpha=2 and a scale parameter
        beta=tau_range^2. tau_range describes the reasonable range of
        tau in the inverse-Gamma prior.
        The data y's at locations x's are assumed to follow Gaussian Process:
        f(x, x') ~ N(0, K(x, x') / 2 tau^2), where K is a kernel
        function defined on x. For n observations, K(x1, x2, ..., xn) is
        an n by n positive definite matrix.
        Given the prior parameter tau_range, number of observations
        n_y, and y_invK_y = y * inv(K) * y',
        the function returns the MAP estimate of tau^2 and
        the log posterior probability of tau^2 at the MAP value:
        log(p(tau^2|tau_range)).
        This function is written primarily for BRSA but can also
        be used elsewhere. y in this case corresponds to the log of
        SNR in each voxel. GBRSA does not rely on this function.
        An alternative form of prior is half-Cauchy prior on tau.
        Inverse-Gamma prior penalizes for both very small and very
        large values of tau, while half-Cauchy prior only penalizes
        for very large values of tau.
        For more information on usage, see description in BRSA class:
        `.BRSA`

        See also: `.prior_GP_var_half_cauchy`

    Parameters
    ----------
    y_invK_y: float
        y * inv(K) * y^T, where y=f(x) is a vector of observations
        of unknown function f at different locations x.
        K is correlation matrix of f between different locations, based
        on a Gaussian Process (GP) describing the smoothness property
        of f. K fully incorporates the form of the kernel
        and the length scale of the GP, but not the variance of the GP
        (the purpose of this function is to estimate the variance).
    n_y: int, number of observations
    tau_range: float,
        The reasonable range of tau, the standard deviation of the
        Gaussian Process imposed on y=f(x). tau_range is parameter
        of the inverse-Gamma prior. Say, if you expect the standard
        deviation of the Gaussian process to be around 3, tau_range
        can be set to 3.
        The smaller it is, the more penalization is imposed
        on large variation of y.
    Returns
    -------
    tau2: The MAP estimation of tau^2 based on the prior on tau
        and y_invK_y.
    log_ptau: log(p(tau)) of the returned tau^2 based on the
        inverse-Gamma prior.
    """
    alpha = 2
    tau2 = (y_invK_y + 2 * tau_range**2) / (alpha * 2 + 2 + n_y)
    log_ptau = scipy.stats.invgamma.logpdf(
        tau2, scale=tau_range**2, a=2)
    return tau2, log_ptau
Imposing an inverse-Gamma prior onto the variance (tau^2)
        parameter of a Gaussian Process, which is in turn a prior
        imposed over an unknown function y = f(x).
        The inverse-Gamma prior of tau^2, tau^2 ~ invgamma(shape, scale)
        is described by a shape parameter alpha=2 and a scale parameter
        beta=tau_range^2. tau_range describes the reasonable range of
        tau in the inverse-Gamma prior.
        The data y's at locations x's are assumed to follow Gaussian Process:
        f(x, x') ~ N(0, K(x, x') / 2 tau^2), where K is a kernel
        function defined on x. For n observations, K(x1, x2, ..., xn) is
        an n by n positive definite matrix.
        Given the prior parameter tau_range, number of observations
        n_y, and y_invK_y = y * inv(K) * y',
        the function returns the MAP estimate of tau^2 and
        the log posterior probability of tau^2 at the MAP value:
        log(p(tau^2|tau_range)).
        This function is written primarily for BRSA but can also
        be used elsewhere. y in this case corresponds to the log of
        SNR in each voxel. GBRSA does not rely on this function.
        An alternative form of prior is half-Cauchy prior on tau.
        Inverse-Gamma prior penalizes for both very small and very
        large values of tau, while half-Cauchy prior only penalizes
        for very large values of tau.
        For more information on usage, see description in BRSA class:
        `.BRSA`

        See also: `.prior_GP_var_half_cauchy`

    Parameters
    ----------
    y_invK_y: float
        y * inv(K) * y^T, where y=f(x) is a vector of observations
        of unknown function f at different locations x.
        K is correlation matrix of f between different locations, based
        on a Gaussian Process (GP) describing the smoothness property
        of f. K fully incorporates the form of the kernel
        and the length scale of the GP, but not the variance of the GP
        (the purpose of this function is to estimate the variance).
    n_y: int, number of observations
    tau_range: float,
        The reasonable range of tau, the standard deviation of the
        Gaussian Process imposed on y=f(x). tau_range is parameter
        of the inverse-Gamma prior. Say, if you expect the standard
        deviation of the Gaussian process to be around 3, tau_range
        can be set to 3.
        The smaller it is, the more penalization is imposed
        on large variation of y.
    Returns
    -------
    tau2: The MAP estimation of tau^2 based on the prior on tau
        and y_invK_y.
    log_ptau: log(p(tau)) of the returned tau^2 based on the
        inverse-Gamma prior.