e-baumer/standard_precip

Parameter Estimation (L-Moments or MLE?)

Closed this issue · 1 comments

Hello, the comments in the fit_distribution section of base_sp.py state that distribution parameters are being fit using Maximum Likelihood Estimation (MLE). However, the function being used is lmom_fit() from the lmoments3 package, which I believe is not, in fact, MLE, but a method of moments approach. Shouldn't the parameter fit be done using the fit function, instead of lmom_fit, available in either the lmoments3 or base stats modules? These functions do in fact use MLE as the estimation process, which is not the case in lmom_fit. Parameter values will be slightly different between (L-moment) lmom_fit and (MLE) fit.

# MLE option using base stats
params = stats.gamma.fit(data, floc=0)

# MLE Option using lmoments3 (same as above)
from lmoments3 import distr
params = distr.gam.fit(data, floc=0)

# As opposed to L-Moments approach using lmoments3
params = distr.gam.lmom_fit(data)

@njdepsky You are absolutely correct. Thank you for catching this. My initial intent was to have the distribution fit using MLE; however, the codes used at NCAR and in the R Cran repository to calculate SPI all use moments to estimate the parameters of the distribution.

In order to match results from these programs, I used the L-moments. I will update the docs and maybe add an option to fit the distribution using MLE.