def METHOD_NAME(X, Y):
    D_X = X.shape[1]
    # sample from horseshoe prior
    lambdas = numpyro.sample("lambdas", dist.HalfCauchy(jnp.ones(D_X)))
    tau = numpyro.sample("tau", dist.HalfCauchy(jnp.ones(1)))
    # note that this reparameterization (i.e. coordinate transformation) improves
    # posterior geometry and makes NUTS sampling more efficient
    unscaled_betas = numpyro.sample("unscaled_betas", dist.Normal(0.0, jnp.ones(D_X)))
    scaled_betas = numpyro.deterministic("betas", tau * lambdas * unscaled_betas)
    # compute mean function using linear coefficients
    mean_function = jnp.dot(X, scaled_betas)
    # observe data
    numpyro.sample("Y", dist.Bernoulli(logits=mean_function), obs=Y)