function [si,mri,sdt,key] = grt_tests(fM,ssp);

% function [si,mri,sdt] = grt_tests(fM,ssp);
%
% si        results of chi-squared tests of (sampling) independence
%               one test for each row of fM
% mri       results of proportion tests of marginal response invariance
% sdt       signal detection theory parameters (d', B) for macroanalyses
%
% fM        frequency (confusion) matrix
% ssp       stimulus level specification (default = [1 1; 1 2; 2 1; 2 2]

% correction for zeros (should be same as model fitting, z-roc analyses)
if sum(sum(sum(fM==0)))>=1
    fM = fM + .25;
end

[rows cols mats] = size(fM);

% key to rows and columns of output arrays
key = cell(1);

if nargin < 2
    ssp = [1 1; 1 2; 2 1; 2 2];
end

% key for chi-square tests
if mats > 1
    key(1) = {'rows = stimuli (aa,ab,ba,bb); cols = chi-square, p-value; array = condition'};
else
    key(1) = {'rows = stimuli (aa,ab,ba,bb); cols = chi-square, p-value'};
end


si = zeros(mats,rows);
for ci = 1:mats
    
    % make this test for SI in each condition...
    
    for chii = 1:4
        obsM = zeros(2,2);
        expM = zeros(2,2);
        for sspi = 1:4
            obsM(ssp(sspi,1),ssp(sspi,2)) = fM(chii,sspi,ci);
        end
        rowmarg = sum(obsM,2);
        colmarg = sum(obsM,1);
        N = sum(sum(obsM));
        for ei = 1:2
            for ej = 1:2
                expM(ei,ej) = rowmarg(ei)*colmarg(ej)/N;
            end
        end
        si(ci,chii) = sum(sum( ( (obsM - expM).^2 ) ./ expM ) );
    end
end
    
% chi-square based tests of SI
bin = 0.01;
%x = bin:bin:10;
%gam = bin*sum(exp(-x).*(x.^-.5));
%x2 = (x.^-.5).*exp(-x./2);%./(gam*(2^.5));
%x2 = x2/(bin*sum(x2));
%x3 = (x.^-.5).*exp(-x./2)./(gam*(2^.5));
%for ci = 1:length(x)
%    cx2(ci) = bin*sum(x2(1:ci));
    %cx3(ci) = bin*sum(x3(1:ci));
%end


% for ci = 1:mats
%     for chii = 1:4
%         if si(chii,1,ci) > x(end)
%             si(chii,2,ci) = abs(1-cx2(end));
%         else
%             [dif, xci] = min(abs(si(chii,1,ci)-x));
%             si(chii,2,ci) = 1-cx2(xci);
%         end
%     end
% end

% z-tests of MRI

key(2,1) = {'rows = tests for mri at a1, a2, b1, b2; cols = conditions'};
mri = zeros(mats,rows);
for ci = 1:mats

    % p( A1 | A1B1 ) ?= p( A1 | A1B2 )
    n_a = sum(fM(1,:,ci));
    p_a = ( fM(1,1,ci) + fM(1,2,ci) )/n_a;
    n_b = sum(fM(2,:,ci));
    p_b = ( fM(2,1,ci) + fM(2,2,ci) )/n_b;
    p_c = ( fM(1,1,ci) + fM(1,2,ci) + fM(2,1,ci) + fM(2,2,ci) )/sum([n_a n_b]);
    z_a1 = ( p_a - p_b ) / sqrt( (p_c*(1-p_c))*( 1/n_a + 1/n_b ) );
    mri(ci,1) = z_a1;
    
    % p( A2 | A2B1 ) ?= p( A2 | A2B2 )
    n_a = sum(fM(3,:,ci));
    p_a = ( fM(3,3,ci) + fM(3,4,ci) )/n_a;
    n_b = sum(fM(4,:,ci));
    p_b = ( fM(4,3,ci) + fM(4,4,ci) )/n_b;
    p_c = ( fM(3,3,ci) + fM(3,4,ci) + fM(4,3,ci) + fM(4,4,ci) )/sum([n_a n_b]);
    z_a2 = ( p_a - p_b ) / sqrt( (p_c*(1-p_c))*( 1/n_a + 1/n_b ) );
    mri(ci,2) = z_a2;

    % p( B1 | A1B1 ) ?= p( B1 | A2B1 )
    n_a = sum(fM(1,:,ci));
    p_a = ( fM(1,1,ci) + fM(1,3,ci) )/n_a;
    n_b = sum(fM(3,:,ci));
    p_b = ( fM(3,1,ci) + fM(3,3,ci) )/n_b;
    p_c = ( fM(1,1,ci) + fM(1,3,ci) + fM(3,1,ci) + fM(3,3,ci) )/sum([n_a n_b]);
    z_a3 = ( p_a - p_b ) / sqrt( (p_c*(1-p_c))*( 1/n_a + 1/n_b ) );
    mri(ci,3) = z_a3;

    % p( B2 | A1B2 ) ?= p( B2 | A2B2 )
    n_a = sum(fM(2,:,ci));
    p_a = ( fM(2,2,ci) + fM(2,4,ci) )/n_a;
    n_b = sum(fM(4,:,ci));
    p_b = ( fM(4,2,ci) + fM(4,4,ci) )/n_b;
    p_c = ( fM(2,2,ci) + fM(2,4,ci) + fM(4,2,ci) + fM(4,4,ci) )/sum([n_a n_b]);
    z_a4 = ( p_a - p_b ) / sqrt( (p_c*(1-p_c))*( 1/n_a + 1/n_b ) );
    mri(ci,4) = z_a4;

end


% sdt based tests of PS, DS

x = -4:bin:4;
x4 = (1/sqrt(2*pi))*exp(-.5*x.^2);
for ci = 1:length(x)
    cx4(ci) = bin*sum(x4(1:ci));
end

pM = zeros(rows,cols,mats);
for ci = 1:mats
    for r = 1:rows
        pM(r,:,ci) = fM(r,:,ci)/sum(fM(r,:,ci));
    end
end

miller = 1;

sdt = zeros(mats,8,2);
for ci = 1:mats
    
    % d'(A,B1)
    p_ht = pM(3,3,ci) + pM(3,4,ci);
    p_fa = pM(1,3,ci) + pM(1,4,ci);
    ht_dif = abs( p_ht - cx4 );
    [ht_min,ht_min_idx] = min(ht_dif);
    fa_dif = abs( p_fa - cx4 );
    [fa_min,fa_min_idx] = min(fa_dif);
    z_ht = x(ht_min_idx);
    z_fa = x(fa_min_idx);
    d_prime = z_ht - z_fa;
    c_star = -z_ht;
    beta = exp( 0.5*(z_fa^2 - c_star^2) );
    log_beta = log(beta);
    c_prime = log_beta/d_prime;
    if miller
        % hits
        N_h = sum(fM(3,:,ci))+1;
        ezph = zeros(N_h,1);
        ezph2 = zeros(N_h,1);
        mu = N_h*p_ht;
        sg = sqrt(N_h*p_ht*(1-p_ht));
        for k = 1:N_h
            kN_dif = abs(((k-1)/N_h)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezph(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezph2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        N_f = sum(fM(1,:,ci))+1;
        ezpf = zeros(N_f,1);
        ezpf2 = zeros(N_f,1);
        mu = N_f*p_fa;
        sg = sqrt(N_f*p_fa*(1-p_fa));
        for k = 1:N_f
            kN_dif = abs(((k-1)/N_f)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezpf(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezpf2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        var_zh = sum(ezph2) - sum(ezph)^2;
        var_zf = sum(ezpf2) - sum(ezpf)^2;
        var_d = var_zh + var_zf;
    else
        var_d = (p_ht*(1-p_ht))/(sum(fM(3,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_ht^2)^2)) + ...
            (p_fa*(1-p_fa))/(sum(fM(1,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_fa^2)^2));
    end
    var_c = .25*var_d;
    sdt(ci,1,1) = d_prime;
    sdt(ci,1,2) = sqrt(var_d);
    %sdt(1,2,ci) = beta;
    %sdt(1,3,ci) = log_beta;
    sdt(ci,3,1) = c_prime;
    sdt(ci,3,2) = sqrt(var_c);
    
    % d'(A,B2)
    p_ht = pM(4,3,ci) + pM(4,4,ci);
    p_fa = pM(2,3,ci) + pM(2,4,ci);
    ht_dif = abs( p_ht - cx4 );
    [ht_min,ht_min_idx] = min(ht_dif);
    fa_dif = abs( p_fa - cx4 );
    [fa_min,fa_min_idx] = min(fa_dif);
    z_ht = x(ht_min_idx);
    z_fa = x(fa_min_idx);
    d_prime = z_ht - z_fa;
    c_star = -z_ht;
    beta = exp( 0.5*(z_fa^2 - c_star^2) );
    log_beta = log(beta);
    c_prime = log_beta/d_prime;
    if miller
        % hits
        N_h = sum(fM(4,:,ci))+1;
        ezph = zeros(N_h,1);
        ezph2 = zeros(N_h,1);
        mu = N_h*p_ht;
        sg = sqrt(N_h*p_ht*(1-p_ht));
        for k = 1:N_h
            kN_dif = abs(((k-1)/N_h)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezph(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezph2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        N_f = sum(fM(2,:,ci))+1;
        ezpf = zeros(N_f,1);
        ezpf2 = zeros(N_f,1);
        mu = N_f*p_fa;
        sg = sqrt(N_f*p_fa*(1-p_fa));
        for k = 1:N_f
            kN_dif = abs(((k-1)/N_f)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezpf(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezpf2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        var_zh = sum(ezph2) - sum(ezph)^2;
        var_zf = sum(ezpf2) - sum(ezpf)^2;
        var_d = var_zh + var_zf;
    else
        var_d = (p_ht*(1-p_ht))/(sum(fM(4,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_ht^2)^2)) + ...
            (p_fa*(1-p_fa))/(sum(fM(2,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_fa^2)^2));
    end
    var_c = .25*var_d;
    sdt(ci,2,1) = d_prime;
    sdt(ci,2,2) = sqrt(var_d);
    %sdt(2,2,ci) = beta;
    %sdt(2,3,ci) = log_beta;
    sdt(ci,4,1) = c_prime;
    sdt(ci,4,2) = sqrt(var_c);
    
    % d'(B,A1)
    p_ht = pM(2,2,ci) + pM(2,4,ci);
    p_fa = pM(1,2,ci) + pM(1,4,ci);
    ht_dif = abs( p_ht - cx4 );
    [ht_min,ht_min_idx] = min(ht_dif);
    fa_dif = abs( p_fa - cx4 );
    [fa_min,fa_min_idx] = min(fa_dif);
    z_ht = x(ht_min_idx);
    z_fa = x(fa_min_idx);
    d_prime = z_ht - z_fa;
    c_star = -z_ht;
    beta = exp( 0.5*(z_fa^2 - c_star^2) );
    log_beta = log(beta);
    c_prime = log_beta/d_prime;
    if miller
        % hits
        N_h = sum(fM(2,:,ci))+1;
        ezph = zeros(N_h,1);
        ezph2 = zeros(N_h,1);
        mu = N_h*p_ht;
        sg = sqrt(N_h*p_ht*(1-p_ht));
        for k = 1:N_h
            kN_dif = abs(((k-1)/N_h)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezph(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezph2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        N_f = sum(fM(1,:,ci))+1;
        ezpf = zeros(N_f,1);
        ezpf2 = zeros(N_f,1);
        mu = N_f*p_fa;
        sg = sqrt(N_f*p_fa*(1-p_fa));
        for k = 1:N_f
            kN_dif = abs(((k-1)/N_f)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezpf(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezpf2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        var_zh = sum(ezph2) - sum(ezph)^2;
        var_zf = sum(ezpf2) - sum(ezpf)^2;
        var_d = var_zh + var_zf;
    else
        var_d = (p_ht*(1-p_ht))/(sum(fM(2,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_ht^2)^2)) + ...
            (p_fa*(1-p_fa))/(sum(fM(1,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_fa^2)^2));
    end
    var_c = .25*var_d;
    sdt(ci,5,1) = d_prime;
    sdt(ci,5,2) = sqrt(var_d);
    %sdt(3,2,ci) = beta;
    %sdt(3,3,ci) = log_beta;
    sdt(ci,7,1) = c_prime;
    sdt(ci,7,2) = sqrt(var_c);

    % d'(B,A2)
    p_ht = pM(4,2,ci) + pM(4,4,ci);
    p_fa = pM(3,2,ci) + pM(3,4,ci);
    ht_dif = abs( p_ht - cx4 );
    [ht_min,ht_min_idx] = min(ht_dif);
    fa_dif = abs( p_fa - cx4 );
    [fa_min,fa_min_idx] = min(fa_dif);
    z_ht = x(ht_min_idx);
    z_fa = x(fa_min_idx);
    d_prime = z_ht - z_fa;
    c_star = -z_ht;
    beta = exp( 0.5*(z_fa^2 - c_star^2) );
    log_beta = log(beta);
    c_prime = log_beta/d_prime;
    if miller
        % hits
        N_h = sum(fM(4,:,ci))+1;
        ezph = zeros(N_h,1);
        ezph2 = zeros(N_h,1);
        mu = N_h*p_ht;
        sg = sqrt(N_h*p_ht*(1-p_ht));
        for k = 1:N_h
            kN_dif = abs(((k-1)/N_h)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezph(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezph2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        N_f = sum(fM(3,:,ci))+1;
        ezpf = zeros(N_f,1);
        ezpf2 = zeros(N_f,1);
        mu = N_f*p_fa;
        sg = sqrt(N_f*p_fa*(1-p_fa));
        for k = 1:N_f
            kN_dif = abs(((k-1)/N_f)-cx4);
            [kN_min,kN_idx] = min(kN_dif);
            z_kN = cx4(kN_idx);
            ezpf(k) = z_kN*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
            ezpf2(k) = (z_kN^2)*((1/(sqrt(2*pi)*sg))*exp(-.5*((k-mu)/sg)^2));
        end
        var_zh = sum(ezph2) - sum(ezph)^2;
        var_zf = sum(ezpf2) - sum(ezpf)^2;
        var_d = var_zh + var_zf;
    else
        var_d = (p_ht*(1-p_ht))/(sum(fM(4,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_ht^2)^2)) + ...
            (p_fa*(1-p_fa))/(sum(fM(3,:,ci))*((1/(sqrt(2*pi)))*exp(-.5*z_fa^2)^2));
    end
    var_c = .25*var_d;
    sdt(ci,6,1) = d_prime;
    sdt(ci,6,2) = sqrt(var_d);
    %sdt(4,2,ci) = beta;
    %sdt(4,3,ci) = log_beta;
    sdt(ci,8,1) = c_prime;    
    sdt(ci,8,2) = sqrt(var_c);
end