Testing goodness of fit for a Zipf distribution (in Matlab)

Question

I have several ranking distributions and would, for each one, like to fit a [Zipf distribution][1], and estimate the goodness of fit relative to some standard benchmark.

With the Matlab code below, I tried to do a sanity check and see if a "textbook" Zipf rank distribution passes the statistical test. Clearly something is wrong, as it does not. If that doesn't, nothing will!

Using the Kolmogorov-Smirnoff test, or the Anderson-Darling test with a custom-built (non-normal) distribution in place of the chi-squared test does not change this.

% Define some empirical frequency distribution
x = 1:10;
freq = randn(1,10); % textbook zipf!

% Define the Zipf distribution
alpha = 1.5; % Shape parameter, 1.5 is apparently a good all-round value to start with
N = sum(freq); % Total number of observations
k = 1:length(x); % Rank of each observation
zipf_dist = N ./ (k.^alpha); % Compute the Zipf distribution

% Plot our empirical frequency distribution alongside the Zipf distribution
figure;
bar(x, freq); % or freq\N
hold on;
plot(x, zipf_dist, 'r--');
xlabel('Rank');
ylabel('Frequency');
legend('Observed', 'Zipf');

% Compute the goodness of fit using the chi-squared test
expected_freq = zipf_dist .* N;
chi_squared = sum((freq - expected_freq).^2 ./ expected_freq);
dof = length(freq) - 1;
p_value = 1 - chi2cdf(chi_squared, dof);

% Display the results
fprintf('Chi-squared statistic = %.4f\n', chi_squared);
fprintf('p-value = %.4f\n', p_value);
if p_value < 0.05
    fprintf('Conclusion: The data is not from a Zipf distribution.\n');
else
    fprintf('Conclusion: The data is from a Zipf distribution.\n');
end

Answer 1

This is a way of doing it with the Kolmogorov-Smirnoff test. To plug in your actual date you'd have to change y_obs. I think this is what you wanted, but please let me know if you have any questions.

% Define the Zipf distribution
N = 100;
alpha = 1;  % Shape parameter, 1.5 is apparently a good all-round value to start with
x = 1:N-1;
y_exp = 1./(x.^alpha); % A correct Zipf distribution

% Create a noisy signal
rng(42)
noise = normrnd(0, 0.01, [1, N-1]);
y_obs = y_exp + noise;

% Plot our empirical frequency distribution alongside the Zipf distribution
plot(x, y_obs, 'r--')
hold on
plot(x, y_exp, 'b--')
xlabel('Rank')
ylabel('Frequency')
legend('Observed', 'Zipf')
hold off

% Compute the goodness of fit using the Kolmogorov-Smirnov test
[H, p_value, D] = kstest2(y_obs, y_exp);

% Display the results
fprintf('Kolmogorov-Smirnov statistic = %.4f\n', D)
fprintf('p-value = %.4f\n', p_value)
if p_value < 0.05
    fprintf('Conclusion: The data is not from a Zipf distribution.\n')
else
    fprintf('Conclusion: The data is from a Zipf distribution.\n')
end

[EDIT]

You could also try to fit an exponential function onto your data. In MATLAB, you can use the anova1 function to perform an analysis of variance (ANOVA) test on the residuals of the exponential regression, which can give you a p-value that indicates whether the model is a good fit for the data.

% Define the exponential function to fit
function y = func(x, p)
    y = p(1) .* x .^ p(2);
end

% Define the dataset
xdata = [1, 2, 3, 4, 5];
ydata = [2.7, 7.5, 16.6, 30.2, 48.6];

% Fit the exponential function to the dataset using lsqcurvefit
p0 = [1, 1]; % Initial parameter values
popt = lsqcurvefit(@(p, xdata) func(xdata, p), p0, xdata, ydata);

% Extract the fitting parameters
a = popt(1);
b = popt(2);

% Calculate the residuals and perform an ANOVA test
yfit = func(xdata, popt);
resid = ydata - yfit;
[pval, tbl, stats] = anova1(resid);

% Print the fitting parameters and p-value
fprintf('a = %f\n', a);
fprintf('b = %f\n', b);
fprintf('p-value = %f\n', pval);

% Check significance against alpha = 0.05
if pval < 0.05
    disp('The exponential regression model is not a good fit for the data.');
else
    disp('The exponential regression model is a good fit for the data.');
end

% Plot the original data and the fitted curve
plot(xdata, ydata, 'ko', 'DisplayName', 'Original data')
hold on
plot(xdata, func(xdata, popt), 'r-', 'DisplayName', 'Fitted curve')
xlabel('x')
ylabel('y')
legend('Location', 'best')
hold off