The following code is adapted from an example invented by @KieranDelamotte.
We try to fit a Gumbel distribution on a sample that is actually obtained from a mixture of a Gumbel and a Fréchet distribution.
GumbelFactory
does not manage to find a sensible set of parameters to try to fit the sample.
WeibullMaxFactory
, on the other hand, manages to find a Weibull distribution that kind of fits, but it has crazy parameters… because it is actually very close to a Gumbel. In fact, when put into a GEV, the GEV \xi parameter is 0: it is a Gumbel distribution.
The bottom line is that WeibullMaxFactory
manages to find a better Gumbel
than GumbelFactory
.
The difference is that WeibullMaxFactory
maximizes (by default) the likelihood whereas GumbelFactory
fits the first and second order moments to the empirical moments of the sample. Likelihood maximization is more robust to (this particular) distribution misspecification than the method of moments.
# Try to fit a WeibullMax on the sample.
# WeibullMaxFactory maximizes the likelihood.
weibullmax = ot.WeibullMaxFactory().build(sample)
print("WeibullMaxFactory result: ", weibullmax)
# Express the estimated WeibullMax as a GEV
gev = ot.GeneralizedExtremeValue(weibullmax)
print("GEV actual distribution: ", gev.getActualDistribution())
# Recompute the GEV actual distribution
# This WeibullMax turns out to be a Gumbel
# because the GEV xi parameter is zero.
print("GEV xi parameter = ", gev.getXi())
gev.setParameter(gev.getParameter())
print("More sensible GEV actual distribution: ", gev.getActualDistribution())
# Try to directly estimate a Gumbel distribution.
# GumbelFactory uses the method of moments.
gumbel = ot.GumbelFactory().build(sample)
print("GumbelFactory result: ", gumbel)
graph = ot.HistogramFactory().build(sample).drawPDF()
graph.add(gev.drawPDF())
graph.add(gumbel.drawPDF())
graph.setColors(["green", "red", "blue"])
graph.setLegends(["Sample", "MLE", "Moments"])
graph.setTitle("Estimated Gumbel distribution")
view = viewer.View(graph)
axes = view.getAxes()
_ = axes[0].set_xlim(-20.0, 20.0)
Output:
WeibullMaxFactory result:
WeibullMax(beta = 4.25531e+06, alpha = 1.45411e+06, gamma = 4.25531e+06)
GEV actual distribution:
WeibullMax(beta = 4.25531e+06, alpha = 1.45411e+06, gamma = 4.25531e+06)
GEV xi parameter = -6.877037385642863e-07
More sensible GEV actual distribution:
Gumbel(beta = 2.92639, gamma = 2.5121)
GumbelFactory result: Gumbel(beta = 15.1936, gamma = -3.88787)