PSFDS Chapter3

Chapter3/myfisher23.m

@@ -0,0 +1,118 @@
+function Pvalue=myfisher23(x)
+%P=MYFISHER23(X)- Fisher's Exact Probability Test on 2x3 matrix.
+%Fisher's exact test of 2x3 contingency tables permits calculation of
+%precise probabilities in situation where, as a consequence of small cell
+%frequencies, the much more rapid normal approximation and chi-square
+%calculations are liable to be inaccurate. The Fisher's exact test involves
+%the computations of several factorials to obtain the probability of the
+%observed and each of the more extreme tables. Factorials growth quickly,
+%so it's necessary use logarithms of factorials. This computations is very
+%easy in Matlab because x!=gamma(x+1) and log(x!)=gammaln(x+1). This
+%function is fully vectorized to speed up the computation.
+%
+% Syntax: 	myfisher23(x)
+%      
+%     Inputs:
+%           X - 2x3 data matrix 
+%     Outputs:
+%           - Three p-values
+%
+%   Example:
+%
+%                A   B   C
+%           -------------------
+%      X         0   3   2
+%           -------------------    
+%      Y         6   5   1
+%           -------------------
+%
+%   Calling on Matlab the function: 
+%             myfisher23([0 3 2; 6 5 1])
+%
+%   Answer is:
+%
+% --------------------------------------------------------------------------------
+% 2x3 matrix Fisher's exact test
+% --------------------------------------------------------------------------------
+%     Tables    two_tails_p_value    Mid_p_correction
+%     ______    _________________    ________________
+% 
+%     18        0.088235             0.074661  
+%
+%           Created by Giuseppe Cardillo
+%           [email protected]
+%
+% To cite this file, this would be an appropriate format:
+% Cardillo G. (2007) MyFisher23: a very compact routine for Fisher's exact
+% test on 2x3 matrix
+% http://www.mathworks.com/matlabcentral/fileexchange/15399
+
+
+%Input error handling
+p = inputParser;
+addRequired(p,'x',@(x) validateattributes(x,{'numeric'},{'real','finite','integer','nonnegative','nonnan','size',[2 3]}));
+parse(p,x);
+clear p
+
+Rs=sum(x,2); %rows sum
+Cs=sum(x); %columns sum
+N=sum(Rs); %Total observations
+
+%If necessed, rearrange matrix
+if ~issorted(Cs)
+    [Cs,ind]=sort(Cs);
+    x=x(:,ind);
+    clear ind
+end
+if ~issorted(Rs)
+    [Rs,ind]=sort(Rs);
+    x=x(ind,:);
+    clear ind
+end
+
+%recall that Fisher's P=[Prod(Rs!)*Prod(Cs!)]/[N!*prod(X(i,j)!)]
+%Log(A*B)=Log(A)+Log(B) and Log(A/B)=Log(A)-Log(B)
+%Construct all possible tables
+%A 2x3 matrix has 2 degrees of freedom...
+A=0:1:min(Rs(1),Cs(1)); %all possible values of X(1,1)
+B=min(Cs(2),Rs(1)-A); %max value of X(1,2) given X(1,1)
+C=max(Rs(1)-A-Cs(3),0); % min value of X(1,2) given X(1,1)
+et=sum(B-C+ones(size(B))); %tables to evaluate
+Tables=zeros(et,6); %Matrix preallocation
+%compute the index
+stop=cumsum(B-C+1);
+start=[1 stop(1:end-1)+1];
+%In the first round of the for cycle, Column 1 assignment should be skipped
+%because it is already zero. So, modify the cycle...
+Tables(start(1):stop(1),2)=C(1):1:B(1); %Put in the Column2 all the possible values of X(1,2) given X(1,1)
+for I=2:length(A)
+    Tables(start(I):stop(I),1)=A(I); %replicate the A(I) value for B(I)-C(I)+1 times
+    %Put in the Column2 all the possible values of X(1,2) given X(1,1)
+    Tables(start(I):stop(I),2)=C(I):1:B(I); 
+end
+clear A B start stop
+%The degrees of freedom are finished, so complete the table...
+%...Put all the possible values of X(1,3) given X(1,1) and X(1,2)
+Tables(:,3)=Rs(1)-sum(Tables(:,1:2),2);
+%Complete the second row given the first row
+Tables(:,4:6)=repmat(Cs,et,1)-Tables(:,1:3);
+
+%Compute log(x!) using the gammaln function
+zf=gammaln(Tables+1); %compute log(x!)
+K=sum(gammaln([Rs' Cs]+1))-gammaln(N+1); %The costant factor K=log(prod(Rs!)*prod(Cs!)/N!)
+np=exp(K-sum(zf,2)); %compute the p-value of each possible matrix
+[~,obt]=ismember(x(1,:),Tables(:,1:3),'rows'); %Find the observed table
+clear zf K tf
+%Finally compute the probability for 2-tailed test
+P=sum(np(np<=np(obt)));
+
+%display results
+tr=repmat('-',1,80);
+disp(tr)
+disp('2x3 matrix Fisher''s exact test')
+disp(tr)
+disp(array2table([et,P,0.5*np(obt)+sum(np(np<np(obt)))],'VariableNames',{'Tables','two_tails_p_value','Mid_p_correction'}));
+
+if nargout
+    Pvalue=P;
+end

Chapter3/perm_func.m

@@ -0,0 +1,9 @@
+function mean_diff = perm_func(x, n1, n2)
+    n = n1 + n2;
+    perm_ind = randperm(n);
+    idx_a = perm_ind(1:n1);
+    idx_b = perm_ind(n1+1:n);
+
+    mean_diff = mean(x(idx_b)) - mean(x(idx_a));
+end
+

Chapter3/prop_test.m

@@ -0,0 +1,62 @@
+function [h,p, chi2stat,df] = prop_test(X , N, correct)
+
+% [h,p, chi2stat,df] = prop_test(X , N, correct)
+% 
+% A simple Chi-square test to compare two proportions
+% It is a 2 sided test with alpha=0.05
+%
+% Input:
+% X = vector with number of success for each sample (e.g. [20 22])
+% N = vector of total counts for each sample (e.g. [48 29])
+% correct = true/false : Yates continuity correction for small samples?
+%
+% Output:
+% h = hypothesis (H1/H0)
+% p = p value
+% chi2stat= Chi-square value
+% df = degrees of freedom (always equal to 1: 2 samples)
+%
+% Needs chi2cdf from the Statistics toolbox
+% Inspired by prop.test() in "R" but much more basic
+%
+% Example: [h,p,chi]=prop_test([20 22],[48 29], true)
+% The above example tests if 20/48 differs from 22/29 using Yate's correction
+
+if (length(X)~= 2)||(length(X)~=length(N))
+    disp('Error: bad vector length')
+elseif (X(1)>N(1))|| (X(2)>N(2))
+    disp('Error: bad counts (X>N)')
+else  
+    df=1; % 2 samples
+
+    % Observed data
+    n1 = X(1);
+    n2 = X(2);
+    N1 = N(1);
+    N2 = N(2);
+
+    % Pooled estimate of proportion
+    p0 = (n1+n2) / (N1+N2);
+
+    % Expected counts under H0 (null hypothesis)
+    n10 = N1 * p0;
+    n20 = N2 * p0;
+    observed = [n1 N1-n1 n2 N2-n2];
+    expected = [n10 N1-n10 n20 N2-n20];
+
+    if correct == false
+        % Standard Chi-square test
+        chi2stat = sum((observed-expected).^2 ./ expected);
+        p = 1 - chi2cdf(chi2stat,1);
+    else
+        % Yates continuity correction        
+        chi2stat = sum((abs(observed - expected) - 0.5).^2 ./ expected);
+        p = 1 - chi2cdf(chi2stat,1);
+    end
+
+    h=0;
+    if p<0.05
+        h=1;
+    end
+end
+end

LICENSE

@@ -1,21 +1,21 @@
-MIT License
-
-Copyright (c) 2020 PSFDS-MATLAB
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
+MIT License
+
+Copyright (c) 2020 PSFDS-MATLAB
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.