K = dims(X)[2]
alpha ~ normal(0, 10)
beta::vector[K] ~ normal(0, 5)
sigma ~ exponential(1)
y ~ normal(alpha + X * beta, sigma)

mu_a ~ normal(0, 5)
sigma_a ~ exponential(1)
sigma_y ~ exponential(1)
alpha::vector[J] ~ normal(mu_a, sigma_a)
y ~ normal(alpha[county], sigma_y)

beta ~ normal(0, 1)
alpha ~ normal(0, 5)
y ~ poisson_log(alpha + beta * x)

sigma ~ exponential(1)
scale = l2_norm(y)
mu ~ normal(0, scale)
y ~ normal(mu, sigma)

sq(x::vector[n]) = begin
  out = rep_vector(0., n)
  for i in 1:n
    out[i] = x[i]^2
  end
  out
end

cube(x::real) = x^3

l2_norm(x::vector[n])::real = sqrt(sum(x .* x))

standardize(x::vector[n]) = (x - mean(x)) / sd(x)

softplus(x::real)::real = log1p(exp(x))

lognormal_qf(p::real, mu::real, sigma::real)::real = exp(mu + sigma * inv_Phi(p))