More support for non-normal maximum likelihood estimation

[jagerber48]

This "idea" may actually be a question if it turns out there are just details about lmfit that I'm missing. This post is motivated by this stack exchange question which Newville and myself have provided answers to.

Historically, I think lmfit arose out of the need for a more unified API to scipy least squares minimization routines so much of the code is oriented towards least squares fitting. However, least squares minimization is just a specific type of maximum likelihood estimation. As I point out in my answer to the question, if the data points are normally distributed, the sum-of-squares-of-residuals cost function is equal to -2 times the log-likelihood function, so minimizing the sum-of-squares-of-residuals corresponds to maximizing the log-likelihood.

---

The idea I'm proposing here is to bake this realization about the relationship between maximum likelihood estimation (MLE) and least squares estimation (LSE) into the lmfit source code in such a way that lmfit becomes more user friendly for users wanting to do maximum likelihood estimation for non-normally distributed data.

Right now, if I wanted to use lmfit to do maximum likelihood estimation on non-normally distributed data I would do the following:

• Write a residual function that returns the negative of the log-likelihood function for my fit parameters given my observed data
• Multiply the covariance matrix I get back from lmfit by 1/2. I think also multiply all standard errors etc. by a factor of 1/sqrt(2).

---

Maybe the easiest way to make lmfit more user friendly for non-normal log-likelihoods would be to create an MLE mode in lmfit where

• The cost function is maximized rather than minimized
• The covariance matrix is calculated as the inverse of the negative Hessian of the cost function (rather than two times the inverse of the Hessian of the cost function). This covariance is used for downstream calculation of errors/confidence intervals.

Additionally, for user-convenience, some common likelihood cost function calculators could be provided. For example, factory functions that

• accept sequences counts data and a counts model depending on parameters theta and returns a log-likelihood function using scipy.poisson.logpmf
• accepts sequences of numbers of successes and number of trials and a binomial probability model that depends on parameters theta and returns log-likelihood function using scipy.binomial.logpfm

These changes are, of course, not strictly necessary because users cal just write these residual functions themselves (just remember to include a negative sign because lmfit only does minimization and won't ever provide the negative sign on its own.

---

As a note, since least squares estimation is a specific case of maximum likelihood estimation, if the MLE machine described above was set up then it could be possible, under the hood, to reconfigure the least squares machinery in lmfit to take advantage of the maximum likelihood code. It would essentially amount to redirecting the least-squares cost function to include a factor -2 so that it become the maximum likelihood cost function.

[jagerber48]

Actually, the easiest way to get this functionality right now without any changes to lmfit would be to just make sure my "log-likelihood" cost function is actually -2 times the actual log-likelihood. It's easy to do, but annoying because it's sort of reverse shoe-horning the general maximum likelihood estimation into the least squares case.

[newville]

@jagerber48 Well, lmfit supports using different functions for "reducing" a residual array ((data-model)*weights) to a scaler "cost function to be minimized.

The default is "sum of squares", but "negentropy" and "’neglogcauchy" are available as built-in options. Or a custom callable can be used.

I'm unsure what to make of your concern about a factor of 2. How is a scale factor an issue?