clouizos
commented
5 years ago Apologies for the late reply; I just noticed the issues here. You can empirically check that the decomposition we have there leads to the same distribution as sampling a half-Cauchy scale for a Gaussian via the following script:
import numpy as np
import matplotlib.pyplot as plt
def ecdf(x):
xs = np.sort(x)
ys = np.arange(1, len(xs) + 1) / float(len(xs))
return xs, ys
n_rand_samples = 100000
k = 1 # global scale of the horseshoe prior
# sample a half-Cauchy r.v. C^+(z; 0, k), i.e. sample a standard Cauchy with
# the given scale 'k' and then take the absolute value
z = np.abs(k * np.random.standard_cauchy(size=(n_rand_samples,)))
# sample the Gaussian N(w; 0, z^2)
w = z * np.random.randn(n_rand_samples)
# calculate and plot the (empirical) CDF of the samples
xw, yw = ecdf(w)
plt.plot(xw, yw, label='Cauchy HS')
# sample the inverse gamma variable IGamma(b; 0.5, 1) by simply taking the inverse
# of the gamma variable Gamma(0.5, 1), since if x ~ Gamma(k, theta) then 1/x ~ IGamma(k, 1/theta)
beta = 1./np.random.gamma(.5, size=(n_rand_samples))
# sample the gamma variable G(a; 0.5, k^2)
alpha = k**2 * np.random.gamma(.5, size=(n_rand_samples))
# multiply and take the square root for z
z = np.sqrt(beta * alpha)
# sample the Gaussian N(w; 0, z^2)
w = z * np.random.randn(n_rand_samples)
# calculate and plot the (empirical) CDF of the samples
xw, yw = ecdf(w)
plt.plot(xw, yw, label='Alternative HS')
plt.legend()
plt.xlim([-10, 10])
plt.show()
@clouizos @KarenUllrich
How is the equation derived according to
