function [dp1,dp2] = kafc_dprime(resp,gt,maxcat)

    % method for computing d-prime for any number of categories. Note that,
    % in this algorithm, we compute N values of d-prime. This method gives 
    % better results than the method mafc_dprime when there are a large 
    % number of responses. This method nonetheless gives results that are
    % extremely similar to mafc_dprime_2 (R^2 = ~1);

    % % % % % % % % % % % %         INPUT       % % % % % % % % % % % %
    
    % resp = column vector containing human responses
    % gt = column vector containing ground-truth labels
    % maxcat = number of categories in category system
    
    % % % % % % % % % % % %         OUTPUT       % % % % % % % % % % % %
    
    % dp1 = estimate of d-prime (weighted according to number of gt cats)
    % dp2 = estimate of d-prime (uniformly weighted)
    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
    
    cats = 1:maxcat;
    nc = length(cats);
    h = zeros(1,nc);
    fa = zeros(1,nc);
    p = zeros(1,nc);
    
    % skip calculation if fixed response. We know this should be d'=0
    if length(unique(resp)) == 1
        dp1 = 0;
        dp2 = 0;
    else

        for i = 1:nc
            hidx = gt == cats(i);        
            p(i) = mean(hidx);

            nh = sum(resp(hidx) == cats(i));
            nht = sum(hidx);

            nfa = sum(resp(~hidx) == cats(i));
            nfat = sum(~hidx);
            
            % correction for no hits: assume random response 
            if ismember(0,[nh,nht]) || isequal(nh,nht) 
                nh = nh+(1/maxcat);
                nht = nht+1;
            end
            
            % correction for no fas: assume random response 
            if ismember(0,[nfa,nfat]) || isequal(nfa,nfat)
                nfa = nfa+(1/maxcat);
                nfat = nfat+1;
            end

            h(i) = nh/nht;
            fa(i) = nfa/nfat;
        end
        
        % weighted by frequency of ground-truth categories 
        dp1 = sum((norminv(h)-norminv(fa)).*p);

        % weighted uniformly
        dp2 = mean(norminv(h) - norminv(fa));
        
    end

end