  We show that discrete distributions on the $d$-dimensional non-negative
integer lattice can be approximated arbitrarily well via the marginals of
stationary distributions for various classes of stochastic chemical reaction
networks. We begin by providing a class of detailed balanced networks and prove
that they can approximate any discrete distribution to any desired accuracy.
However, these detailed balanced constructions rely on the ability to
initialize a system precisely, and are therefore susceptible to perturbations
in the initial conditions. We therefore provide another construction based on
the ability to approximate point mass distributions and prove that this
construction is capable of approximating arbitrary discrete distributions for
any choice of initial condition. In particular, the developed models are
ergodic, so their limit distributions are robust to a finite number of
perturbations over time in the counts of molecules.
