2
$\begingroup$

I am trying to use the package ‘fitdistrplus’ in R to fit one non standard distribution to my data set. I am trying to copy the methods the package creators used in their tutorial for specifying a Gumbel distribution but unfortunately I am not able to do it.

(The package documentation can be found in this pdf.)

So my problem is that I have my CDF function and my density function which are not standard :

$ \begin{align*} H_{G}(x)&=\frac{1}{P(1)}\left[ \frac{\Phi^{3}\left( \frac{x-\mu}{\sigma}\right)}{6}-\frac{\delta}{2} \Phi^{2}\left( \frac{x-\mu}{\sigma}\right)+\left( \frac{\delta^{2}}{2} + \beta\right) \Phi\left( \frac{x-\mu}{\sigma}\right)\right]\\ h_{G}(x)&= \frac{1}{\sigma P(1)}\left[ \frac{\left[ \Phi^{2}\left( \frac{x-\mu}{\sigma}\right)\right] }{2}-\delta \Phi\left( \frac{x-\mu}{\sigma}\right)+\left( \frac{\delta^{2}}{2} + \beta\right) \right] \phi\left( \frac{x-\mu}{\sigma}\right) \end{align*}$

with $ P(1)= \left[ \frac{1}{6}-\frac{\delta}{2}+\left( \frac{\delta^{2}}{2} + \beta\right) x\right] $ and \begin{equation*} \Phi(\dfrac{x-\mu}{\sigma}) =\frac{1}{\sqrt{2\pi}}\int^{x}_{-\infty}e^{-\frac{\left(t-\mu\right) ^{2}}{2\sigma^{2}}}dt \qquad \text{and} \qquad \phi(\dfrac{x-\mu}{\sigma}) =\frac{1}{\sqrt{2\pi}}e^{-\frac{\left(x-\mu\right) ^{2}}{2\sigma^{2}}}. \end{equation*}

And here is some of my code:

###### I define my density as in the gumbel ######

ds <- function(u,a,e,mu,sigma){
  P = function(a,e){ (( (1/6)*(1^3) ) - ((a/2)*(1^2)) + (((((a)^2)/2) + e)*1)) }
  D = function(u,mu,sigma){ (1/((sigma)*sqrt(2*pi)))*exp(-(((u-mu)^2)/(2*(sigma^2)))) }
  K = function(u,a,e){ (((1/2)*(u^2))- (a*u) +(((a^2)/2)+e)) }
  H = function(u,mu,sigma){ pnorm(u,mu,sigma) }
  Fprim = function(u,a,e,mu,sigma){ (1/P(a,e))*(D(u,mu,sigma))*(K(H(u,mu,sigma),a,e)) }
  return(Fprim)
} 

###### I define my CDF function as in the gumbel ######

ps <- function(u,a,e,mu,sigma) {
  S   = function(u) (((1/6)*(u^3)) - ((a/2)*(u^2)) + 
                    ( ((((a)^2)/2)+e)*u))/(((1/6)*(1^3)) - ((a/2)*(1^2)) + 
                    (((((a)^2)/2)+e)*1))
  cdf = S(pnorm(u, mu, sigma))
  return(cdf)
}

###### I define my quantile function as in the gumbel ######
qs <- function(u,a,e,mu,sigma) optimize(function(z) (ps(z)-u)^2, c(-10,10))$minimum


###### The MLE estimation using fitdist and mledist ######
fgu <- fitdist(X, "s", start=list(a=0.035, e=0.005, mu=-0.52, sigma=1))

mledist(X, "s", start=list(a=0.035, e=0.005, mu=-0.52, sigma=1))
$\endgroup$
  • $\begingroup$ Are these functions from some parametric distribution? Or how did you get them? And what happens/what error do you get when you run your code? $\endgroup$ – HFBrowning Apr 9 '14 at 21:52
  • $\begingroup$ My density is from composition of two distribution a polynomial of degre 3 and a normal, this is not in any familly of distribution. $\endgroup$ – Lea Apr 9 '14 at 22:15
  • $\begingroup$ The erro : [1] "Error in log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)))) : \n non-numeric argument to mathematical function\n"attr(,"class")[1] "try-error" attr(,"condition") <simpleError in log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)))): non-numeric argument to mathematical function> Error in fitdist(X, "s", start = list(a = 0.035, e = 0.005, mu = -0.52, : the function mle failed to estimate the parameters, with the error code 100 $\endgroup$ – Lea Apr 9 '14 at 22:17
  • $\begingroup$ In the tutorial there is one example with the gumbel distribution dgumbel <- function(x,a,b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b)) pgumbel <- function(q,a,b) exp(-exp((a-q)/b)) qgumbel <- function(p,a,b) a-b*log(-log(p)) f1c <- fitdist(x1,"gumbel",start=list(a=10,b=5)) print(f1c) plot(f1c) $\endgroup$ – Lea Apr 9 '14 at 22:20
  • $\begingroup$ Okay, a couple of things. The package is designed to fit data to parametric distributions so while there may be some way to coerce your functions into working, I'm not sure what it would be. (Maybe someone else knows but I do think it is outside the intent of the package.) A bigger question is, however: why do you think your data fit this distribution? Is there any reason why your empirical distribution wouldn't do just as well? Can you provide a graph maybe? $\endgroup$ – HFBrowning Apr 9 '14 at 22:31
1
$\begingroup$

The primary issue seems to be in the definition of ds. You would like to return a density value, but right now you return a function. The change needed in ds is to just remove that internal function definition:

Fprim = (1/P(a,e))*(D(u,mu,sigma))*(K(H(u,mu,sigma),a,e))

There is also an issue in qs, where the ps function needs to be passed the distribution parameters. The change needed in qs is to add the distribution parameters to the function call:

qs <- function(u,a,e,mu,sigma) optimize(function(z) (ps(z,a,e,mu,sigma)-u)^2, c(-10,10))$minimum

Next, there is a potential issue of having starting values for optimization that produce values that are too unlikely. If values end up getting rounded to 0 or 1 along the way, then optim grinds to a halt. So, the start values need to be good.

For example, with those changes, this seems to work:

set.seed(12)
X <- rnorm(100, 0, 1)

###### The MLE estimation using fitdist and mledist ######
fgu <- fitdist(X, "s", start=list(a=0.035, e=0.005, mu=-0.52, sigma=1))

mledist(X, "s", start=list(a=0.035, e=0.005, mu=-0.52, sigma=1))

It looks like there is a lot of sensitivity to the start values, in terms of whether or how well fitdist will converge.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.