To incorporate the heteroscedasticity, we need to rely on the heteroscedastic Gaussian process regression (implemented as VHGPR following [1]} in this work). In VHGPR, we are able to approximate (REF ) as (given \(\mathcal {D}\) )
\(f(x) & \sim \mathcal {GP}(\mu _f(x), {\rm {cov}}_f(x,x^{\prime })),  \\R(x, \omega ) & \sim \mathcal {N}(0, e^{g(x)}),  \\g(x) & \sim \mathcal {GP} (\mu _g(x), {\rm {cov}}_g(x,x^{\prime })) ,\) 
