catstats.m

function [POD, FAR, POFD, FBI, ACC, CSI, HRand, HSS] = catstats(ref, pred)

%------------------------------------------------------------------------------
% catstats computes categorical statistics between two 3D datasets
%------------------------------------------------------------------------------
% 
% This function calculates various categorical statistics to compare two 
% datasets (e.g., observations vs. model predictions) on a pixel-wise basis.
% It computes:
%   - Probability of Detection (POD)
%   - False Alarm Ratio (FAR)
%   - Probability of False Detection (POFD)
%   - Frequency Bias Index (FBI)
%   - Accuracy (ACC)
%   - Critical Success Index (CSI)
%   - Hanssen and Kuipers Discriminant (HRand)
%   - Heidke Skill Score (HSS)
%
% Inputs:
%   ref: 3D numeric matrix (latitude x longitude x time) [Reference data]
%   pred: 3D numeric matrix (latitude x longitude x time) [Predicted data]
%
% Outputs:
%   POD   - Probability of Detection (Hits / Hits + Misses)
%   FAR   - False Alarm Ratio (False Alarms / Hits + False Alarms)
%   POFD  - Probability of False Detection (False Alarms / False Alarms + Correct Rejections)
%   FBI   - Frequency Bias Index ((Hits + False Alarms) / (Hits + Misses))
%   ACC   - Accuracy ((Hits + Correct Rejections) / Total)
%   CSI   - Critical Success Index (Hits / Hits + Misses + False Alarms)
%   HRand - Hanssen and Kuipers Discriminant ((Hits + Misses) * (Hits + False Alarms) / Total)
%   HSS   - Heidke Skill Score
%
%------------------------------------------------------------------------------

    arguments
        ref {mustBeNumeric, mustBe3D}
        pred {mustBeNumeric, mustBe3D, mustBeSameSize(ref, pred)}
    end

    % Vectorized computation
    a = sum(ref(:) > 0 & pred(:) > 0);
    b = sum(ref(:) == 0 & pred(:) > 0);
    c = sum(ref(:) > 0 & pred(:) == 0);
    d = sum(ref(:) == 0 & pred(:) == 0);

    % Compute statistics
    POD = a / (a + c);
    FAR = b / (a + b);
    POFD = b / (b + d);
    FBI = (a + b) / (a + c);
    ACC = (a + d) / (a + b + c + d);
    CSI = a / (a + b + c);
    HRand = (a + c) * (a + b) / (a + b + c + d);
    HSS = 2 * (a * d - b * c) / ((a + c) * (c + d) + (a + b) * (b + d));

end

% Custom validation function to check if input is 3D
function mustBe3D(data)
    if ndims(data) ~= 3
        error("Input must be a 3D matrix (latitude x longitude x time)");
    end
end

% Custom validation function to check if inputs are the same size
function mustBeSameSize(ref, pred)
    if ~isequal(size(ref), size(pred))
        error("Both datasets must have the same dimensions.");
    end
end