How to apply loop on following case?

1 view (last 30 days)
Ahmed
Ahmed on 1 Mar 2024
Commented: Ahmed on 2 Mar 2024
This is my code where I computed Prediction interval coverage probability for IP_OPT and now want to compute for IS_OPT and RH_OPT (line 2 and 3).
One way to write a separate code for IS_OPT and RH_OPT as I did for IP_OPT which looks not a good way to make code. How can make a loop here for three all three IS_OPT, IS_OPT and RH_OPT to get three separate figures as shown below.
IP_p = prctile(IP_OPT,[0.5 99.5],2);
IS_p = prctile(IS_OPT,[2.5 97.5],2);
RH_p = prctile(RH_OPT,[2.5 97.5],2);
% Assuming have a dataset with n_total_points
n_total_points = size(IP_OPT, 1);
% Assuming you are using a certain percentage for training (e.g., 100%)
train_percentage = 1;
n_train_points = round(train_percentage * n_total_points);
d_obs1 = ip;
q = IP_OPT;
[P50,P1,P99,P] = CI_prob(q);
PICP_train = [];
PICP_pred = [];
for i = 1:10
ind1 = 21 + (i - 1) * (-1);
upper = P(:, ind1);
lower = P(:, i);
cp_indx = (upper>=d_obs1)&(d_obs1>=lower);
% Calculate the coverage probability for the training data
CP_train = sum(cp_indx(1:n_train_points)) / n_train_points;
PICP_train = [PICP_train CP_train];
end
ref_line = 0:0.1:1;
CP_theo = [0.99 0.9:-0.1:0.1];
PICP_linear_train_without = PICP_train;
h2 = figure();
plot(CP_theo,PICP_linear_train_without,'d');
hold on;
plot(ref_line,ref_line,'-.r','linewidth',1.0);
grid on;
xlabel('Theoretical Coverage Probability');
ylabel('Actual Coverage Probability (Training)');
title('Prediction Interval Coverage Probability');
legend('Actual Training Data', 'Ideal Linear Behavior');
%% Function
function [P50,P1,P99,P] = CI_prob(q)
Y = prctile(q,[0.5 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 99.5],2);
P50 = Y(:,11);
P1 = Y(:,end);
P99 = Y(:,1);
P = Y;
% ind = [1,11,21];
% P(:,ind) = [];
end
  11 Comments
VBBV
VBBV on 2 Mar 2024
Where do you use h2 in your code ? elswhere in other code snippets ?
h2 = figure()
Ahmed
Ahmed on 2 Mar 2024
@VBBV I think the problem is here
d_obs1 = ip;
and
cp_indx = (upper>=d_obs1)&(d_obs1>=lower);
loop should also work on it e.g., for ip is and rho e.g.,
d_obs1 = [ip is rho] for every q, d_obs1 should be corresponding one like for IP_OPT is ip, for IS_OPT is is, and for RH_OPT is rho.

Sign in to comment.

Accepted Answer

VBBV
VBBV on 2 Mar 2024
Try using the cell array for the testdata you gave
data = load('datatest.mat');
X_p = {data.IP_OPT, data.IS_OPT, data.RH_OPT} % make a cell array
X_p = 1x3 cell array
{101x500 double} {101x500 double} {101x500 double}
D_p = [0.5 99.5;2.5 97.5;2.5 97.5]; % make a vector
for k = 1:3 % just to be sure of elements are string scalars
F_p = prctile(X_p{k},[D_p(k,:)],2); % use a cell array
n_total_points = size(F_p, 1);
% Assuming you are using a certain percentage for training (e.g., 100%)
train_percentage = 1;
n_train_points = round(train_percentage * n_total_points);
d_obs1 = ip; % what is ip here thats not define anywhere
q = X_p{k}; % use a cell array
[P50,P1,P99,P] = CI_prob(q);
%
PICP_train = [];
PICP_pred = [];
for i = 1:10
ind1 = 21 + (i - 1) * (-1);
upper = P(:, ind1);
lower = P(:, i);
cp_indx = (upper>=d_obs1)&(d_obs1>=lower);
% Calculate the coverage probability for the training data
CP_train = sum(cp_indx(1:n_train_points)) / n_train_points;
PICP_train = [PICP_train CP_train];
end
% the below lines need to be modified similar to function calls
ref_line = 0:0.1:1;
CP_theo = [0.99 0.9:-0.1:0.1];
PICP_linear_train_without = PICP_train;
h2 = figure(k);
plot(CP_theo,PICP_linear_train_without,'d');
hold on;
plot(ref_line,ref_line,'-.r','linewidth',1.0);
grid on;
% update figure labels and legends similarly for each figure if neeeded
xlabel('Theoretical Coverage Probability');
ylabel('Actual Coverage Probability (Training)');
title('Prediction Interval Coverage Probability');
legend('Actual Training Data', 'Ideal Linear Behavior');
end
Unrecognized function or variable 'ip'.
%% Function
function [P50,P1,P99,P] = CI_prob(q)
Y = prctile(q,[0.5 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 99.5],2);
P50 = Y(:,11);
P1 = Y(:,end);
P99 = Y(:,1);
P = Y;
% ind = [1,11,21];
% P(:,ind) = [];
end
  3 Comments
VBBV
VBBV on 2 Mar 2024
data = load('datatest.mat');
iprho = load('ipisrho.mat')
iprho = struct with fields:
ip: [101×1 double] is: [101×1 double] rho: [101×1 double]
X_p = {data.IP_OPT, data.IS_OPT, data.RH_OPT} % make a cell array
X_p = 1×3 cell array
{101×500 double} {101×500 double} {101×500 double}
ipisrho = {iprho.ip iprho.is iprho.rho};
D_p = [0.5 99.5;2.5 97.5;2.5 97.5]; % make a vector
for k = 1:3 % just to be sure of elements are string scalars
F_p = prctile(X_p{k},[D_p(k,:)],2); % use a cell array
n_total_points = size(F_p, 1);
% Assuming you are using a certain percentage for training (e.g., 100%)
train_percentage = 1;
n_train_points = round(train_percentage * n_total_points);
d_obs1 = ipisrho{k}; % use cell array like before
q = X_p{k}; % use a cell array
[P50,P1,P99,P] = CI_prob(q);
%
PICP_train = [];
PICP_pred = [];
for i = 1:10
ind1 = 21 + (i - 1) * (-1);
upper = P(:, ind1);
lower = P(:, i);
cp_indx = (upper>=d_obs1)&(d_obs1>=lower);
% Calculate the coverage probability for the training data
CP_train = sum(cp_indx(1:n_train_points)) / n_train_points;
PICP_train = [PICP_train CP_train];
end
% the below lines need to be modified similar to function calls
ref_line = 0:0.1:1;
CP_theo = [0.99 0.9:-0.1:0.1];
PICP_linear_train_without = PICP_train;
h2 = figure(k);
plot(CP_theo,PICP_linear_train_without,'d');
hold on;
plot(ref_line,ref_line,'-.r','linewidth',1.0);
grid on;
% update figure labels and legends similarly for each figure if neeeded
xlabel('Theoretical Coverage Probability');
ylabel('Actual Coverage Probability (Training)');
title('Prediction Interval Coverage Probability');
legend('Actual Training Data', 'Ideal Linear Behavior');
end
%% Function
function [P50,P1,P99,P] = CI_prob(q)
Y = prctile(q,[0.5 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 99.5],2);
P50 = Y(:,11);
P1 = Y(:,end);
P99 = Y(:,1);
P = Y;
% ind = [1,11,21];
% P(:,ind) = [];
end
Ahmed
Ahmed on 2 Mar 2024
@VBBV Thank you dear, its fine now, one last thing is if I want plot all the three figues in loop as a subplot (e.g., subplot(131) subplot(132) and so on. How I can amend the code?

Sign in to comment.

More Answers (1)

Torsten
Torsten on 1 Mar 2024
Edited: Torsten on 1 Mar 2024
Make a function of the part of the code that you have to run through three times and call this function for IP_OPT, IS_OPT and RH_OPT.
Consider to make the plotting in the calling script part.
  2 Comments
Ahmed
Ahmed on 1 Mar 2024
did not understnd you answer.
Torsten
Torsten on 1 Mar 2024
Edited: Torsten on 1 Mar 2024
Small example:
You want to compute x^2 for x = 1,2,3:
x = [1 2 3];
for i = 1:numel(x)
f(i) = square(x(i));
end
f
f = 1×3
1 4 9
function f = square(x)
f = x^2;
end
Now imagine x = IP_OPT, IS_OPT and RH_OPT.
Put the necessary commands in a function (like the square function above) such that it returns PICP_train, PISP_train and PRHP_train when called with IP_OPT, IS_OPT and RH_OPT.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!