
% function takes a contingeny matrix and returns
% chi-square and p-values 

function [chi2, pval] = chi_square_pval (cont_mat)
%compute row, column, and population sums
row_sum = sum(cont_mat, 1);
col_sum = sum(cont_mat, 2);
pop_sum = sum(sum(cont_mat, 1), 2);

%duplicate row/column sums for all rows/columns 
row_mat = repmat(row_sum,  size(cont_mat, 1), 1);
col_mat = repmat(col_sum,  1, size(cont_mat, 2));

% or simply use meshgrid
[row_mat, col_mat] = meshgrid(row_sum, col_sum);

%compute expected observations

expect  = row_mat .* col_mat / pop_sum;

%compute chi square values and sum matrix 
chi2_mat = (cont_mat - expect).^ 2 ./ expect;
chi2 = sum(sum(chi2_mat, 1), 2);

% compute the p-value using df = (r-1)*(c-1)
r = size(cont_mat, 1);
c = size(cont_mat, 2);
pval = 1.0 - chi2cdf(chi2, (r - 1) * (c - 1));
