function [pvalue Zvalue] = stat_proportion(g0, g1);
%function [pvalue Zvalue] = stat_proportion(g0, g1);
%
% Computes pvalue, Zvalue of testing equilvance of proportions in
% 2D biniry images. If g1 is empty, if will simply determine the
% significance of g0 only. 
% 
% g0 = 3D matrix (image dimension x image dimension x # of images)
% g1 = 3D matrix (image dimension x image dimension x # of images)
%
% Instead of using normcdf, we use erf. 
% normcdf(x) = 1/sqrt(2pi) * \int_{-infinity}^x exp(-t^2/2) dt
% erf(x) = 2/sqrrt(pi) \int_0^x exp(-t^2) dt
% By the change variable  t = s/sqrt(2), we have
% erf(x) = 2/sqrt(2pi) \int_0^{sqrt(2)x} exp(-s^2/2) dt
% erf(x/sqrt(2))/2 = 1/sqrt(2pi) \int_0^x exp(-t^2/2) dt
% Hence, normcdf(x) = 0.5 + erf(x/sqrt(2))/2.
% ex.  0.5 = normcdf(0.5) = 0.5 + erf(0/sqrt(2))/2
%       0.975 = normcdf(1.96) = 0.5 + erf(1.96/sqrt(2))/2
%       0.025 = normcdf(-1.96) = 0.5 + erf(-1.96/sqrt(2))/2
%
% (C) 2012 Moo K. Chung
% University of Wisconsin-Madison
% mkchung@wisc.edu
%
% Version 1.0 May 12, 2012 created
% Version 1.2 June 27, 2013 typo fixed

    p0=mean(g0, 3); % mean proportion
    n0=size(g0,3); % number of images
    [m1 m2] =size(g0(:,:,1)); % image dimension

if ~isempty(g1)
    % two sample test
    p1=mean(g1, 3); % mean proportion
    n1 = size(g1, 3); % number of images
    
    maxc = max(max(max(p0)), max(max(p1)));  %max proportion
    
    figure;
    subplot(2,2,1); imagesc(p0); axis off; caxis([0 maxc]); colorbar
    title('Group 1')
    
    
    subplot(2,2,2); imagesc(p1); axis off; caxis([0 maxc]); colorbar
    title('Group 2')
    
    %Z-statistic on group difference
    %Zvalue = (p0-p1)./sqrt( p0.*(1-p0)/n0 + p1.*(1-p1)); OLD CODE
    Zvalue = (p0-p1)./sqrt( p0.*(1-p0)/n0 + p1.*(1-p1)/n1);
    
    subplot(2,2,3); imagesc(Zvalue); axis off; colorbar
    title('Z-statistic (Group 1- Group2)')
    
    %p-value on group difference
    pvalue =zeros(m1*m2,1);
    Zvalue=reshape(Zvalue, m1*m2,1);
    ind =find(Zvalue>= 0);
    pvalue(ind) = 1-normcdf(Zvalue(ind));
    
    ind =find(Zvalue< 0);
    %pvalue(ind) = normcdf(Zvalue(ind));  OLD CODE
    x=Zvalue(ind);
    pvalue(ind) = 0.5 + erf(x/sqrt(2))/2;
    pvalue=reshape(pvalue, m1,m2);
    Zvalue=reshape(Zvalue, m1,m2);
       
    % only region where pvalue is less 0.05 (uncorrected) is shown
    ind = find(pvalue > 0.1);
    pvalue(ind)=0.1;
    
    ind = find(pvalue == 0);
    pvalue(ind)=0.1;
  
    subplot(2,2,4); imagesc(pvalue); axis off; colorbar
    title('p-value');
    
    colormap('hot');
    
else
    % one sample test
    figure;
    subplot(2,2,1); imagesc(p0); axis square; axis off;  colorbar
    title('Group 2')
      
    %Z-statistic computation
    % it is only testing H0: p = 0.5
    pc=0.1; 
    Zvalue = (p0-pc)./sqrt( pc.*(1-pc)/n0);
    
    subplot(2,2, 2); imagesc(Zvalue); axis square; axis off; colorbar
    title('Z-statistic')
    
    %p-value computation
    pvalue =zeros(m1*m2,1);
    Zvalue=reshape(Zvalue, m1*m2,1);
    ind =find(Zvalue>= 0);
    %pvalue(ind) = 1-normcdf(Zvalue(ind)); OLD CODE
    x=Zvalue(ind);
    pvalue(ind) =0.5 - erf(x/sqrt(2))/2;
    
        
    Zvalue=reshape(Zvalue, m1,m2);
    
    ind = find(pvalue > 0.1);
    pvalue(ind)=0.1;
    ind = find(pvalue == 0);
    pvalue(ind)=0.1;
    pvalue=reshape(pvalue, m1,m2);
  
    subplot(2,2,3); imagesc(pvalue);  axis square; axis off; colorbar
    title('p-value')
    
    colormap('hot');
end;








