Гамма-распределение по данным измерений

[cgem]

close all; clear all;

measurements = dlmread('measurements.log',',');

% Compute histogram
upperb = ceil(max(measurements));
lowerb = floor(min(measurements));
x = linspace(lowerb, upperb, max(upperb + 1, 100))';
n = histc(measurements, x);

% Compute theoretical Gamma PDF
ex = mean(measurements);
dx = var(measurements);
theta = dx / ex;
k = ex / theta;
gp = gampdf(x, theta, k);

h = figure;
hold on;
stairs(x, n, 'k-', 'LineWidth', 2);
plot(x, length(measurements) * gp, 'r-', 'LineWidth', 2);
legend('real distribution', 'theoretical distribution')
axis([х(1) x(end) 0 1.1 * max(n)])
xlabel('Time, sec')
ylabel('Count, number')
print(h,'-dpng','-color','distribution.png')