-
Notifications
You must be signed in to change notification settings - Fork 801
/
gammaRainfallDemo.m
52 lines (43 loc) · 1.15 KB
/
gammaRainfallDemo.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
%% Fit a Gamma distribution to the rainfall data by MLE and moment matching
% Rice (1995) p383
%%
% This file is from pmtk3.googlecode.com
function gammaRainfallDemo()
X = loadData('rainfall');
X = X'; X = X(:); % concatenate across rows, not columns
%X = X(1:end-5); % removing trailing 0s
[a(1), b(1)] = gamMOM(X);
[a(2), b(2)] = gamMLE(X);
[v, binc] = hist(X);
h = binc(2)-binc(1);
N = length(X);
areaH = h*N;
figure;bar(binc, v/areaH);hold on
xs = linspace(0.05, binc(end));
linestyles = {'k:', 'r-'};
for i=2
ps = exp(gammaLogprob(struct('a', a(i), 'b', b(i)), xs));
h(i) = plot(xs, ps, linestyles{i}, 'linewidth', 3);
end
printPmtkFigure('rainfallDemo');
figure;bar(binc, v/areaH);hold on
for i=1:2
ps = exp(gammaLogprob(struct('a', a(i), 'b', b(i)), xs));
h(i) = plot(xs, ps, linestyles{i}, 'linewidth', 3);
end
legend(h, 'MoM', 'MLE')
printPmtkFigure('rainfallDemoMOMMLE');
end
function [a,b] = gamMLE(X)
%% MLE for Gamma a = shape, b= rate (not scale)
[a, s] = gamma_fit(X);
b = 1/s;
end
function [a,b] = gamMOM(X)
%% method of moments estimate for Gamma
% a = shape, b= rate
xbar = mean(X);
s2hat = var(X);
a = xbar^2/s2hat;
b = xbar/s2hat;
end