How to propagate measurement uncertainty in predictors *and* responses for multidimensional, non-parametric regression (and software to do it)?

Question

Background

Errors-in-variables models are defined as:

regression models that account for measurement errors in the independent variables. In contrast, standard regression models assume that those regressors have been measured exactly, or observed without error; as such, those models account only for errors in the dependent variables, or responses.

Whether data is experimental or computational (e.g. via physics-based simulation), uncertainty in the independent variables is common. This is particularly relevant in materials science where uncertainty exists in the processing, structure, and property measurements. Additionally, this uncertainty is not necessarily constant for all parameters (a temperature sensor might have higher uncertainty at the limits of the temperature range, various aspects of a material's structure may have been determined by separate instruments, property measurements are data-mined from several sources, etc.). See also Uncertainty Quantification. I'm sure this is relevant for other disciplines too (e.g. physics, chemistry, biology, engineering).

Existing Methods

What I'm not looking for

• linear regression models (e.g. orthogonal regression/total least squares OrthogonalLineFit[], linortfitn(), odr(), odregress()).
• response-uncertainty-only regression models
  • Nonlinear parametric fitting with a specified vector of response uncertainties using NonlinearModelFit[]
  • Non-parametric fitting with a specified scalar response uncertainty using fitrgp(). This allows for the specification of a vector of "smoothness lengths", one per predictor, but not for specifying uncertainty of predictors

What I am looking for

A method that caters to multidimensional, non-parametric regression with propagated measurement uncertainty in predictors and responses (i.e. uncertainty propagation, not just weighting the points) and preferably software that goes along with it (Mathematica, MATLAB, Python, R, Stan, etc.). Multidimensional refers to the predictors (i.e. responses are scalar).

Ideally, it would be something that is explained and can be performed via:

mdl = some_function(X, Xsd, y, ysd)
ypred, ysd2, cov = mdl(X2, Xsd2)

Inputs

1. Matrix/Array of predictors (X $X$)
2. Matrix/Array of predictor uncertainties (Xsd $\sigma_X$)
3. Vector of response values (y $y$)
4. Vector of response uncertainties (ysd $\sigma$)
5. Matrix/Array of new predictors (X2 $X_2$)
6. (Optional) Matrix/Array of new predictor uncertainties (Xsd2 $\sigma_{X2}$)

Outputs

1. Vector of new responses (ypred $y_{pred}$)
2. Vector of new response uncertainties (ysd2 $\sigma_2$)
3. (Optional) predictor covariance matrix (cov $\Sigma_2$)

Assume uncertainty means standard deviation of a normal distribution. I'm fine with starting out with simpler assumptions before moving onto e.g. other distributions.

Updates

• This seems to be supported through BoTorch, a package based in Python (see tutorial), and to later be exposed in Meta's Adaptive Experimentation platform. May eventually try this and flesh it out into an answer. If someone else would like to try it out and provide a MWE, this is also welcome.

Related

• SE: Non-linear fitting with uncertainty in dependent and independent variable
• Errors-in-variables multivariate polynomial regression (R)
• Including model uncertainty in non-linear least squares minimization
• How do I specify the uncertainty of a parameter when using NonlinearModelFit?
• How to add uncertainty to your neural network
• BLiTZ — A Bayesian Neural Network library for PyTorch

Answer 1

One of the more interesting choices in R is rstan, where you could code this up yourself in the Stan modeling language (which tends to be amazing in that it can produce inference for models that we used to be unable to do for a long time). However, getting started can be a little challenging and it sounds like you'd like a higher level interface.

That could be the brms package, that uses the R model syntax and in the background generates Stan code. Via that detour (generate the Stan code via brms in R and then use the generated Stan code with pystan - or any of the other Stan-tie-ins in other languages such as MathematicaStan or MatlabStan - you can then also use it in Stan). In brms there's the me() function for predictors with measurement errors and a quite nice range of modeling options for non-linear models (it depends a bit of what exactly you are thinking about whether that's covered) and it also supports Gaussian processes, but I'm not sure from what you described to what extent you can fit it all together the way you want (if not you may be best off looking at the - usually quite well-chosen / efficient - Stan code that brms generates and having to fit it together exactly the way you want manually in Stan).