Is the regularised horseshoe prior example correct?

#1

In the example model page, the regularised horseshoe prior is given as follows as per Piironen & Vehtari (2017) Electronic Journal of Statistics:

regularized_horseshoe <- function (tau = 1,  c = 1, dim = NULL) {
  stopifnot(c > 0)
  lambda <- cauchy(0, 1, truncation = c(0, Inf), dim = dim)
  lambda_tilde <- (c^2 * lambda^2) / (c^2 + tau^2 * lambda^2)
  sd <- tau ^ 2 * lambda_tilde ^ 2
  normal(0, sd, dim = dim)
}

However, Piironen & Vehtari (2017) seems to denote tau ^ 2 * lambda_tilde ^ 2 as variance sd^2. If this is true, should the penultimate line be sd <- tau * lambda_tilde instead?

#2

Yeah, that does like like it is misspecified - thanks for spotting that!
In fact judging by equation 2.8, lambda_tilde looks like it should actually be called lambda_tilde_squared, so the square root of that needs to be taken before multiplying by tau to get sd

We should check the definition of the non-regularised one too. I’ve opened an issue for this here: https://github.com/greta-dev/greta/issues/406

Tagging @dirmeier who I think looked at these a while back.

#3

Ah, thanks for pointing out another line that needs change. Looking forward to the update! Meanwhile, below is a quick patch for myself, which hopefully is correct.

regularized_horseshoe <- function (tau = 1,  c = 1, dim = NULL) {
    stopifnot(c > 0)
    lambda <- cauchy(0, 1, truncation = c(0, Inf), dim = dim)
    lambda_tilde_squared <- (c^2 * lambda^2) / (c^2 + tau^2 * lambda^2)
    # sd <- sqrt(tau ^ 2 * lambda_tilde_squared)
    sd <- tau * sqrt(lambda_tilde_squared)   # same as above
    normal(0, sd, dim = dim)
}
#4

A colleague suggested not to square c and lambda immediately before taking their square root for numerical reasons. We have changed our version to:

regularized_horseshoe <- function (tau = 1,  c = 1, dim = NULL) {
    stopifnot(c > 0)
    stopifnot(tau > 0)
    lambda <- cauchy(0, 1, truncation = c(0, Inf), dim = dim)
    sd <- tau * c * lambda / sqrt(c^2 + tau^2 * lambda^2)
    normal(0, sd, dim = dim)
}

Hopefully we got it right.

All the best towards v0.4!!!

1 Like
#5

Thanks! We are looking forward to getting a soft release of 0.4.0 out soon!