The function uses the method of moments to estimate distribution parameters by minimizing
the weighted sum of squared relative errors between theoretical and sample moments
(orders 1 through 5). The optimization employs the Nelder-Mead simplex algorithm,
which is derivative-free and particularly robust for this problem.
Key implementation features: logarithmic calculations for numerical stability,
adaptive numerical integration using Simpson's rule with fallback to trapezoidal rule,
multiple random starting points to avoid local minima, decreasing weights for
higher-order moments (1.0, 0.8, 0.6, 0.4, 0.2), and automatic parameter constraint
enforcement.
Parameter Constraints:
All parameters are constrained to positive values. Additionally, family-specific
constraints are enforced: alpha and beta in (0.1, 50.0), gamma in (0.1, 10.0) for
GKw-related families or (0.1, 50.0) for Beta, delta in (0.01, 10.0), and lambda in
(0.1, 20.0).
The function will issue warnings for empty input vectors, sample sizes less than 10
(unreliable estimation), or failure to find valid parameter estimates (returns defaults).