Variances of dependent parameters are wrong

[henrikjacobsenfys]

All parameters are treated as independent of each other and themselves. This means that dependency expressions like a*a*a get the wrong variance, because it treats the three a as three independent parameters.

[damskii9992]

Well, how is this done in Scipp? Is this really a Scipp issue?

[henrikjacobsenfys]

Scipp does it correctly for e.g. a+a, but not for a+a.copy() (sorry if syntax is wrong, you get the point). The way we're resolving dependent parameters is via a lot of copy, so we're often not giving scipp the information to do it correctly

[damskii9992]

I don't think I understand, we do a+a with scipp under the hood . . .
We don't copy anything . . .

[henrikjacobsenfys]

We do, when we make new parameters out of old ones when e.g. adding them

[damskii9992]

That is not copying, that is creating a new parameter. And scipp does the same thing. Scipp also creates a new Scipp scalar when adding two scalars. So I still don't understand the difference . . .

[henrikjacobsenfys]

Sorry, I need to be more specific. addition usually works. But e.g. dependency expressions get evaluated term by term if I understood correctly, and here parameters are created as intermediaries, and that's where the problem can occur

[damskii9992]

But again, Scipp does EXACTLY the same thing. This is a Python thing. Any arithmetic operation is done term by term in Python. We are literally doing the same thing as scipp for variances/errors.

[damskii9992]

Okay, I've conferred with the Scipp team, specifically with @jokasimr, the variances will still be wrong if you write your dependency expression wrong. Writing a dependency as "a+a+a" instead of "3*a" is WRONG, but it SHOULD give the same result.

In Scipp it does, because of a simple and naive check they have. We can replicate the Scipp functionality by checking if the other in an arithmetic operation is itself, and then skipping the "_convert_unit" call, which messes up with this Scipp check.

However, if a dependency is written as "A + (1 - A)" we will still get the wrong result, as will Scipp, due to a fundamental limitation of Python.

In conclusion:
We can fix this "bug" by adding a small check in the arithmetic operations code, which will fully replicate Scipp functionality, HOWEVER this does not mean that we treat variances correctly for non-reduced equations. So we SHOULD add a note to dependent parameters specifiying that they should ALWAYS fully reduce their expressions.

[jokasimr]

I don't think you will get the right result if you write a + a + a because a + a will be written to a new variable and then the final addition new_var + a will not recognize that new_var is correlated with a.

This might be what you meant already, I just wanted to clarify to there are no misunderstandings.

Generally I think it is best to assume that scipp operations will not handle correlations correctly.
But in some rare and simple cases they do.

[jokasimr]

If local linear uncertainty propagation is insufficient for you, and you need "global" linear uncertainty propagation, and you can evaluate the expressions you need the variances of symbolically, it might be worth looking into sympy to do that.

Maybe using something like:

import sympy as sp


def linearized_variance_function(expr, symbols):
    variances = [
        sp.Symbol(f"var_{x}", nonnegative=True)
        for x in symbols
    ]

    var_expr = sum(
        sp.diff(expr, x)**2 * v
        for x, v in zip(symbols, variances)
    )

    return sp.lambdify(
        (*symbols, *variances),
        var_expr,
        modules="numpy",
    )


x, y, z = sp.symbols("x y z")

f_var = linearized_variance_function(
    x * y + z,
    [y, x],  # z is assumed exact
)

f_var(y, x, var_y, var_x)