How to deal with heteroscedasticity in GLM regression


I’m currently working on developing a Bayesian model to estimate certain linear regression parameters, and I’ve encountered a challenge regarding how to model the variance of the dependent variable (y) as it changes with the independent variable (x). I understand that it’s common practice to use relative residuals by normalizing the observed data within the likelihood function. However, I’m facing a difficulty in PyMC where I can’t directly normalize the observed data and then return it to its original form afterward, right?

Do you have any suggestions on how I could address the heteroscedasticity using PyMC?

These are my previous questions, provided for context:

Question 1
Question 2
I still working in the same model (I know, long time working in the samething…)

Thank you!

Your variance can also be a function of the predictors such as in y = pm.Normal("y", mu=..., sigma=pm.math.exp(b0 + b1 * x), observed=data), just like the mean is usually.