Estimate the Parameters of the 4-Parameter Asymmetric Exponential Power Distribution

Share:

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.

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.

Usage

1
2
3
paraep4(lmom, checklmom=TRUE, method=c("A", "DG", "ADG"),
        sqrt.t3t4=TRUE, eps=1e-4, checkbounds=TRUE, kapapproved=TRUE,
        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 τ_4 and τ_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 τ_3 and τ_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 τ_3 and τ_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 τ_4 be verified and if sample \hatτ_3 and \hatτ_4 are outside of these bounds, then NA are returned for the solutions.

kapapproved

Should the Kappa distribution be fit by parkap if \hatτ_4 is below the lower bounds of τ_4? This fitting is only possible if checkbounds is true. The Kappa and AEP4 overlap partially. The AEP4 extends τ_4 above Generalized Logistic and Kappa extends τ_4 below the lower bounds of τ_4 for AEP4 and extends all the way to the theoretical limits as used within are.lmom.valid.

A.guess

A user specified guess of the α parameter to provide to the optimization of any of the methods. This argument just superceeds the simple initial guess of α = 1.

K.guess

A user specified guess of the κ 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.

Value

An R 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 R 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 λ_2, λ_3, and λ_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 R 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 τ_3 and τ_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 τ_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

 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)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.