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.
