Description Usage Arguments Value Author(s) References See Also Examples

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) argue for numerical methods to use the following objective function

*ε(α, κ, h) = \log(1 + ∑_{r=2}^4 (\hatλ_r - λ_r)^2)\mbox{,}*

and subsequently solve directly for *ξ*. 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 *α* is may be lost according to the observation that this parameter can be analytically computed given *λ_2*, *κ*, and *h*.

It is established practice in L-moment theory of four (and similarly three) parameter distributions to see expressions for *τ_3* and *τ_4* used for numerical optimization to obtain the two higher parameters (*α* and *h*) first and then see analytical expressions directly compute the two lower parameters (*ξ* and *α*). The author made various exploratory studies by optimizing on *τ_3* and *τ_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 *\{α, h\}* and *\{τ_3, τ_4\}*.

Another optimization, referred to here as the “A” (Asquith) method, is available for parameter estimation using the following objective function

*ε(κ, h) = √{(\hatτ_3 - τ_3)^2 + (\hatτ_4 - τ_4)^2}\mbox{,}*

and subsequently solve directly for *α* and then *ξ*. The “A” method appears to perform better in *κ* and *h* estimation and quite a bit better in *α* and and *ξ* as seemingly expected because these last two are analytically computed (Asquith, 2014). The objective function of the “A” method defaults to use of the *√{x}* but this can be removed by setting `sqrt.t3t4=FALSE`

.

The initial guesses for the *κ* and *h* parameters derives from a hashed environment in in file

‘sysdata.rda’ (.lmomcohash$AEPkh2lmrTable) in which the *\{κ, h\}* pair having the minimum *ε(κ, h)* in which *τ_3* and *τ_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

*κ \mapsto \{-3 ≤ \log(κ) ≤ 3, Δ\log(κ)=0.05\}\mbox{,}*

and

*h \mapsto \{-3 ≤ \log(h) ≤ 3, Δ\log(h)=0.05\}\mbox{.}*

The function will not return parameters if the following lower bounds of *τ_4* are not met:

*τ_4 ≥ 0.77555(|τ_3|) - 3.3355(|τ_3|)^2 + 14.196(|τ_3|)^3 - 29.909(|τ_3|)^4 + 37.214(|τ_3|)^5 - 24.741(|τ_3|)^6 + 6.7998(|τ_3|)^7*. For this polynomial, the residual standard error is RSE = 0.0003125 and the maximum absolute error for *τ_3{:}[0,1] < 0.0015*. The actual coefficients in `paraep4`

have additional significant figures.

1 2 3 |

`lmom` |
An L-moment object created by |

`checklmom` |
Should the L-moments be checked for validity using the |

`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 |

`sqrt.t3t4` |
If true and the method is “A”, then the square root of the sum of square errors in |

`eps` |
A small term or threshold for which the square root of the sum of square errors in |

`checkbounds` |
Should the lower bounds of |

`kapapproved` |
Should the Kappa distribution be fit by |

`A.guess` |
A user specified guess of the |

`K.guess` |
A user specified guess of the |

`H.guess` |
A user specified guess of the |

`...` |
Other arguments to pass. |

An **R** `list`

is returned.

`type` |
The type of distribution: |

`para` |
The parameters of the distribution. |

`source` |
The source of the parameters: “paraep4”. |

`method` |
The method as specified by the |

`ifail` |
A numeric failure code. |

`ifailtext` |
A text message for the failure code. |

`L234` |
Optional and dependent on method “DG” or “ADG”. Another |

`T34` |
Optional and dependent on method “A” or “ADG”. Another |

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 *τ_4* polynomial shown here.

Additional and self explanatory elements on the returned list will be present if the Kappa distribution was fit instead.

W.H. Asquith

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.

`lmomaep4`

, `cdfaep4`

, `pdfaep4`

, `quaaep4`

, `quaaep4kapmix`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 | ```
# As of version 1.6.2, it is felt that in spirit of CRAN CPU
# reduction that the intensive operations of paraep4() should
# be kept a bay. However, the following examples are useful.
## 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-Goria failed
# and the ifail attribute in PARdg is TRUE.
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-Goria solution is obviously off.
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)
``` |

lmomco documentation built on May 29, 2017, 6:34 p.m.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.