Just for the record, there is an option to the
chisq.test function in R that allows to simulate the p-value by randomly generate a given number of independent tables instead of deriving it from the chi-squared distribution :
chisq.test(x, y, simulate.p.value=TRUE, B=2000)
From the function help page :
If ‘simulate.p.value’ is ‘FALSE’, the p-value is computed from the
asymptotic chi-squared distribution of the test statistic;
continuity correction is only used in the 2-by-2 case (if
‘correct’ is ‘TRUE’, the default). Otherwise the p-value is
computed for a Monte Carlo test (Hope, 1968) with ‘B’ replicates.
In the contingency table case simulation is done by random
sampling from the set of all contingency tables with given
marginals, and works only if the marginals are strictly positive.
(A C translation of the algorithm of Patefield (1981) is used.)
Continuity correction is never used, and the statistic is quoted
without it. Note that this is not the usual sampling situation
assumed for the chi-squared test but rather that for Fisher's
That could be a way to estimate a p-value without worrying about expected counts.