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

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;

Answer 1

[Paras Gupta]

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.