Package ‘ﬁtdistrplus’ I am trying to use the package  ‘ﬁtdistrplus’ 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))

 A: 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.
