code: def prior_GP_var_half_cauchy(y_invK_y, n_y, tau_range):
      
      tau2 = (y_invK_y - n_y * tau_range**2
              + np.sqrt(n_y**2 * tau_range**4 + (2 * n_y + 8)
                        * tau_range**2 * y_invK_y + y_invK_y**2))\
          / 2 / (n_y + 2)
      log_ptau = scipy.stats.halfcauchy.logpdf(
          tau2**0.5, scale=tau_range)
      return tau2, log_ptau
docs: Imposing a half-Cauchy prior onto the standard deviation (tau)
          of the Gaussian Process which is in turn a prior imposed over
          a function y = f(x).
          The scale parameter of the half-Cauchy prior is tau_range.
          The function returns the MAP estimate of tau^2 and
          log(p(tau|tau_range)) for the MAP value of tau^2,
          where tau_range describes the reasonable range of tau
          in the half-Cauchy prior.
          An alternative form of prior is inverse-Gamma prior on tau^2.
          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`