OK, w/o a lot of effort a first pass of the basic idea might be sotoo...
rain=importdata('rain.dat');
rain=rain.data(:,4);
daymo=diff(datenum(1980,1:12*6+1,1,0,0,0));
r=rain./daymo;
r=reshape(r,12,[]);
format bank
[r mean(r,2) var(r,[],2)]
figure, plot(r), xlim([1 12])
OK, I learned the exponential with parameters eta and lambda where eta is the Matlab a and lambda --> 1/b
For gamma, the mean is eta/lambda and variance is mean/lambda or eta/lambda*2. Solving for those by month leads us by some simple algebra to
lamb=mean(r,2)./var(r,[],2);
eta=mean(r,2).*lamb;
x = gaminv((0.005:0.01:0.995),eta(1),1/lamb(1));
y = gampdf(x,eta(1),1/lamb(1));
figure,plot(x,y)
rgam=gamrnd(eta(1),1./lamb(1),31,1);
OK, here's some results...
>> [min(rgam) max(rgam) mean(rgam) sum(rgam)]
ans =
0.30 1.92 0.85 26.45
>>
That's at least in the realm of reasonableness altho you'll likely need to do something to capture the real variability a little more on an annual basis...
The alternative is, as said, to decide what the extremes in your observed values might represent in terms of the realistic upper limits and use the cdf percentage points and back solve for parameters that match that and, say, the mean and see what that leads to.
You've got any number of possibilities given the lack of actual data structure...
