# Estimate the Parameters of the 4-Parameter Asymmetric Exponential Power Distribution 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. If the $\tau_3$ of the distribution is zero (symmetrical), then the distribution is known as the Exponential Power (see `lmrdia46`). Delicado and Goria (2008) argue for numerical methods to use the following objective function [REMOVE_ME]$$ \epsilon(\alpha, \kappa, h) = \log(1 + \sum_{r=2}^4 (\hat\lambda_r - \lambda_r)^2)\mbox{,} [REMOVE_ME_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 (Asquith) method, is available for parameter estimation using the following objective function [REMOVE_ME]$$ \epsilon(\kappa, h) = \sqrt{(\hat\tau_3 - \tau_3)^2 + (\hat\tau_4 - \tau_4)^2}\mbox{,} [REMOVE_ME_2]$$ and subsequently solve directly for $\alpha$ and then $\xi$. The A method appears to perform 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 (Asquith, 2014). 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 SysDataBuilder01.R provides additional technical details on how the `AEPkh2lmrTable` was generated. The table represents a systematic double-loop sweep through `lmomaep4` for [REMOVE_ME]$$ \kappa \mapsto \{-3 \le \log(\kappa) \le 3, \Delta\log(\kappa)=0.05\}\mbox{,} [REMOVE_ME_2]$$ and [REMOVE_ME]$$ h \mapsto \{-3 \le \log(h) \le 3, \Delta\log(h)=0.05\}\mbox{.} [REMOVE_ME_2]$$ The function will not return parameters if the following lower (estimated) bounds of $\tau_4$ are 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. However, the argument `snap.tau4`, if set, will set $\tau_4$ equal to the prediction from the polynomial. This value of $\tau_4$ should be close enough numerically to the boundary because the optimization is made using a log-transformation to ensure that $\alpha$, $\kappa$, and $h$ remain in the positive domain---though the argument `nudge.tau4` is provided to offset $\tau_4$ upward just incase of optimization problems. ## Description 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. If the $\tau_3$ of the distribution is zero (symmetrical), then the distribution is known as the Exponential Power (see `lmrdia46`). Delicado and Goria (2008) 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)\mbox{,} $$ 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 (Asquith) 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}\mbox{,} $$ and subsequently solve directly for $\alpha$ and then $\xi$. The A method appears to perform 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 (Asquith, 2014). 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 SysDataBuilder01.R provides additional technical details on how the `AEPkh2lmrTable` was generated. The table represents a systematic double-loop sweep through `lmomaep4` for $$ \kappa \mapsto \{-3 \le \log(\kappa) \le 3, \Delta\log(\kappa)=0.05\}\mbox{,} $$ and $$ h \mapsto \{-3 \le \log(h) \le 3, \Delta\log(h)=0.05\}\mbox{.} $$ The function will not return parameters if the following lower (estimated) bounds of $\tau_4$ are 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. However, the argument `snap.tau4`, if set, will set $\tau_4$ equal to the prediction from the polynomial. This value of $\tau_4$ should be close enough numerically to the boundary because the optimization is made using a log-transformation to ensure that $\alpha$, $\kappa$, and $h$ remain in the positive domain---though the argument `nudge.tau4` is provided to offset $\tau_4$ upward just incase of optimization problems. ```r paraep4(lmom, checklmom=TRUE, method=c("A", "DG", "ADG"), sqrt.t3t4=TRUE, eps=1e-4, checkbounds=TRUE, kapapproved=TRUE, snap.tau4=FALSE, nudge.tau4=0, A.guess=NULL, K.guess=NULL, H.guess=NULL, ...) ``` ## Arguments - `lmom`: An L-moment object created by `lmoms` or `vec2lmom`. - `checklmom`: Should the L-moments be checked for validity using the `are.lmom.valid` function. Normally this should be left as the default and it is very unlikely that the L-moments will not be viable (particularly in the $\tau_4$ and $\tau_3$ inequality). However, for some circumstances or large simulation exercises then one might want to bypass this check. - `method`: Which method for parameter estimation should be used. The A or DG methods. The ADG method will run both methods and retains the salient optimization results of each but the official parameters in `para` are those from the A method. Lastly, all minimization is based on the `optim` function using the Nelder--Mead method and default arguments. - `sqrt.t3t4`: If true and the method is A , then the square root of the sum of square errors in $\tau_3$ and $\tau_4$ are used instead of sum of square differences alone. - `eps`: A small term or threshold for which the square root of the sum of square errors in $\tau_3$ and $\tau_4$ is compared to to judge good enough for the alogrithm to set the `ifail` on return in addition to convergence flags coming from the `optim` function. Note that `eps` is only used if the A or ADG methods are triggered because the other method uses the scale parameter which in reality could be quite large relative to the other two shape parameters, and a reasonable default for such a secondary error threshold check would be ambiguous. - `checkbounds`: Should the lower bounds of $\tau_4$ be verified and if sample $\hat\tau_3$ and $\hat\tau_4$ are outside of these bounds, then `NA` are returned for the solutions. - `kapapproved`: Should the Kappa distribution be fit by `parkap` if $\hat\tau_4$ is below the lower bounds of $\tau_4$? This fitting is only possible if `checkbounds` is true. The Kappa and AEP4 overlap partially. The AEP4 extends $\tau_4$ above Generalized Logistic and Kappa extends $\tau_4$ below the lower bounds of $\tau_4$ for AEP4 and extends all the way to the theoretical limits as used within `are.lmom.valid`. - `snap.tau4`: A logical to snap the $\tau_4$ upwards to the lower boundary if the given $\tau_4$ is lower than the boundary described in the polynomial. - `nudge.tau4`: An offset to the snapping of $\tau_4$ intended to move $\tau_4$ just above the lower bounds in case of optimization problems. (The absolute value of the nudge is made internally to ensure only upward adjustment by an addition operation.) - `A.guess`: A user specified guess of the $\alpha$ parameter to provide to the optimization of any of the methods. This argument just superceeds the simple initial guess of $\alpha = 1$. - `K.guess`: A user specified guess of the $\kappa$ parameter to supercede that derived from the `.lmomcohash$AEPkh2lmrTable` in file sysdata.rda . - `H.guess`: A user specified guess of the $h$ parameter to supercede that derived from the `.lmomcohash$AEPkh2lmrTable` in file sysdata.rda . - `...`: Other arguments to pass. ## Returns An `list` is returned. - **type**: The type of distribution: `aep4`. - **para**: The parameters of the distribution. - **source**: The source of the parameters: paraep4 . - **method**: The method as specified by the `method`. - **ifail**: A numeric failure code. - **ifailtext**: A text message for the failure code. - **L234**: Optional and dependent on method DG or ADG . Another `list` containing the optimization details by the DG method along with the estimated parameters in `para_L234`. The _234 is to signify that optimization is made using $\lambda_2$, $\lambda_3$, and $\lambda_4$. The parameter values in `para` are those only when the DG method is used. - **T34**: Optional and dependent on method A or ADG . Another `list` containing the optimization details by the A method along with the estimated parameters in `para_T34`. The _T34 is to signify that opimization is being conducted using $\tau_3$ and $\tau_4$ only. The parameter values in `para` are those by the A method. The values for `ifail` or produced by three mechanisms. First, the convergence number emanating from the `optim` function itself. Second, the integer 1 is used when the failure is attributable to the `optim` function. Third, the interger 2 is a general attempt to have a singular failure by sometype of `eps` outside of `optim`. Fourth, the integer 3 is used to show that the parameters fail against a parameter validity check in `are.paraep4.valid`. And fifth, the integer 4 is used to show that the sample L-moments are below the lower bounds of the $\tau_4$ polynomial shown here. Additional and self explanatory elements on the returned list will be present if the Kappa distribution was fit instead. ## Author(s) W.H. Asquith ## References Asquith, W.H., 2014, Parameter estimation for the 4-parameter asymmetric exponential power distribution by the method of L-moments using R: Computational Statistics and Data Analysis, v. 71, pp. 955--970. Delicado, P., and Goria, M.N., 2008, A small sample comparison of maximum likelihood, moments and L-moments methods for the asymmetric exponential power distribution: Computational Statistics and Data Analysis, v. 52, no. 3, pp. 1661--1673. ## See Also `lmomaep4`, `cdfaep4`, `pdfaep4`, `quaaep4`, `quaaep4kapmix` ## Examples ```r # As a general rule AEP4 optimization can be CPU intensive ## Not run: lmr <- vec2lmom(c(305, 263, 0.815, 0.631)) plotlmrdia(lmrdia()); points(lmr$ratios[3], lmr$ratios[4], pch=16, cex=3) PAR <- paraep4(lmr, snap.tau4=TRUE) # will just miss the default eps FF <- nonexceeds(sig6=TRUE) plot(FF, quaaep4(FF, PAR), type="l", log="y") lmomaep4(PAR) # 305, 263, 0.8150952, 0.6602706 (compare to those in lmr) ## End(Not run) ## Not run: PAR <- list(para=c(100, 1000, 1.7, 1.4), type="aep4") lmr <- lmomaep4(PAR) aep4 <- paraep4(lmr, method="ADG") print(aep4) # ## End(Not run) ## Not run: PARdg <- paraep4(lmr, method="DG") PARasq <- paraep4(lmr, method="A") print(PARdg) print(PARasq) F <- c(0.001, 0.005, seq(0.01,0.99, by=0.01), 0.995, 0.999) qF <- qnorm(F) ylim <- range( quaaep4(F, PAR), quaaep4(F, PARdg), quaaep4(F, PARasq) ) plot(qF, quaaep4(F, PARdg), type="n", ylim=ylim, xlab="STANDARD NORMAL VARIATE", ylab="QUANTILE") lines(qF, quaaep4(F, PAR), col=8, lwd=10) # the true curve lines(qF, quaaep4(F, PARdg), col=2, lwd=3) lines(qF, quaaep4(F, PARasq), col=3, lwd=2, lty=2) # See how the red curve deviates, Delicado and Goria failed # and the ifail attribute in PARdg is TRUE. Note for lmomco 2.3.1+ # that after movement to log-exp transform to the parameters during # optimization that this "error" as described does not appear to occur. print(PAR$para) print(PARdg$para) print(PARasq$para) ePAR1dg <- abs((PAR$para[1] - PARdg$para[1])/PAR$para[1]) ePAR2dg <- abs((PAR$para[2] - PARdg$para[2])/PAR$para[2]) ePAR3dg <- abs((PAR$para[3] - PARdg$para[3])/PAR$para[3]) ePAR4dg <- abs((PAR$para[4] - PARdg$para[4])/PAR$para[4]) ePAR1asq <- abs((PAR$para[1] - PARasq$para[1])/PAR$para[1]) ePAR2asq <- abs((PAR$para[2] - PARasq$para[2])/PAR$para[2]) ePAR3asq <- abs((PAR$para[3] - PARasq$para[3])/PAR$para[3]) ePAR4asq <- abs((PAR$para[4] - PARasq$para[4])/PAR$para[4]) MADdg <- mean(ePAR1dg, ePAR2dg, ePAR3dg, ePAR4dg) MADasq <- mean(ePAR1asq, ePAR2asq, ePAR3asq, ePAR4asq) # We see that the Asquith method performs better for the example # parameters in PAR and inspection of the graphic will show that # the Delicado and Goria solution is obviously off. (See Note above) print(MADdg) print(MADasq) # Repeat the above with this change in parameter to # PAR <- list(para=c(100, 1000, .7, 1.4), type="aep4") # and the user will see that all three methods converged on the # correct values. ## End(Not run) ```

paraep4 function