I have the following model with bayesian lasso prior, where y
is of dimension Nx1
, X
is Nx3
and beta
is 3x1
.
with pm.Model() as blasso:
# Blasso parameters
lam = 2
sig_1 = pm.math.sqrt(pm.InverseGamma(name="sig_1", alpha=3, beta=1, shape=()))
sig_2 = pm.math.sqrt(pm.InverseGamma(name="sig_2", alpha=3, beta=1, shape=()))
sig_3 = pm.math.sqrt(pm.InverseGamma(name="sig_3", alpha=3, beta=1, shape=()))
# Slope
beta_1= pm.Laplace(name="beta_1", mu=0.0, b=sig_1/lam, shape=())
beta_2= pm.Laplace(name="beta_2", mu=0.0, b=sig_2/lam, shape=())
beta_3= pm.Laplace(name="beta_3", mu=0.0, b=sig_3/lam, shape=())
# Model error
eps = pm.Gamma(name="eps", alpha=9.0, beta=4.0, shape=())
# Model mean
y_hat = pm.Deterministic(name="y_hat", var=beta_1*X[:, 0] + beta_2*X[:, 1] + beta_3*X[:, 2])
# Likelihood
y_like = pm.Normal(name="y_like", mu=y_hat, sigma=eps, observed=y)
However, this is very inefficiently written, as I have to specify each parameter separately.
Can I write the model in a more efficient way, something like this, without changing the implications?:
with pm.Model() as blasso:
# Blasso parameters
lam = 2
sig = pm.math.sqrt(pm.InverseGamma(name="sig", alpha=3, beta=1, shape=3))
# Slope
beta= pm.Laplace(name="beta", mu=0.0, b=sig/lam, shape=3)
# Model error
eps = pm.Gamma(name="eps", alpha=9.0, beta=4.0, shape=())
# Model mean
y_hat = pm.Deterministic(name="y_hat", var=pm.math.dot(X, beta))
# Likelihood
y_like = pm.Normal(name="y_like", mu=y_hat, sigma=eps, observed=y)
I am not sure whether this is the same as the algorithm takes longer for sampling and does not converge with the same number of draws.