data {
    int N;
    int p;
    vector[N] y;
    matrix[N,p] x;
}
parameters {
    real b0;
    vector[p] d_beta;
    vector[p] d_a;
    vector<lower=0>[p] d_b;
    real<lower=0, upper=1> tau;
    real<lower=0> sigma;
    real<lower=0> nu;
}
transformed parameters {
    vector[p] beta;
    beta = d_beta .* d_a .* sqrt(d_b) * tau;
}
model {
    y ~ normal(b0 + x * beta, sigma);
    // HORSESHOE PRIOR:
    d_beta ~ normal(0, 1);
    d_a ~ normal(0, 1);
    d_b ~ inv_gamma(0.5 * nu, 0.5 * nu);
    // tau ~ uniform(0, 1); // uniform
    // priors on noise distribution:
    sigma ~ normal(0, 10);
    // prior on shape parameter for horseshoe prior
    nu ~ normal(0, 20);
}

CXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function