<p>This is going to be a long story, possibly more suited to <a href="https://stats.stackexchange.com/">https://stats.stackexchange.com/</a>.</p>

<p>====== Part 1 -- The problem ======</p>

<p>This is the sequence generating the error:</p>

<pre><code>library(fExtremes)
samp &lt;- read.csv("optimdata.csv")[ ,2]
## does not converge
para &lt;- gevFit(samp, type = "mle")
</code></pre>

<p>We are facing the typical cause of lack-of-convergence when using <code>optim()</code> and friends: inadequate starting values for the optimisation.</p>

<p>To see what goes wrong, let us use the PWM estimator (<a href="http://arxiv.org/abs/1310.3222" rel="nofollow noreferrer">http://arxiv.org/abs/1310.3222</a>); this consists of an analytical formula, hence it does not incur into convergence problems, since it makes no use of <code>optim()</code>:</p>

<pre><code>para &lt;- gevFit(samp, type = "pwm")
fitpwm&lt;- attr(para, "fit")
fitpwm$par.ests
</code></pre>

<p>The estimated tail parameter <code>xi</code> is negative, corresponding to a bounded upper tail; in fact the fitted distribution displays even more "upper tail boundedness" than the sample data, as you can see from the "leveling off" of the quantile-quantile graph at the right:</p>

<pre><code>qqgevplot &lt;- function(samp, params){
  probs &lt;- seq(0.1,0.99,by=0.01)
  qqempir &lt;- quantile(samp, probs)
  qqtheor &lt;-  qgev(probs, xi=params["xi"], mu=params["mu"], beta=params["beta"])
  rang &lt;- range(qqempir,qqtheor)
  plot(qqempir, qqtheor, xlim=rang, ylim=rang,
     xlab="empirical", ylab="theoretical",
     main="Quantile-quantile plot")
  abline(a=0,b=1, col=2)
}
qqgevplot(samp, fitpwm$par.ests)
</code></pre>

<p>For <code>xi&lt;0.5</code> the MLE estimator is not regular (<a href="http://arxiv.org/abs/1301.5611" rel="nofollow noreferrer">http://arxiv.org/abs/1301.5611</a>): the value of  -0.46 estimated by PWM for <code>xi</code> is very close to that. Now the PWM estimates are used internally by <code>gevFit()</code> as starting values for <code>optim()</code>: you can see this if you print out the code for the function <code>gevFit()</code>:</p>

<pre><code>print(gevFit)
print(.gevFit)
print(.gevmleFit)
</code></pre>

<p>The starting value for optim is <code>theta</code>, obtained by PWM. For the specific data at hand, this starting value is not adequate, in that it leads to non-convergence of  <code>optim()</code>.</p>

<p>====== Part 2 -- solutions? ======</p>

<p>Solution 1 is to use <code>para &lt;- gevFit(samp, type = "pwm")</code> as above. If you'd like to use ML, then you have to specify good starting values for <code>optim()</code>. Unfortunately, the <code>fExtremes</code> package does not make it easy to do so. You can then re-define your own version of <code>.gevmleFit</code> to include those, e.g.</p>

<pre><code>.gevmleFit &lt;- function (data, block = NA, start.param, ...) 
{
  data = as.numeric(data)
  n = length(data)
  if(missing(start.param)){
    theta = .gevpwmFit(data)$par.ests
  }else{
    theta = start.param
  }
  fit = optim(theta, .gevLLH, hessian = TRUE, ..., tmp = data)
  if (fit$convergence) 
    warning("optimization may not have succeeded")
  par.ests = fit$par
  varcov = solve(fit$hessian)
  par.ses = sqrt(diag(varcov))
  ans = list(n = n, data = data, par.ests = par.ests, par.ses = par.ses, 
             varcov = varcov, converged = fit$convergence, nllh.final = fit$value)
  class(ans) = "gev"
  ans
}
## diverges, just as above
.gevmleFit(samp)
## diverges, just as above
startp &lt;- fitpwm$par.ests
.gevmleFit(samp, start.param=startp)
## converges
startp &lt;- structure(c(-0.1, 1, 1), names=names(fitpwm$par.ests))
.gevmleFit(samp, start.param=startp)$par.ests
</code></pre>

<p>Now check this out: the <code>beta</code> estimated by PWM is 0.1245; by changing this to a tiny amount, the MLE is made to converge:</p>

<pre><code>startp &lt;- fitpwm$par.ests
startp["beta"]
startp["beta"] &lt;- 0.13
.gevmleFit(samp, start.param=startp)$par.ests
</code></pre>

<p>This hopefully clearly illustrates that to blindly <code>optim()</code>ise works until it doesn't and might then turn into a quite delicate endeavour indeed. For this reason, it might be useful to leave this reply here, rather than to migrate to CrossValidated.</p>
R optim(){fExtremes} gets 0 hessian matrix
<p>I am using R {fExtremes} to find best parameters of GEV distribution for my data (a vector). but get the following error message</p>

<blockquote>
  <p>Error in solve.default(fit$hessian) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0</p>
</blockquote>

<p>I traced back to fit$hessian, found my hessian matrix is a sigular matrix, all of the elements are 0s. The source code (<a href="https://github.com/cran/fExtremes/blob/master/R/GevFit.R" rel="nofollow">https://github.com/cran/fExtremes/blob/master/R/GevFit.R</a>) of gevFit()  shows fit$hessian is calculated by optim(). The output parameters are exactly the same value as the initial parameters. I am wondering what could be the problems of my data that cause this problem? I copied my code here </p>

<pre><code>&gt; min(sample);
[1] 5.240909

&gt; max(sample)
[1] 175.8677

&gt; length(sample)
[1] 6789

&gt; mean(sample)
[1] 78.04107

&gt;para&lt;-gevFit(sample, type = "mle")
Error in solve.default(fit$hessian) : 
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

fit = optim(theta, .gumLLH, hessian = TRUE, ..., tmp = data)
&gt; fit

   $par

xi   -0.3129225
mu   72.5542497 
beta  16.4450897 

$value
[1] 1e+06

$counts
function gradient 
       4       NA 

$convergence
[1] 0

$message
NULL



$hessian

     xi  mu beta

xi    0    0     0

mu    0    0      0

beta  0     0      0
</code></pre>

<p>I updated my dataset on google docs: 
<a href="https://docs.google.com/spreadsheets/d/1IRRpjmdrrJPhNmfiLism_P0efV_Ot4HlEsa6kwMnljc/edit?usp=sharing" rel="nofollow">https://docs.google.com/spreadsheets/d/1IRRpjmdrrJPhNmfiLism_P0efV_Ot4HlEsa6kwMnljc/edit?usp=sharing</a></p>