I want to create a QQ & PP plot using L-Moments to fit the GEV but get this error: unable to perform assignment b/c the left & right sides have a different number of elements

5 views (last 30 days)
Here is my code:
% Step A: Read data and extract columns
data = readmatrix('Wabash_River_MaxFlow.xls');
Error using readmatrix
Unable to find or open 'Wabash_River_MaxFlow.xls'. Check the path and filename or file permissions.
wateryear = data(:, 1);
streamflow = data(:, 2);
% Step B: Compute L-Moments
n = length(data);
b0 = mean(data);
run_sum = 0;
for i = 1:(n-1)
run_sum = run_sum + (i+1-1)*data(i+1);
end
b1 = run_sum/(n*(n-1));
run_sum = 0;
for i = 2:(n-1)
run_sum = run_sum + (i+1-1)*(i+1-2)*data(i+1);
end
b2 = run_sum/(n*(n-1)*(n-2));
run_sum = 0;
for i = 3:(n-1)
run_sum = run_sum + (i+1-1)*(i+1-2)*(i+1-3)*data(i+1);
end
b3 = run_sum/(n*(n-1)*(n-2)*(n-3));
run_sum = 0;
for i = 4:(n-1)
run_sum = run_sum + (i+1-1)*(i+1-2)*(i+1-3)*(i+1-4)*data(i+1);
end
b4 = run_sum/(n*(n-1)*(n-2)*(n-3)*(n-4));
l1 = b0;
l2 = 2*b1-b0;
l3 = 6*b2-6*b1+b0;
l4 = 20*b3-30*b2+12*b1-b0;
L_CV = l2/l1;
L_SK = l3/l2;
L_KU = l4/l2;
% Step C: Generate QQ plot
X_plot = zeros(1, 80);
pp_plot = linspace(1/81, 1, 80); % Update: generate pp_plot directly
for i = 1:80
X_plot(i) = psi + alpha*(1-(-1*log(pp_plot(i)))^k)/k;
end
figure;
scatter(pp_plot, streamflow);
hold on;
plot(pp_plot, X_plot, '-r');
xlabel('Non-Exceedance Prob.');
ylabel('Observations');
title('QQ Plot');
legend('Observations', 'GEV Fit');
% Step D: Generate PP plot
pp_data = linspace(1/(n+1), n/(n+1), n); % Update: generate pp_data directly
pp_theory = psi + alpha*(1-(-1*log(pp_data))^k)/k;
figure;
plot(linspace(0, 1, 100), linspace(0, 1, 100), '-r');
hold on;
scatter(pp_theory, pp_data);
xlabel('Theoretical Non-Exceedance Prob.');
ylabel('Fitted Non-Exceedance Prob.');
title('PP Plot');
legend('1:1 Line', 'GEV Fit');
Here is the error I am getting:
Unable to perform assignment because the left and right sides have a different
number of elements.
Error in CEE204_HW2_Problem1 (line 47)
X_plot(i) = psi + alpha*(1-(-1*log(pp_plot(i)))^k)/k;
Here is what my data looks like:

Answers (1)

Jeff Miller
Jeff Miller on 6 Apr 2023
Assign the problematic result to a new variable so that you can see what elements it has, like this:
temp_var = psi + alpha*(1-(-1*log(pp_plot(i)))^k)/k;
Probably you will see that temp_var does not have the number of elements you expected. For example, maybe instead you need:
X_plot(i) = psi + alpha*(1-(-1*log(pp_plot(i))).^k)/k;

Products


Release

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!