Towards `bayem`

[TTitscher]

Some time ago, we decided to pack this VB algorithm into a pypi packaged named bayem.

For the usage, I have the following API in mind that is strongly inspired by scipy.optimize

import bayem

n0 = bayem.Gamma(...)
p0 = bayem.MVN(...)

result = bayem.solve(fun=model, param0=p0, noise0=n0, jac=None) # uses central differences Jacobian
# or 
result = bayem.solve(fun=model, param0=p0, noise0=n0, jac=model.jacobian, tolerance=..., callback=...)

bayem.pair_plot(result, show=True)
bayem.trace_plot(result, labels=...)

Also, I'd like to split the code (at least) into

• distributions.py contains MVN and Gamma
• vb.py contains the actual algorithm
• visualzation.py contains some plots

and then automatically load those when import bayem, so that you use bayem.Gamma (see above) instead of bayem.distributions.Gamma.

• Is the name OK?
• API OK? bayem.vba instead of solve?
• Split into files OK?
• I'd like to slowly step back from the project and like to also mention another author (or multiple) here. Who is interested? @ajafarihub @aradermacher @joergfunger ?

[aradermacher]

I agree to all :) and I would be honored to step in here.

[ajafarihub]

Hi @TTitscher , thanks for coming up with the nice API.

Just a couple of detailed stuff that I think are quite needed to be handled: 1- keys for noises and parameters and model outputs should be consistent. Maybe a check is needed somewhere. 2- which noise models are prescribed and which are to be identified as well. So, maybe a _flag input argument for that is needed. 3- as you also suggested, maybe it is good to specify which kind of Bayesian inference is used (like via 'bayem.vba'). So far, it is mostly only VB and also a specific solver behind it (e.g. no ARD, no ML, ... ). But we will probably extend it to other slight alternatives in terms of the algorithm and assumptions.

So, the following script could be a way to go.

n0 = bayem.Gamma(key='noise_model_1', to_identify=False, ...)
## ...
bayes_setup = VBSetup(...) # can be subclass of a more general one like BayesianSetup
result = bayem.solve(fun=model, param0=p0, noise0=n0, jac=model.jacobian, setup=bayes_setup) # in there, we first check consistency of inputs and then execute the proper algorithm based on bayes_setup object.
## ...

[ajafarihub]

I also would like to stress a general remark. IMO, it is generally more transparent if we distinguish btw. model and forward model. To my eye, "model" has nothing to do with any calibration/inference stuff, so, it is just a (numerical) model which represents a physical phenomenon for us, thus, having its own inputs and outputs. The model can be developed by someone who has totally no attention to any calibration and data. On the other side, "forward model" is something specifically designed for integrating the data into the model at some point. This incorporation of data is very often at the very end (when we say: model_error = model_response - data), but it can be the case that we integrate data into a forward model as input; e.g. the FEMU-F setup with identification of dU instead of U. So, I would be more happy if we make this distinction in the code (even in the user-level).

[TTitscher]

Thanks for your comments! I really like that to_indentify flag directly at the gamma distribution.

n0 = bayem.Gamma(key='noise_model_1', to_identify=False, ...)

Regarding

IMO, it is generally more transparent if we distinguish btw. model and forward model.

I like to see this repository as the implementation of a mathematical algorithm that basically performs something like a statistical least squares fit. The minimal information for the algorithm to work is the model error (and its derivative) and that is what the user should provide. So I actually like to stress it as a feature (of this very special implementation) that we do not distinguish between model and model error. IMO your concerns are still valid and addressed in probeye.

[TTitscher]

Further refinements of the interface:

result = bayem.solve(fun=f, x0=prior, noise0=noise0...) # solves with CDF jacobian
result = bayem.solve(fun=f, x0=prior, jac=j) # solves with provided "j"
result = bayem.solve(fun=f, x0=prior, jac=True) # indicates that f returns a tuple (f, jac)

Then, the format of f. The input is always a (numpy) vector, but IMO practical output choices would be:

• np.ndarray of shape (N, 1)
• list[np.ndarray] where each entry corresponds to a noise group
• dict[some_key: np.ndarray] where each entry corresponds to a noise group

Note that the format of f, jac and noise0 must be consistent.

Any thoughts on that?

[ajafarihub]

As far as I see, some solvers like the current version of VB needs the model-error (you called it, f, I think) as a dictionary, while, others might need it as a vector. My suggestion would be to put a concatenate=False flag in the model-error (default value of Flase implies outputing as dictionary) and adjust it depending on which solver is called. For example:

def concatenate_dict_me(dict_e):
    err = np.array([])
    for e in dict_e.values():
        err = np.concatenate((err, e.flatten()))
    return err