'IndexError: tuple index out of range' error when using Confidence via fit.est_errors method

[DougBurke]

The following sequence, fitting a simple 1D dataset with errors using chi-square statistics, fails when it comes to using Confidence to estimate the errors. It appears to be due to the two-phase approach (fit with c1 frozen then with c1 thawed). I've put details and test scripts in
https://gist.github.com/DougBurke/2e5000ed458d261537907b569f2afc02

create Data1D
create Polynom1D
fit c0 and c2 parameters
fit c0, c1, c2 parameters
call est_errors with Confidence set as the estmethod

This results in the error:

IndexError: tuple index out of range

This is not a recent change, as I see the same behaviour in CIAO 4.8 (with Python 2.7) and master (with Python 3.5). It also appears not to happen using the sherpa.ui layer.

Note that the error case (python fail.py 1 in the gist) reports this when evaluating the errors before it errors out

   mdl.c0         -9.38489    
   mdl.c1         -2.41692    
   mdl.c2         0.478273    
mdl.c0 lower bound: -2.91751
mdl.c0 upper bound: 2.91751
mdl.c2 -: WARNING: The confidence level lies within (-2.419416e+00, -2.416915e+00)
mdl.c2 lower bound: -0.00125014
mdl.c2 +: WARNING: The confidence level lies within (-2.416911e+00, -2.416915e+00)
mdl.c2 upper bound: 2.44169e-06

whereas the successful run reports

   mdl.c0         -9.38489    
   mdl.c1         -2.41692    
   mdl.c2         0.478273    
mdl.c0 lower bound: -2.91751
mdl.c0 upper bound: 2.91751
mdl.c1 lower bound: -0.250889
mdl.c2 lower bound: -0.0312677
mdl.c1 upper bound: 0.250889
mdl.c2 upper bound: 0.0312677

What is interesting here is that the c2 bounds are very different in the two runs: ~0.03 in the second case but 0.001 and lower in the failing case. My guess is that this is part of the problem that leads to the IndexError (i.e. something internal is really messed up).

[DougBurke]

I think the issue might be that the fit object only thinks there's two thawed parameters, even though there are three. Just before calling est_errors (which will fail) we have

In [1]: f.thaw_indices
Out[1]: (0, 2)

In [2]: [p.name for p in f.model.pars if not p.frozen]
Out[2]: ['c0', 'c1', 'c2']

If I call f.est_errors() now, it errors out. However, if I manually change the thaw_indices to include 0, 1, and 2 then the est_errors call succeeds:

In [1]: f.thaw_indices = (0, 1, 2)

In [2]: conferrs = f.est_errors()
mdl.c0 lower bound:	-2.91751
mdl.c1 lower bound:	-0.250889
mdl.c2 lower bound:	-0.0312677
mdl.c0 upper bound:	2.91751
mdl.c1 upper bound:	0.250889
mdl.c2 upper bound:	0.0312677

In [3]: print(conferrs.format())
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   mdl.c0           -9.38489     -2.91751      2.91751
   mdl.c1           -2.41692    -0.250889     0.250889
   mdl.c2           0.478273   -0.0312677    0.0312677

[DougBurke]

So, if you get this, one work around appears to be

f.thaw_indices = tuple(i for i, p in enumerate(f.model.pars) if not p.frozen)