human
Solve the following python programming problem given the following code:
# Construct a variational linear model
class VariationalLinearModel(hb.model.Model):
    def setUp(self):
        # data should be stored in hb.param.Data class
        self.X = hb.param.Data(X)
        self.Y = hb.param.Data(Y)
        
        # Variational parameters 
        self.a = hb.variationals.Normal(shape=[1])
        self.b = hb.variationals.Normal(shape=[1])
        self.c = hb.variationals.Normal(shape=[1])
        
        # To access the posterior samples, 
        # we define op that returns samples of the fit.
        # As described in Henbun_structure.ipynb, tf_mode is used here.
        with self.tf_mode():
            self.fit = self.a + self.b*self.X
        
    @hb.model.AutoOptimize()
    def elbo(self):
        """
        This method returns the sum of log-likelihood and log-prior.
        """
        # log likelihood
        log_lik = hb.densities.gaussian(self.Y, self.fit, tf.exp(self.c))
        # KL() methods automatically gather KL for all the variational parameters.
        KL = self.KL()
        return tf.reduce_sum(log_lik) - KL

vlinear_model = VariationalLinearModel()
Given the above code, generate some markdown instructions for it.
assistant
#Here is the Python program you requested.

Variational linear model.
Here, we construct a linear model within Bayesian framework.
The likelihood is assumed as Gaussian,
$$
p(y_i | a, b, x_i) = \mathcal{N}(y_i | a + b*x_i, c)
$$
where $\sigma$ is the variance parameter.
We assume weak prior for $a$, $b$ and $\sigma$,
$$
p(a) = \mathcal{N}(a|0,1) \
p(b) = \mathcal{N}(b|0,1) \
p(c) = \mathcal{N}(c|0,1) \
$$
In this model, we approximate the posterior distribution by independent Gaussian, i.e.
$$
p(a|\mathbf{x},\mathbf{y}) \sim q(a) = \mathcal{N}(a| \mu_a, \sigma_a) \
p(b|\mathbf{x},\mathbf{y}) \sim q(b) = \mathcal{N}(b| \mu_b, \sigma_b) \
p(c|\mathbf{x},\mathbf{y}) \sim q(c) = \mathcal{N}(c| \mu_c, \sigma_c) \
$$
where $\mu_a, \mu_b, \mu_c$ and $\sigma_a, \sigma_b, \sigma_c$ are the variational parameters to be optimized.
By this approximation, the free energy (evidence lower bound: ELBO) to be maximized is
$$
\mathrm{ELBO}(a, b, c) = 
\int q(a)q(b)q(c) \log p(\mathbf{y}|\mathbf{x}, a, b, c)\; \mathrm{d}a\;\mathrm{d}b\;\mathrm{d}c \
- \mathrm{KL}[q(a)||p(a)]  \
- \mathrm{KL}[q(b)||p(b)]  \
- \mathrm{KL}[q(c)||p(c)]
$$
Construct the variational model