title: R fitting Frechet distribution: replicating the answer
body: <p>I am doing a similar work to estimate the parameters of Frechet distribution given the guides provided here. The resource provided here (<a href="https://stackoverflow.com/a/28011584/7308826">https://stackoverflow.com/a/28011584/7308826</a>) which was very useful. But, after doing the steps I got the following error, "could not find function "checkparam". Any idea what is going wrong in the code? I am adding the modifeid code below:</p>
  
  #
  
  <h1>Copyright (c) 2010 Marie Laure Delignette-Muller, Christophe Dutang</h1>
  
  <pre><code>    foo_mmedist &lt;- function (data, distr, order, memp, start=NULL, fix.arg=NULL,
                       optim.method="default", lower=-Inf, upper=Inf, custom.optim=NULL, 
                       weights=NULL, silent=TRUE, gradient=NULL, ...) 
  {
    if (!is.character(distr)) 
      stop("distr must be a character string naming a distribution")
    else 
      distname &lt;- distr
  
    if (is.element(distname, c("norm", "lnorm", "pois", "exp", "gamma", "nbinom", "geom",
                               "beta", "unif", "logis")))
      meth &lt;- "closed formula"
    else
      meth &lt;- optim.method
  
    mdistname &lt;- paste("m", distname, sep="")
    ddistname &lt;- paste("d", distname, sep="")
    if(is.null(custom.optim))
      optim.method &lt;- match.arg(optim.method, c("default", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"))
  
    if(!is.null(weights))
    {
      if(any(weights &lt; 0))
        stop("weights should be a vector of integers greater than 0")
      if(!is.allint.w(weights))
        stop("weights should be a vector of (strictly) positive integers")
      if(length(weights) != NROW(data))
        stop("weights should be a vector with a length equal to the observation number")
      warning("weights are not taken into account in the default initial values")
    }
  
    if(meth != "closed formula")
    {
  
      if (!exists(as.character(expression(memp)), mode = "function"))
        stop(paste0("The moment ", mdistname, " function must be defined."))     
      # mdistname contains the good name of the theoretical moment function    
    }
    if (!(is.numeric(data) &amp; length(data)&gt;1)) 
      stop("data must be a numeric vector of length greater than 1")
  
    if(meth == "closed formula")
    {
      n &lt;- length(data)
      if(is.null(weights))
      {
        m &lt;- mean(data)
        v &lt;- (n - 1)/n*var(data)
      }else #weighted version from util-wtdstat.R
      {
        m &lt;- wtd.mean(data, weights=weights)
        v &lt;- wtd.var(data, weights=weights)
      }
  
      if (!is.null(fix.arg))
        stop("argument fix.arg cannot be used when a closed formula is used.")
      # Fitting by matching moments
      if (!(is.vector(data) &amp; is.numeric(data) &amp; length(data)&gt;1))
        stop("data must be a numeric vector of length greater than 1")
      if (distname == "norm") {
        estimate &lt;- c(mean=m, sd=sqrt(v))
        order &lt;- 1:2
      }
      if (distname == "lnorm") {
        if (any(data &lt;= 0)) 
          stop("values must be positive to fit a lognormal distribution")
        sd2 &lt;- log(1+v/m^2)
        estimate &lt;- c(meanlog=log(m) - sd2/2, sdlog=sqrt(sd2))
        order &lt;- 1:2            
      }
  
  
  
      if (distname == "pois") {
        estimate &lt;- c(lambda=m)
        order &lt;- 1          
      }
      if (distname == "exp") {
        estimate &lt;- c(rate=1/m)
        order &lt;- 1          
      }
      if (distname == "gamma" ) {
        shape &lt;- m^2/v
        rate &lt;- m/v
        estimate&lt;-c(shape=shape, rate=rate)
        order &lt;- 1:2            
      }
      if (distname == "nbinom" ) {
        size &lt;- if (v &gt; m) m^2/(v - m)
        else NaN
        estimate&lt;-c(size=size, mu=m)
        order &lt;- 1:2           
      }
      if (distname == "geom" ) {
        prob&lt;-if (m&gt;0) 1/(1+m)
        else NaN
        estimate&lt;-c(prob=prob)
        order &lt;- 1         
      }
      if (distname == "beta" ) {
        if (any(data &lt; 0) | any(data &gt; 1)) 
          stop("values must be in [0-1] to fit a beta distribution")
        aux&lt;-m*(1-m)/v - 1
        shape1 &lt;- m*aux
        shape2 &lt;- (1-m)*aux
        estimate&lt;-c(shape1=shape1, shape2=shape2)
        order &lt;- 1:2            
      }
      if (distname == "unif" ) {
        min1 &lt;- m-sqrt(3*v)
        max1 &lt;- m+sqrt(3*v)
        estimate&lt;-c(min1,max1)
        order &lt;- 1:2            
      }
      if (distname == "logis" ) {
        scale &lt;- sqrt(3*v)/pi
        estimate&lt;-c(location=m, scale=scale)
        order &lt;- 1:2            
      }
      res &lt;- list(estimate=estimate, convergence=0, value=NULL, hessian=NULL,
                  optim.function=NULL, order=order, memp=NULL, counts=NULL,
                  optim.message=NULL)
      opt.meth &lt;- fix.arg.fun &lt;- NULL
  
    }else #an optimimisation has to be done, where fix.arg and start can be a function
    {
      start.arg &lt;- start #to avoid confusion with the start() function of stats pkg (check is done lines 87-100)
      if(is.vector(start.arg)) #backward compatibility
        start.arg &lt;- as.list(start.arg)
  
      # definition of starting/fixed values values
      argmdistname &lt;- names(formals(mdistname))
  
      chfixstt &lt;- checkparam (start.arg=start.arg, fix.arg=fix.arg, argdistname=argmdistname, errtxt=NULL, data10=head(data, 10), distname=distname)
  
      if(!chfixstt$ok)
        stop(chfixstt$txt)
      #unlist starting values as needed in optim()
      if(is.function(chfixstt$start.arg))
        vstart &lt;- unlist(chfixstt$start.arg(data))
      else
        vstart &lt;- unlist(chfixstt$start.arg)
      if(is.function(fix.arg)) #function
      { 
        fix.arg.fun &lt;- fix.arg
        fix.arg &lt;- fix.arg(data)
      }else
        fix.arg.fun &lt;- NULL
      #otherwise fix.arg is a named list or NULL
  
      # end of the definition of starting/fixed values
  
      if(length(vstart) != length(order))
        stop("wrong dimension for the moment order to match")
      if(missing(memp)) 
        stop("the empirical moment function must be defined")
      #backward compatibility when memp is the name of the function and not the function itself
      if(is.character(memp))
        memp &lt;- get0(memp, envir=pos.to.env(1))
  
  
      #check the memp function
      if(!is.function(memp)) 
        stop("the empirical moment must be defined as a function")
      if(is.null(weights))
      {
        txt &lt;- "the empirical moment function must be a two-argument function of 'x', 'order'"
        if(length(formals(memp)) != 2)
          stop(txt)
        if(any(names(formals(memp)) != c("x", "order")))
          stop(txt)
      }else
      {
        txt &lt;- "the empirical moment function must be a three-argument function of 'x', 'order', 'weights'"
        if(length(formals(memp)) != 3)
          stop(txt)
        if(any(names(formals(memp)) != c("x", "order", "weights")))
          stop(txt)
      }
  
  
      ############# MME fit using optim or custom.optim ##########
  
      # definition of the function to minimize : least square (Cramer - von Mises type)
      if(is.null(weights))
      {
        DIFF2 &lt;- function(par, fix.arg, order, obs, mdistnam, memp, weights)
        {
          momtheo &lt;- do.call(mdistnam, c(as.list(order), as.list(par), as.list(fix.arg)) )
          momemp &lt;- as.numeric(memp(obs, order))
          (momemp - momtheo)^2
        }
        fnobj &lt;- function(par, fix.arg, obs, mdistnam, memp, weights)
          sum( sapply(order, function(o) DIFF2(par, fix.arg, o, obs, mdistnam, memp)) )
      }else
      {
        DIFF2 &lt;- function(par, fix.arg, order, obs, mdistnam, memp, weights)
        {
          momtheo &lt;- do.call(mdistnam, c(as.list(order), as.list(par), as.list(fix.arg)) )
          momemp &lt;- as.numeric(memp(obs, order, weights))
          (momemp - momtheo)^2
        }
        fnobj &lt;- function(par, fix.arg, obs, mdistnam, memp, weights)
          sum( sapply(order, function(o) DIFF2(par, fix.arg, o, obs, mdistnam, memp, weights)) )
      }
  
      cens &lt;- FALSE
      if(cens)
        stop("Moment matching estimation for censored data is not yet available.")
  
  
      owarn &lt;- getOption("warn")
      # Try to minimize the stat distance using the base R optim function
      if(is.null(custom.optim))
      {
        hasbound &lt;- any(is.finite(lower) | is.finite(upper))
  
        # Choice of the optimization method  
        if (optim.method == "default")
        {
          opt.meth &lt;- ifelse(length(vstart) &gt; 1, "Nelder-Mead", "BFGS") 
        }else
          opt.meth &lt;- optim.method
  
        if(opt.meth == "BFGS" &amp;&amp; hasbound &amp;&amp; is.null(gradient))
        {
          opt.meth &lt;- "L-BFGS-B"
          txt1 &lt;- "The BFGS method cannot be used with bounds without provided the gradient."
          txt2 &lt;- "The method is changed to L-BFGS-B."
          warning(paste(txt1, txt2))
        }
  
        options(warn=ifelse(silent, -1, 0))
        #select optim or constrOptim
        if(hasbound) #finite bounds are provided
        {
          if(!is.null(gradient))
          {
            opt.fun &lt;- "constrOptim"
          }else #gradient == NULL
          {
            if(opt.meth == "Nelder-Mead")
              opt.fun &lt;- "constrOptim"
            else if(opt.meth %in% c("L-BFGS-B", "Brent"))
              opt.fun &lt;- "optim"
            else
            {
              txt1 &lt;- paste("The method", opt.meth, "cannot be used by constrOptim() nor optim() without gradient and bounds.")
              txt2 &lt;- "Only optimization methods L-BFGS-B, Brent and Nelder-Mead can be used in such case."
              stop(paste(txt1, txt2))
            }
          }
          if(opt.fun == "constrOptim")
          {
            #recycle parameters
            npar &lt;- length(vstart) #as in optim() line 34
            lower &lt;- as.double(rep_len(lower, npar)) #as in optim() line 64
            upper &lt;- as.double(rep_len(upper, npar))
  
            # constraints are : Mat %*% theta &gt;= Bnd, i.e. 
            # +1 * theta[i] &gt;= lower[i]; 
            # -1 * theta[i] &gt;= -upper[i]
  
            #select rows from the identity matrix
            haslow &lt;- is.finite(lower)
            Mat &lt;- diag(npar)[haslow, ]
            #select rows from the opposite of the identity matrix
            hasupp &lt;- is.finite(upper)
            Mat &lt;- rbind(Mat, -diag(npar)[hasupp, ])
            colnames(Mat) &lt;- names(vstart)
            rownames(Mat) &lt;- paste0("constr", 1:NROW(Mat))
  
            #select the bounds
            Bnd &lt;- c(lower[is.finite(lower)], -upper[is.finite(upper)])
            names(Bnd) &lt;- paste0("constr", 1:length(Bnd))
  
            initconstr &lt;- Mat %*% vstart - Bnd
            if(any(initconstr &lt; 0))
              stop("Starting values must be in the feasible region.")
  
            opttryerror &lt;- try(opt &lt;- constrOptim(theta=vstart, f=fnobj, ui=Mat, ci=Bnd, grad=gradient,
                                                  fix.arg=fix.arg, obs=data, mdistnam=mdistname, memp=memp, 
                                                  hessian=!is.null(gradient), method=opt.meth, weights=weights, 
                                                  ...), silent=TRUE)
            if(!inherits(opttryerror, "try-error"))
              if(length(opt$counts) == 1) #appears when the initial point is a solution
                opt$counts &lt;- c(opt$counts, NA)
  
          }else #opt.fun == "optim"
          {
            opttryerror &lt;- try(opt &lt;- optim(par=vstart, fn=fnobj, fix.arg=fix.arg, obs=data, gr=gradient,
                                            mdistnam=mdistname, memp=memp, hessian=TRUE, method=opt.meth, 
                                            lower=lower, upper=upper, weights=weights, ...), silent=TRUE)       
          }
  
        }else #hasbound == FALSE
        {
          opt.fun &lt;- "optim"
          opttryerror &lt;- try(opt &lt;- optim(par=vstart, fn=fnobj, fix.arg=fix.arg, obs=data, gr=gradient,
                                          mdistnam=mdistname, memp=memp, hessian=TRUE, method=opt.meth, lower=lower, 
                                          upper=upper, ...), silent=TRUE)       
        }
        options(warn=owarn)
  
        if (inherits(opttryerror,"try-error"))
        {
          warnings("The function optim encountered an error and stopped.")
          if(getOption("show.error.messages")) print(attr(opttryerror, "condition"))              
          return(list(estimate = rep(NA,length(vstart)), convergence = 100, value = NA, 
                      hessian = NA))
        }
  
        if (opt$convergence&gt;0) {
          warnings("The function optim failed to converge, with the error code ",
                   opt$convergence)
        }
        if(is.null(names(opt$par)))
          names(opt$par) &lt;- names(vstart)
        res &lt;- list(estimate = opt$par, convergence = opt$convergence, value = opt$value, 
                    hessian = opt$hessian, optim.function=opt.fun, order=order, 
                    memp=memp, counts=opt$counts, optim.message=opt$message)  
  
      }else # Try to minimize the stat distance using a user-supplied optim function 
      {
        opt.meth &lt;- NULL
        if (!cens)
        {
          options(warn=ifelse(silent, -1, 0))
          opttryerror &lt;- try(opt &lt;- custom.optim(fn=fnobj, fix.arg=fix.arg, obs=data, mdistnam=mdistname, 
                                                 memp=memp, par=vstart, weights=weights, ...), silent=TRUE)
          options(warn=owarn)
        }else
          stop("Moment matching estimation for censored data is not yet available.")
  
        if (inherits(opttryerror,"try-error"))
        {
          warnings("The customized optimization function encountered an error and stopped.")
          if(getOption("show.error.messages")) print(attr(opttryerror, "condition"))              
          return(list(estimate = rep(NA,length(vstart)), convergence = 100, value = NA, 
                      hessian = NA))
        }
  
        if (opt$convergence&gt;0) {
          warnings("The customized optimization function failed to converge, with the error code ",
                   opt$convergence)
        }
  
        res &lt;- list(estimate = opt$par, convergence = opt$convergence, value = opt$value, 
                    hessian = opt$hessian, optim.function=custom.optim, order=order, 
                    memp=memp, counts=opt$counts, optim.message=opt$message)  
  
      }   
  
    }
  
    if(is.null(weights))
    {  
      loglik &lt;- function(par, fix.arg, obs, ddistnam) 
        sum(log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)) ) ) )
    }else
    {
      loglik &lt;- function(par, fix.arg, obs, ddistnam) 
        sum(weights * log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)) ) ) )
    }
    if(exists(ddistname))
      loglik &lt;- loglik(res$estimate, fix.arg, data, ddistname)
    else
      loglik &lt;- NULL
  
    res &lt;- c(res, fix.arg=fix.arg, loglik=loglik, method=meth, optim.method=opt.meth, 
             fix.arg.fun=fix.arg.fun, list(weights = weights))
  
    return(res)
  }
  
  ## old function with previous name 
  momdist&lt;-function (data, distr) 
  {
    stop("the name \"momdist\" for matching moments function is NO MORE used and is replaced by \"mmedist\".")
  }
  </code></pre>