  Small-scale density fluctuations can significantly affect reionization, but
are typically modelled quite crudely. Unresolved fluctuations in numerical
simulations and analytical calculations are included using a gas
\textit{clumping factor}, typically assumed to be independent of the local
environment. In \cite{Mao2019} (Paper I) we presented an improved, local
density-dependent model for the sub-grid gas clumping. Here we extend this
using an \textit{empirical stochastic model} based on the results from
high-resolution numerical simulations which fully resolve all relevant
fluctuations. Our model reproduces well both the mean density-clumping relation
and its scatter. We applied our stochastic model, along with the mean clumping
one and the Paper I deterministic model, to create a large-volume realisations
of the clumping field, and used these in radiative transfer simulations of
cosmic reionization. Our results show that the simplistic mean clumping model
delays reionization compared to local density-dependent models, despite
producing fewer recombinations overall. This is due to the very different
spatial distribution of clumping, resulting in much higher photoionization
rates in the latter cases. The mean clumping model produces smaller H~II
regions throughout most of reionization, but those percolate faster at late
times. It also causes significant delay in the 21-cm fluctuations peak and
yields lower non-Gaussianity and many fewer bright pixels in the PDF
distribution. The stochastic density-dependent model shows relatively minor
differences from the deterministic one, mostly concentrated around overlap,
where it significantly suppresses the 21-cm fluctuations, and at the bright
tail of the 21-cm PDFs, where it produces noticeably more bright pixels.
