This function estimates the parameters of the 4-parameter Asymmetric Exponential Power distribution given the L-moments of the data in an L-moment object such as that returned by lmoms. The relation between distribution parameters and L-moments is seen under lmomaep4. Relatively straightforward, but difficult to numerically achieve, optimization is needed to extract the parameters from the L-moments.Delicado and Goria (2008) make argue for numerical methods to use the following objective function
$$\epsilon(\alpha, \kappa, h) = \log(1 + \sum_{r=2}^4 (\hat\lambda_r - \lambda_r)^2)$$
and subsequently solve directly for $\xi$.
This objective function was chosen by Delicado and Goria because the solution surface can become quite flat for away from the minimum. The author of lmomco agrees with the findings of those authors from limited exploratory analysis and the development of the algorithms used here under the rubic of the DG method. This exploration resulted in an alternative algorithm using tabulated initial guesses described below. An evident drawback of the Delicado-Goria algorithm, is that precision in $\alpha$ is may be lost according to the observation that this parameter can be analytically computed given $\lambda_2$, $\kappa$, and $h$.
It is established practice in L-moment theory of four (and similarly three) parameter distributions to see expressions for $\tau_3$ and $\tau_4$ used for numerical optimization to obtain the two higher parameters ($\alpha$ and $h$) first and then see analytical expressions directly compute the two lower parameters ($\xi$ and $\alpha$). The author made various exploratory studies by optimizing on $\tau_3$ and $\tau_4$ through a least squares objective function. Such a practice seems to perform acceptably when compared to that recommended by Delicado and Goria (2008) when the initial guesses for the parameters are drawn from pretabulation of the relation between ${\alpha, h}$ and ${\tau_3, \tau_4}$.
Another optimization, referred to here as the A method, is available for parameter estimation using the following objective function
$$\epsilon(\kappa, h) = \sqrt{(\hat\tau_3 - \tau_3)^2 + (\hat\tau_4 - \tau_4)^2}$$
and subsequently solve directly for $\alpha$ and then $\xi$. The A method appears to perform slightly better in $\kappa$ and $h$ estimation and quite a bit better in $\alpha$ and and $\xi$ as seemingly expected because these last two are analytically computed. The objective function of the A method defaults to use of the $\sqrt{x}$ but this can be removed by setting sqrt.t3t4=FALSE.
The initial guesses for the $\kappa$ and $h$ parameters derives from a hashed environment in in file
sysdata.rda (.lmomcohash$AEPkh2lmrTable) in which the ${\kappa, h}$ pair having the minimum $\epsilon(\kappa, h)$ in which $\tau_3$ and $\tau_4$ derive from the table as well. The file SysDataBuilder.R provides additional technical details on how the AEPkh2lmrTable was generated.
The table represents a systematic double-loop sweep through lmomaep4 for
$$-3 \le \log(\kappa) \le 3, \Delta=0.05$$
and
$$-3 \le \log(h) \le 3, \Delta=0.05$$
The function will not return parameters if the following lower bounds of $\tau_4$ is not met:
$\tau_4 \ge 0.77555(|\tau_3|) - 3.3355(|\tau_3|)^2 + 14.196(|\tau_3|)^3 - 29.909(|\tau_3|)^4 + 37.214(|\tau_3|)^5 - 24.741(|\tau_3|)^6 + 6.7998(|\tau_3|)^7$. For this polynomial, the residual standard error is RSE = 0.0003125 and the maximum absolute error for $\tau_3{:}[0,1] < 0.0015$. The actual coefficients in paraep4 have additional significant figures.