Hello,

I am trying to fit a multivariate linear regression (Y = alpha + beta1 * X1 + beta2 * X2) wherein my 2 input variables (X1, X2) are correlated.

I have uncertainties for all the variables (Y, X1, X2). I am using the following code to estimate the values for the intercept and the slopes:

```
with pm.Model() as model:
# Priors for unknown model parameters
alpha = pm.Normal("alpha", mu=0, sigma=10)
beta1 = pm.Normal("beta1", mu=0, sigma=10)
beta2 = pm.Normal("beta2", mu=0, sigma=10)
# Prior for covariance
chol, corr, stds = pm.LKJCholeskyCov(
"chol", n=2, eta=2.0,
sd_dist=pm.Exponential.dist(1.0), compute_corr=True
)
μ = pm.Normal("μ", 0., 10., shape=2, testval=test_values)
obs = pm.MvNormal("obs", μ, chol=chol, observed=observed_df)
# Likelihood
mu = alpha + beta1 * obs[:, 0] + beta2 * obs[:, 1]
Y_obs = pm.StudentT("Y_obs", mu=mu, sigma=y_error,
observed=y, nu=5)
```

`y`

and `y_error`

are the dependent variable ‘Y’ and its uncertainty, respectively.

`observed_df`

is a DataFrame with X1 and X2 as columns.

How can I incorporate the uncertainties on X1 and X2 in my model?