Clear Filters
Clear Filters

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.

Categories

Find more on Vibration Analysis in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!