How to apply loop on following case?

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

[VBBV]

if CI_prob is your core function that is called several times, then you can do something like this using a loop

X_p = [IP_p, IS_p, RH_p];
D_p = [0.5 99.5;2.5 97.5;2.5 97.5];

for k = 1:numel(X_p)
    F_p = prctile(X_p(k,[D_p(k,:)],2); 
    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;
    q = X_p(k);
    [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]

@VBBV it is giving following error

Index in position 2 is invalid. Array indices must be positive integers or logical values.

Error in uncertainity_part_2 (line 42)

F_p = prctile(X_p(k,[D_p(k,:)],2));

[VBBV]

F_p = prctile(X_p(k),[D_p(k,:)],2);   % sorry a typo error 
                  %^

Please check now.

[VBBV]

Also, what are [IP_p, IS_p, RH_p]; ?

Are they numeric arrays ? or strings ? if they are strings, then use a cell array

X_p = {'IP_p', 'IS_p', 'RH_p'};

then code will become something like

X_p = [IP_p, IS_p, RH_p];
D_p = [0.5 99.5;2.5 97.5;2.5 97.5];

for k = 1:3  % just to be sure of elements are string scalars
    F_p = prctile(X_p(k),[D_p(k,:)],2);   %corrected typo
    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;
    q = X_p(k);
    [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]

Also, what are [IP_p, IS_p, RH_p]; These represents the confidence interval (upper and lower limits each size is 101 1)

[Ahmed]

@VBBV code is runing but the results are incorrect as below, figure should be like in the question:

[VBBV]

Ok, fif they are arrays of 101 x 1 size then do the following change to this line

F_p = prctile(X_p(:,k),[D_p(k,:)],2);   

[Ahmed]

@VBBV This my data attached and following is code, it running but the outcome - results in the figures are wrong.

IP_p = prctile(IP_OPT,[0.5 99.5],2); %IP0_p = prctile(ip0in,[0.5 99.5],2);

IS_p = prctile(IS_OPT,[0.5 99.5],2); %IS0_p = prctile(is0in,[0.5 99.5],2);

RH_p = prctile(RH_OPT,[0.5 99.5],2); %RH0_p = prctile(rho0in,[0.5 99.5],2);

X_p = [IP_p, IS_p, RH_p];

D_p = [0.5 99.5;0.5 99.5;0.5 99.5];

for k = 1:3

F_p = prctile(X_p(:,k),[D_p(k,:)],2);

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;

q = X_p(k);

[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

[VBBV]

q = X_p(k);

Also, change the above line as

q = X_p(:,k);

i forgot this line calls each of the vectors of [IP_p, IS_p, RH_p] from X_p variable. Check whether it gives correct results now

[VBBV]

Where do you use h2 in your code ? elswhere in other code snippets ?

h2 = figure()

[Ahmed]

@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.

Answer 1

[VBBV]

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

[Ahmed]

This is working but only for first figure (in the case of IP_OPT and ip) but not for (IS_OPT and is, RH_OPT and rho). Ths because, d-obs1 is equal to only ip i.e.,

d_obs1 = ip;

How to put another loop or include d-obs1 and

cp_indx = (upper>=d_obs1)&(d_obs1>=lower);

for each e.g., d_obs1 = [ip is rho]

[VBBV]

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]

@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?

Answer 2

[Torsten]

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.

[Ahmed]

did not understnd you answer.

[Torsten]

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.