CODEWITHSUNDEEP
rdistributiongamma-distribution
Pablo Sanhueza
Pablo Sanhueza

Reputation: 21

Problems fitting a log pearson III with negative scale

I'd like to perform a log Pearson III fit to some data points I have in R. I followed this guidelines link. But i encounter a problem when my skewness (g) is negative (and of course the parameter "scale" is negative too, because the "sign(g)" in the computation of scale). The distribution from "FAdist" does not work with negative scale parameter, i need this for starting values for the "fitdist" (from fitdistrplus). Some pages say that the parameters shape and scale is only positive in a pearson III (or generalize gamma distribution) and other dont, i run out of ideas. The code is:

library(fitdistrplus)
library(FAdist)
library(e1071) 
#data
df <-(92.8, 53.2, 112.0, 164.0, 132.0,  69.9, 140.0,  48.3, 123.0 ,24.6, 179.0,  55.1,  31.3,  17.0, 111.0,  35.4, 133.0, 505.0, 303.0, 121.5, 203.0, 198.0, 250.0, 232.0, 185.0, 222.0, 191.0, 238.0,  53.0, 121.0, 106.4, 347.3, 186.4,  89.1, 131.4 ,53.2, 252.6)
# log of df
df<-log(df)
#Pearson 3 Sample moments
 m <- mean(df)
v <- var(df)
s <- sd(df)
g <- e1071::skewness(df, type=1)
n <- length(df)
#Pearson 3 Parameter estimation
my.shape <- (2/g)^2
my.scale <- sqrt(v)/sqrt(my.shape)*sign(g) # modified as recommended by Carl Schwarz, this is negative =(
my.thres <- m-(my.shape*my.scale)
# All parameter together
my.param <- list(shape=my.shape, scale=my.scale, thres=my.thres )
# fit dist from the "fitdistrplus" and "lgamma3" from "FAdist"
q.p3 <- fitdist(data = caudales,distr = "lgamma3", start = my.param)

Give me the following

Error in fitdist(data = df, distr = "lgamma3", start = my.param, : the function mle failed to estimate the parameters, with the error code 100

Upvotes: 1

Views: 572

Answers (1)

Pablo Sanhueza
Pablo Sanhueza

Reputation: 21

I already solve the problem. I used the package "PearsonDS" instead of "FAdist" to compute log pearson III, why? because the latter accept negative values of skewness (g) and in consecuence negative values of "scale" parameter (if u look at the documentation of PearsonDS in Pearson III, they add an absolute value in the argument). I had to do a little changes to the equation of Pearson III to work right with package "fitdistrplus" (to fit distribution to sample data):

The final code is right this (IMPORTANT : DATA MUST BE THE ORIGINAL VALUES NOT THE LOG(DATA)), BUT MEAN,VAR AND G THEY NEED TO BE COMPUTED WITH LOG(DATA):

library(FAdist); library(fitdistrplus)
library("e1071")
library(PearsonDS)
load("data")
## PP for precipitation
pp <- as.numeric(data$pp)
## IMPORTANT : DATA MUST BE THE ORIGINAL VALUES NOT THE LOG(DATA)
data <- (pp)
m <- mean(log(data))
v <- var(log(data))
s <- sd(log(data))
g <- e1071::skewness(log(data), type=1)
n <- length(data)
#g <- g*(sqrt(n*(n-1))/(n-2))*(1+8.5/n)
my.shape <- (2/g)^2
my.scale <- (sqrt(v/my.shape))*sign(g) # modified as recommended by Carl Schwarz
my.thres <- m-my.shape*my.scale
#load functions
dlPIII<-function(x, shape, location, scale) {PearsonDS::dpearsonIII(log(x), shape, location, scale, log=FALSE)/x}
plPIII<-function(q, shape, location, scale) {PearsonDS::ppearsonIII(log(q), shape, location, scale, lower.tail = TRUE, log.p = FALSE)}
qlPIII<-function(p, shape, location, scale) {exp(PearsonDS::qpearsonIII(p, shape, location, scale, lower.tail = TRUE, log.p = FALSE))}
# Fit the distribution
f.pl3 <- fitdist(data, distr="lPIII", method="mle", start=list(shape=my.shape,location=my.thres ,scale=my.scale )

Upvotes: 1

Related Questions