Problems fitting Log Pearson III in R

Question

I'd like to perform a log Pearson III fit to some data points I have. However, everytime I try to it I get errors messages I don't really know what to do about. I should perhaps add that I'm only using R since a couple of days ago, so, I'm not an expert in it.

The important code part, the part without the import stuff and so so on is this:

pIIIpars<-list(shape=1, location=1, scale=1) 

dPIII<-function(x, shape, location, scale) PearsonDS::dpearsonIII(x, shape=1, location=1, scale=1, params=pIIIpars, log=FALSE)

pPIII<-function(q, shape, location, scale) PearsonDS::ppearsonIII(q, shape=1, location=1, scale=1, params=pIIIpars, lower.tail = TRUE, log.p = FALSE)

qPIII<-function(p, shape, location, scale) PearsonDS::qpearsonIII(p, shape=1, location=1, scale=1, params=pIIIpars, lower.tail = TRUE, log.p = FALSE)

fitPIII<-fitdistrplus::fitdist(flowdata3$OEP, distr="PIII", method="mle", start=list("shape"=5000, "location"=5000, "scale"=5000))

summary(fitPIII)

plot(fitPIII)

I'm using the PearsonDS package for the definition of the Log Pearson III distribution and fitdistrplus to do the fit.

The error message I always get is this:

[1] "Error in optim(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data,  : \n  function cannot be evaluated at initial parameters\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in optim(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data,     ddistnam = ddistname, hessian = TRUE, method = meth, lower = lower,     upper = upper, ...): function cannot be evaluated at initial parameters>
Error in fitdistrplus::fitdist(flowdata3$OEP, distr = "PIII", method = "mle",  : 
  the function mle failed to estimate the parameters, 
                with the error code 100

I do unterstand the error message, it's just; if that's not the correct way to pass initial values, what is? Anyone have an idea?

Cheers, Robert

Answer 1

The following sample follows Kite (2004) and matches his results.

# Annual maximum discharge data for the St Mary River at Stillwater Nova Scotia (Kite, 2004)
# PIII fit to the logs of the discharges

StMary <- c(565,294,303,569,232,405,228,232,394,238,524,368,464,411,368,487,394,
            337,385,351,518,365,515,280,289,255,334,456,479,334,394,348,428,337,
            311,453,328,564,527,510,371,824,292,345,442,360,371,544,552,651,190,
            202,405,583,725,232,974,456,289,348,564,479,303,603,514,377,318,342,
            593,378,255,292)

LStMary <- log(StMary)

m <- mean(LStMary)
v <- var(LStMary)
s <- sd(LStMary)
g <- e1071::skewness(LStMary, type=1)

# Correct the sample skew for bias using the recommendation of 
# Bobee, B. and R. Robitaille (1977). "The use of the Pearson Type 3 and Log Pearson Type 3 distributions revisited." 
# Water Resources Reseach 13(2): 427-443, as used by Kite

n <- length(StMary)
g <- g*(sqrt(n*(n-1))/(n-2))*(1+8.5/n)

# We will use method of moment estimates as starting values for the MLE search

my.shape <- (2/g)^2
my.scale <- sqrt(v)/sqrt(my.shape)*sign(g) # modified as recommended by Carl Schwarz
my.location <- m-sqrt(v * my.shape)

my.param <- list(shape=my.shape, scale=my.scale, location=my.location)


dPIII<-function(x, shape, location, scale) PearsonDS::dpearsonIII(x, shape, location, scale, log=FALSE)
pPIII<-function(q, shape, location, scale) PearsonDS::ppearsonIII(q, shape, location, scale, lower.tail = TRUE, log.p = FALSE)
qPIII<-function(p, shape, location, scale) PearsonDS::qpearsonIII(p, shape, location, scale, lower.tail = TRUE, log.p = FALSE)

fitdistrplus::fitdist(LStMary, distr="PIII", method="mle", start=my.param)

Also note that the MLE estimates may not always be applicable. See Kite (2004, p119). He quotes Matalas and Wallis (1973) who note that if the sample skew is small then a solution is my not be possible. You can see that in the method of moments estimators because the shape parameter will blow up as g goes to zero.

Kite, G. W. (2004) Frequency and risk analyses in hydrology. Water Resources Publications

Matalas, N. C. and J. R. Wallis (1973). "Eureka! It fits a Pearson Type 3 Distribution." Water Resources Research 9(2): 281-289.

Answer 2

The above code will work if the data are positively skewed, i.e. a long right tail. But if you want to fit the data to a negatively skewed distribution, you need to multiply the my.scale by the sign of the skewness coefficient.

Refer to equations (8), (9), and (10) on page 24 of

Guidelines for determining flood flow and frequency Bulletin 17

available at https://acwi.gov/hydrology/Frequency/b17c/

So modify the line for my.scale to

my.scale <- sqrt(v)/sqrt(my.shape)*sign(g)

If you don't modify the scale, the starting values often won't lead to convergence of the fitting function.