Estimating parameters of Log Pearson III distribution in R

Question

How to estimate the Parameters of "Log Pearson III" distribution and how to generate the random nos. using these estimated parameters in R?

Kindly guide

Akshata

Answer 1

The Pearson type III distribution is a shifted gamma distribution. A random variable $X$ has log-Pearson III distribution if $\log(X)$ has Pearson III distribution.

Then, the simulation part is simple, simulate from a gamma distribution, shift it (translate it) and then take the exp() of these values.

Now, in order to estimate the parameters, you can either obtain the density of the log-Pearson III distribution or to use its relationship with the shifted-gamma. This is, if you have a sample $(x_1,...,x_n)\sim \mbox{log-Pearson III}$, then $(\log(x_1),...,\log(x_n))\sim \mbox{Pearson III}$ which is also a shifted-gamma and this is easier to estimate.

The following R code implements both, the simulation and the estimation of the parameters of this distribution.

rm(list=ls())
# Simulation from a log-Pearson III
set.seed(1000)
rlogpearson <- function(n,a,b,c) return( exp(rgamma(n,shape=a,rate=b) - c) )

data<- rlogpearson(1000,3,3,5)

hist(data)

# Transformation of the data to obtain a shifted gamma

datat <- log(data)

# - log-likelihoood

ll <- function(par){
if(par[1]>0 &  par[2]>0 & par[3]> -min(datat)) return( -sum(dgamma(datat+par[3],shape=par[1],rate=par[2],log=TRUE))  )
else return(Inf)
}

# optimisation step

optim(c(3,3,5),ll)

# MLE

optim(c(3,3,5),ll)$par

As you can see, the estimators are close to the theoretical values.

Edit

The following code shows a function that calculates the quantiles of this distribution. It uses the relationship with the gamma distribution.

qlogpearson <- function(p,a,b,c) return( exp(qgamma(p,shape=a,rate=b) - c) )

qlogpearson(0.5,3,3,5)

If you want to use the estimated parameters, just plug them in this function

param <- optim(c(3,3,5),ll)$par
qlogpearson(0.5,param[1],param[2],param[3])

Another approximation can be obtained by simulating a large sample and calculating the empirical quantiles as follows

samp<- rlogpearson(10000,param[1],param[2],param[3])

quantile(samp,0.5)

Edit II

The CDF can also be easily calculated using that

$$F(x)=P(X<x)= P(\log(X)+c<\log(x)+c),$$

and using that $\log(X)+c$ is a gamma random variable. In R

plogpearson <- function(x,a,b,c) return(pgamma(log(x)+c,shape=a,rate=b))