Possible issue with weighting in `fit_linear_model`

[pabloreque]

Description

While testing the likelihood scaling for a spectroscopic observation with small flux uncertainties, I noticed that the fitted scaling coefficient seemed to be driven much more strongly than expected by points with very small err_obs values. I traced this back to the generic linear scaling routine used in the model/observation comparison, and I wanted to ask whether the current weighting in fit_linear_model() is intentional.

I noticed that fit_linear_model() applies the observational errors (here) as:

weights = 1.0 / err_obs**2

and then builds the weighted least-squares system (here) as:

A_w = A * weights[:, None]
b_w = flx_obs * weights

The system is then solved with either:

res = optimize.lsq_linear(A_free, b_w, bounds=(low, high))

or:

coeffs_free, *_ = np.linalg.lstsq(A_free, b_w, rcond=None)

Since lsq_linear / lstsq minimize the squared residual internally,

$$ \min_c \sum_i \left(A_{w,i} c - b_{w,i}\right)^2, $$

using weights = 1 / err_obs**2 in both A_w and b_w seems to make the effective weighting proportional to:

$$ \frac{1}{\mathrm{err}_{\mathrm{obs}}^4} $$

rather than the usual chi-square weighting:

$$ \frac{1}{\mathrm{err}_{\mathrm{obs}}^2}. $$

Explanation

For a single model component, the current formulation is equivalent to:

$$ A_w = \frac{\mathrm{model}}{\mathrm{err}_{\mathrm{obs}}^2}, \qquad b_w = \frac{\mathrm{obs}}{\mathrm{err}_{\mathrm{obs}}^2}. $$

Therefore, the solver minimizes:

$$ \sum_i \left(A_{w,i} c - b_{w,i}\right)^2, $$

which becomes:

$$ \sum_i \left( \frac{c , \mathrm{model}_i - \mathrm{obs}_i} {\mathrm{err}_{\mathrm{obs},i}^2} \right)^2, $$

or equivalently:

$$ \sum_i \frac{ \left(c , \mathrm{model}_i - \mathrm{obs}_i\right)^2 }{ \mathrm{err}_{\mathrm{obs},i}^4 }. $$

Why this may matter

If err_obs represents the usual 1-sigma uncertainty, the standard Gaussian chi-square for a scaled model is normally:

$$ \chi^2 = \sum_i \left( \frac{ \mathrm{obs}_i - c , \mathrm{model}_i }{ \mathrm{err}_{\mathrm{obs},i} } \right)^2. $$

To express this as a least-squares problem, the weighted system would usually be built as:

$$ A_w = \frac{ \mathrm{model} }{ \mathrm{err}_{\mathrm{obs}} }, \qquad b_w = \frac{ \mathrm{obs} }{ \mathrm{err}_{\mathrm{obs}} }, $$

which means that

weights = 1.0 / err_obs**2

should actually be

weights = 1.0 / err_obs

and therefore the usual $1 / \mathrm{err}_{\mathrm{obs}}^2$ chi-square weighting.

With the current 1 / err_obs**2 applied before the least-squares solver, points with small uncertainties receive much more weight than expected.

Question

Could you confirm whether this effective $1 / \mathrm{err}_{\mathrm{obs}}^4$ weighting is intentional in fit_linear_model() or if there is somewhere where I am mistaken?