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
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Here is my code:
% Step A: Read data and extract columns
data = readmatrix('Wabash_River_MaxFlow.xls');
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:
0 commentaires
Réponses (1)
Jeff Miller
le 6 Avr 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;
0 commentaires
Voir également
Catégories
En savoir plus sur Generalized Extreme Value Distribution dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!