alternate fitting statistics with Model interface

[ggggggggg]

I often fit data from single photon sensing spectrometers with few counts, meaning that I often have bins that have zero or 1 count in them along with bins that have 100s or thousands of counts. The sherpa fitting package for chandra offers the "cash statistic" for this case: https://cxc.cfa.harvard.edu/sherpa/ahelp/cash.html and offers a general method for modifying the fit statistic: https://cxc.cfa.harvard.edu/sherpa/statistics/

As far as I can tell, the closest thing in lmfit is the reduce_fcn argument to Minimizer. The docs for reduce_fcn say:

Function to convert a residual array to a scalar value for the scalar minimizers. Optional values are (where r is the residual array):

None : sum-of-squares of residual (default)

= (r*r).sum()

’negentropy’ : neg entropy, using normal distribution

= rho*log(rho).sum()`, where rho = exp(-r*r/2)/(sqrt(2*pi))

’neglogcauchy’ : neg log likelihood, using Cauchy distribution

= -log(1/(pi*(1+r*r))).sum()

callable : must take one argument (r) and return a float.

So reduce_fcn allows the same modifications to fit statistics, but appears limited in the following ways:

1. Not available through the very desirable Model interface.
2. Only works with scalar minimizers.

I think that allowing the use of reduce_fcn or similar through the Model interface and with non-scalar minimizers would be a valuable addition to LMfit. If you like the idea, what would be a good path to do so?

[newville]

@ggggggggg Hm, this suggests that our GitHub Issue Templates are not being used, and so need to be updated.

We are happy to have discussions about lmfit. And we want to use GitHub Issues only and precisely as a bug tracker.

Converting to a discussion.

[newville]

@ggggggggg I think that allowing reduce_fcn in Model.fit (and ModelResult.fit) would help here and be pretty easy to do -- it would just get forwarded to minimize.

But: That would push the residual to something other than Least-Squares. And using those "reduce" functions are painful because they remove the array-ness of the residuals and replace it with a scalar. That makes the fit much harder.

If I understand the Cash statistic from the Sherpa docs (and maybe I am not, so please correct if wrong), we would really want to minimize

2*(data*log(model) - model).

I think it would be better to be able to add an option to Model._residual to apply that "log(model)" bias. But, I'll also say that I do not know how common that is -- it seems sort of strange to me to weight the data by the model, and not also weight the model itself. I could understand

(data-model)*log(model)

but (data*log(model) - model) seems like it has "dimensional analysis" problems to me. But, maybe I am missing something.....