0
$\begingroup$

I am fitting a GEV distribution to some rainfall data, but the software I am using (Matlab and Easyfit) are giving a distribution which includes negative numbers (i.e. negative rainfall). Is there a way of forcing this distribution to be non-negative in Matlab or Easyfit?

I tend to get this issue when I have lots of zeros or small values in the data, such as in the CDF below.

GEV CDF

$\endgroup$
  • 3
    $\begingroup$ The GEV is in the first instance a distribution for extremes, not all values. If your zeros really are extremes, then at least some of your data are for dry seasons and/or extremely arid places. Perhaps you are mixing data in some way. So, as always, the question of whether the data are suitable for the distribution also arises. $\endgroup$ – Nick Cox Apr 22 '15 at 2:23
  • $\begingroup$ Out of curiosity: how there can be negative rainfall..? $\endgroup$ – Tim Jul 1 '15 at 17:02
  • $\begingroup$ @Tim, there can't be negative rainfall, that is my problem. My data doesn't have negative rainfall (which you can see on the empirical CDF above). When I fit the distribution to my data, it doesn't force the CDF to be positive, hence my question. $\endgroup$ – emmalgale Jul 2 '15 at 10:15
3
$\begingroup$

First question would be if your problem maps appropriately to the GEV? Are your rainfall measurements actually maxima of blocks (months, years, etc)? If not, there's no reason that the fit should be meaningful.

At what level are you calculating your maxima? For example if you're calculating the maximum of daily readings each week, can you instead aggregate at the monthly level, making zero (and other low values) less likely, and thus pulling the distribution away from negative values?

You should also consider using all of your data and using a POT (Points or Peaks Over Threshold, which uses the Generalized Pareto Distribution) approach, which would avoid zeros all together. This approach uses more data and can take into account the fact that maxima may occur in clusters. (R's evd package provides functions to do this.)

Second question could be, if your measurements contain 0, have you considered using a Hurdle Model? (That is a two-stage model where the first stage predicts zeros, and the second stage applies to non-zeros.) This may not actually be appropriate or helpful, but it might. Speaking of which, are all of your 0's actually 0's, or might some of them actually be NA's -- there was no measurement during that time period?

Third question would be: why do you care about negative values? In general, with a GEV you are looking at things like the N-year return, which is the $1 - {1 \over N}$ quantile, so you'd be looking at the 95 percentile for a 20-year return, the 99 percentile for the 100-year return, etc, so does the 1 or 2 percentile really matter? (I know this is a slippery slope that can do things like "justify" using a Normal approximation for values that cannot be negative, but if you're actually doing an extreme value analysis perhaps it doesn't matter as much what happens in the left tail of the extreme values.)

Is it at all acceptable to consider values less than zero to be zero? This is ad hoc, but it could turn out that if you have excess zeros it actually works. Not a fan of this idea, just throwing it out there.

Last, if you're able to force the shape parameter to be > 0, or equivalently if you can specify a Frechet distribution, the support of the Frechet is positive. This would exclude 0, which could be a problem, and I'm not sure that this would make sense in terms of your subject matter at any rate, but mathematically it could help. Again, just throwing the idea out there.

$\endgroup$
0
$\begingroup$

The GEV has a lower limit, which is $\mu-\frac{\sigma}{k}$, in Matlab notation. What you would need is to fix the parameters in the subspace $\mu=\frac{\sigma}{k}$, which I think is not possible in an easy way in Matlab. What I would do (excluding to avoid Matlab) is to use the function negloglik (see here): it gives you the log-likelihood of the distribution you want for the data you have, and then making a simple loop in which you optimize this function, where the parameter $\mu$ follows the constraint you are giving, and you propose $\sigma$ and $k$ freely.

slightly off topic: I see that you are fitting precipitation data. Doesn't it has a threshold anyways?

$\endgroup$
  • $\begingroup$ If it makes it any easier, I don't have to use Matlab to do this. Is it easier to do it in another program? My data has a lower limit of zero (obviously as you can't get negative rainfall), but if I fit the GEV in Matlab, it gives the above distribution which doesn't have a lower bound of zero. $\endgroup$ – emmalgale Dec 15 '14 at 14:31
0
$\begingroup$

In MATLAB, you can use the MLE function to fix the subspace as chuse suggested. Here's an example. First an unconstrained fit.

>> rng(0)
>> x = gevrnd(.3,.3,1,1000,1); ecdf(x)
>> p = gevfit(x)
p =
    0.2947    0.2785    0.9901
>> gevinv(0,p(1),p(2),p(3))
ans =
    0.0450

Using the unconstrained fit, the lower bound is positive, which I know is not a problem for you but that's the way it turned out. Let's make it zero. Just optimize over k and sigma, setting mu=sigma/k.

>> p2 = mle(x,'pdf',@(x,k,sigma)gevpdf(x,k,sigma,sigma/k),'start',[1 1])
p2 =
    0.2818    0.2795
>> gevinv(0,p2(1),p2(2),p2(2)/p2(1))
ans =
     0

However, if the pdf is asymptotic to the x axis at the lower bound, I suspect the presence of exact zeros in your data is going to cause problems. Even without that, it's sometimes better to estimate log(k) and log(sigma) so that MLE won't try values outside the valid range of the parameters.

$\endgroup$
  • $\begingroup$ I have tested the above code on some of my data (which doesn't contain any zeros in it, but the GEV fit does produce negative rainfall in the lower tail), and I get an error: "Error using mlecustom>llf_pdfcdf (line 440). The PDF function returned negative or zero values". Any ideas? $\endgroup$ – emmalgale Dec 19 '14 at 14:44
  • $\begingroup$ Okay, I cooked up an example that had this issue, but it went away if I estimated on the log scale (sorry I don't know how to format code well in a comment): p2 = mle(x,'pdf',@(x,logk,logsigma)gevpdf(x,exp(logk),exp(logsigma),exp(logsigma-logk)),'start',[0 0]); p2 = exp(p2); gevinv(0,p2(1),p2(2),p2(2)/p2(1)) $\endgroup$ – Tom Lane Dec 20 '14 at 15:59

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.