How to define your own distribution for fitdistr function in R using lmomco function
I would like to define my own distributions to use with fitdistrplus to now match my monthly rainfall data indicated as "month". I am using the "lmomco" function to help me identify distributions, but I cannot get it to work. For example, I define the distribution of the generalized extreme value (gev) as follows:
dgev<-pdfgev #functions which are included in lmomco
pgev<-cdfgev
qgev<-quagev
Since "fitdistrplus" needs a "start" argument, which consists of the initial parameter values for the desired distribution, I estimate these initial values as follows:
lmom=lmoms(month,nmom=5) #from lmomco package
para=pargev(lmom, checklmom=TRUE)
Now I will finally try to use the "fitdist" function to put the "month" in the gev distribution like:
fitgev <- fitdist(month, "gev", start=para[2]) #fitdistrplus
I am getting an error similar to the one below. It doesn't matter which distribution I define with "lmomco", I get the same error. Can anyone give me a hint as to what I am doing wrong? Thank!
fitgev <- fitdist(month, "gev", start=para[2])
[1] "Error in dgev(c(27.6, 97.9, 100.6, 107.3, 108.5, 109, 112.4, 120.9, 137.8, : \n unused arguments (para.xi = 196.19347977195, para.alpha = 91.9579520442104, para.kappa = -0.00762962879097294)\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in dgev(c(27.6, 97.9, 100.6, 107.3, 108.5, 109, 112.4, 120.9, 137.8, 138.4, 144.7, 156.8, 163.1, 168.9, 169.1, 171.4, 176.1, 177.1, 178.8, 178.9, 187.2, 190.2, 190.5, 190.8, 191.2, 193.1, 195.2, 198.5, 199.8, 201.7, 206.9, 213.4, 220.7, 240, 253.5, 254.5, 256.1, 256.4, 257.5, 258.3, 261.5, 263.7, 264.7, 279.1, 284.2, 313.1, 314.7, 319.4, 321.6, 328.9, 330.1, 332.2, 358.3, 366.8, 367.9, 403.5, 424.1, 425.9, 457.3, 459.7, 468, 497.1, 508.5, 547.1), para.xi = 196.19347977195, para.alpha = 91.9579520442104, para.kappa = -0.00762962879097294): unused arguments (para.xi = 196.19347977195, para.alpha = 91.9579520442104, para.kappa = -0.00762962879097294)>
Error in fitdist(month, "gev", start = para[2]) :
the function mle failed to estimate the parameters,
with the error code 100
source to share
fitdist
expects a density / distribution function with named arguments.
library("lmomco")
library("fitdistrplus")
## reproducible:
month <- c(27.6, 97.9, 100.6, 107.3, 108.5,
109, 112.4, 120.9, 137.8)
Setting:
lmom <- lmoms(month,nmom=5) #from lmomco package
para <- pargev(lmom, checklmom=TRUE)
dgev <- pdfgev #functions which are included in lmomco
pgev <- cdfgev
fitgev <- fitdist(month, "gev", start=para[[2]])
## Error in mledist(data, distname, start, fix.arg, ...) :
## 'start' must specify names which are arguments to 'distr'
It turns out we need to override a dgev
few extra bits of plumbing to make us happy fitdist
and pdfgev
:
dgev <- function(x,xi,alpha,kappa) {
pdfgev(x,list(type="gev",para=c(xi,alpha,kappa),source="pargev"))
}
fitgev <- fitdist(month, "gev", start=para[[2]])
## Fitting of the distribution ' gev ' by maximum likelihood
## Parameters:
## estimate Std. Error
## xi -25.587734 NA
## alpha 75.009842 NA
## kappa 1.593815 NA
source to share