# Calculating confidence interval of return period

I have a time series of annual rainfall data and I wish to calculate the return period curve for this data. To do this, I have used Matlab to fit a gamma distribution, and then used the Lilliefors test to test the distribution fit. If I cannot reject the null hypothesis at the 5% level, I will use the fitted gamma distribution to calculate the return period. I have done this fine, but what I need to do is calculate the confidence interval for this return period curve.

How do I go about doing this? I am not a statistician, so I am struggling to understand a lot of the information I have found online about how to do this. If anyone can explain it to me quite simply, that would be great! I appreciate there may be inbuilt ways of doing this in Matlab, but I need to understand the process behind it so I can write it up in my thesis.

This image is an example of what I want to produce.

This is done using the extRemes package for R, but I want to recreate this myself as the package doesn't use the gamma distribution to do the fitting. I think this figures uses the GEV distribution. I have currently reproduced the black line (the return period curve calculated after fitting the gamma distribution) from my data, but I don't know how to do the blue lines (the confidence interval).

• Can you explain further? The graph you show looks like a regression fit, maybe polynomial. Your text describes a distribution fit. I don't understand. – Tom Lane Nov 26 '14 at 15:52
• This graph is produced by fitting a GEV distribution to the underlying data, and then calculating the return period from the GEV CDF, where return period = 1/(1-CDF). After doing some more research, the above confidence intervals are calculated using the profile-likelihood method, which is simple for the GEV distribution (it is described in Stuart Coles (2001)). It is next to impossible to use this method for gamma, so I am now doing some sort of bootstrapping - still working on it. – emmalgale Nov 27 '14 at 11:10

x = gamrnd(100,1.3,100,1);