Clear Filters
Clear Filters

where is error in my code to ompute Prediction interval coverage probability (PICP)

2 views (last 30 days)
Following is my code and I have also attached the data. I want to compute Prediction interval coverage probability (PICP) and the graph has zero values for the y axis. Where I am doing wrong.
IP_p = prctile(IP_OPT,[2.5 97.5],2);
%%
theoreticalCoverage = linspace(0, 1, 100); % Adjust the number of points as needed
actualCoverage = zeros(size(theoreticalCoverage));
trueModel = ip;
numRealizations = 500;
% Loop through each realization
for i = 1:numRealizations
% Compute prediction interval for the ith realization
lowerLimit = IP_p(:, 1); % Access the first column (lower limit)
upperLimit = IP_p(:, 2); % Access the second column (upper limit)
% Check if true model falls within the interval
withinInterval = (trueModel >= lowerLimit(i)) & (trueModel <= upperLimit(i));
% Update actual coverage probability
actualCoverage = actualCoverage + withinInterval;
end
% Normalize actual coverage probability
actualCoverageProbability = actualCoverage / numRealizations;
% Plotting
figure;
plot(theoreticalCoverage, actualCoverageProbability, 'o-', 'LineWidth', 2);
xlabel('Theoretical Coverage Probability');
ylabel('Actual Coverage Probability');
title('Prediction Interval Coverage Probability');
ylim([0 1])
grid on;

Answers (1)

Paras Gupta
Paras Gupta on 22 Jun 2024
Hi Ahmed,
I understand you want to calculate the Prediction Interval Coverage Probability (PICP) and ensure the graph is plotted correctly. Below is the revised code that addresses the issues you encountered:
IP_p = prctile(IP_OPT, [2.5 97.5], 2);
theoreticalCoverage = linspace(0, 1, 101); % Adjust the number of points as needed
actualCoverage = zeros(size(theoreticalCoverage));
trueModel = ip;
numRealizations = 500;
% Loop through each realization
for i = 1:numRealizations
% Compute prediction interval for the ith realization
lowerLimit = IP_p(:, 1); % Access the first column (lower limit)
upperLimit = IP_p(:, 2); % Access the second column (upper limit)
% Check if true model falls within the interval for each prediction point
withinInterval = (trueModel >= lowerLimit) & (trueModel <= upperLimit);
% Update actual coverage probability
actualCoverage = actualCoverage + withinInterval';
end
% Normalize actual coverage probability
actualCoverageProbability = actualCoverage / numRealizations;
% Plotting
figure;
plot(theoreticalCoverage, actualCoverageProbability, 'o-', 'LineWidth', 2);
xlabel('Theoretical Coverage Probability');
ylabel('Actual Coverage Probability');
title('Prediction Interval Coverage Probability');
ylim([0 1]);
grid on;
The following changes were made when compared with the initial code:
  • Adjusted the initialization of 'theoreticalCoverage' to ensure correct size of (1, 101).
  • Took the transpose of 'withinInterval' to ensure that it is of the same size as 'actualCoverage'.
  • Corrected the loop to avoid indexing 'lowerLimit' and 'upperLimit' by 'i', as they do not change across realizations.
Moreover, since the lower and upper bounds are the same across all realizations and we are only interested in whether the true values fall within these bounds, the for loop in the above code can be removed and the computation can be done only for one iteration.
Hope this helps.

Community Treasure Hunt

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

Start Hunting!