Function to help estimate fit parameter uncertainties

[RichardWaiteSTFC]

Description

New function in ndbase to calculate the hessian matrix given a cost function and parameters. The hessian matrix evaluated at the minimum can be used to calculate the covariance matrix and hence the uncertainties on fit parameters (this is useful because not all ndbase optimizers return uncertainties).

To-Do

• [x] Use in powder fitting (need to add test)
• [x] Optionally output cost function at parameters as well as hessian

To Test

(1) Run this to compare the hessian to that calculated using the jacobian of the residuals from lsqnonlin

% make data
rng default;
fit_fun = @(x, params) params(1)*x + params(2)./x + params(3);
pars = [5,6,5];
x = linspace(0.1,2);
y = fit_fun(x, pars);
e = sqrt(y);
y = y + e.*randn(size(y));

% fit data
fcost = @(params) sum(((y-fit_fun(x, params))./e).^2);  % chisq
[pars_fit, chisq, ~] = ndbase.simplex([],fcost,pars);

figure('color', 'white');
ax =  subplot(1, 1, 1);
hold on; box on;
errorbar(ax, x, y, e, 'ok');
plot(ax, x, fit_fun(x, pars_fit), '-r');

% estimate errors on fit parameters
hess = ndbase.estimate_hessian(fcost, pars_fit);
ndof = (numel(y) - numel(pars_fit))
cov = inv(hess) * 2.0 * chisq/ndof;
perr = sqrt(diag(cov)) % errors on best fit parameters

% compare with jacobian from lsqnonlin

fresid = @(params) ((y-fit_fun(x, params))./abs(e));
[x, cost, resid_sq, exitflag, ~,~,jacobian] = lsqnonlin(fresid,pars);
reduced_cost = (cost/ndof);
cov_lsqnonlin = reduced_cost * inv(jacobian' * jacobian);
perr_lsqnonlin = sqrt(diag(cov_lsqnonlin))   

[codecov-commenter]

:warning: Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

Project coverage is 42.65%. Comparing base (1b52407) to head (479e28f).

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #198 +/- ## ========================================== + Coverage 42.37% 42.65% +0.27% ========================================== Files 240 241 +1 Lines 16184 16255 +71 ========================================== + Hits 6858 6933 +75 + Misses 9326 9322 -4 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[github-actions[bot]]

Test Results

    4 files  ± 0    116 suites  +4   4m 55s :stopwatch: +19s   731 tests +11    713 :white_check_mark: +11  18 :zzz: ±0  0 :x: ±0  2 096 runs  +44  2 060 :white_check_mark: +44  36 :zzz: ±0  0 :x: ±0 

Results for commit 479e28f8. ± Comparison against base commit 1b524075.

:recycle: This comment has been updated with latest results.

[RichardWaiteSTFC]

@mducle Thanks for the review, I have done most of the changes requested, in addition in https://github.com/SpinW/spinw/pull/198/commits/479e28f898bde4ced589a207ba617d9655049475 I now output parameter errors as a vector the same size as the input parameters (with 0 error for parameters which were fixed). I thought this was more standard/expected behaviour.

I think there are only 2 points still to be ironed out:

• Do you want me to have a smarter damping ratio for step-size?
• Do you want me to copy the parameter vector in the loop over parameters rather than altering in-place?

[mducle]

@RichardWaiteSTFC Thanks for this.

Let's keep the current auto step-size algorithm. If we find when we use it that it's too slow then we can come back and revisit this and maybe update it.

I'm happy to leave the in-place changing of the parameter vector.