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)
%
% (C) 2012 Moo K. Chung
% University of Wisconsin-Madison
% mkchung@wisc.edu
%

    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 0')
    
    
    subplot(2,2,2); imagesc(p1); axis off; caxis([0 maxc]); colorbar
    title('Group 1')
    
    %Z-statistic on group difference
    Zvalue = (p0-p1)./sqrt( p0.*(1-p0)/n0 + p1.*(1-p1));
    
    subplot(2,2,3); imagesc(Zvalue); axis off; colorbar
    title('Z-statistic (Group 0- Group1)')
    
    %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));  
    pvalue=reshape(pvalue, m1,m2);
    
    Zvalue=reshape(Zvalue, m1,m2);
       
    % only region where pvalue is less 0.05 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(1,3,1); imagesc(p0); axis off;  colorbar
    title('Group 0')
      
    %Z-statistic 
    % it is only testing H0: p = 0.5
    pc=0.1; 
    Zvalue = (p0-pc)./sqrt( pc.*(1-pc)/n0);
    
    subplot(1,3, 2); imagesc(Zvalue); axis off; colorbar
    title('Z-statistic')
    
    %p-value 
    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));  
    
    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(1, 3,3); imagesc(pvalue);  axis off; colorbar
    title('p-value')
    
    colormap('hot');
end;






