  After reviewing a large body of literature on the modeling of bivariate
discrete distributions with finite support, \cite{Gee20} made a compelling case
for the use of $I$-projections in the sense of \cite{Csi75} as a sound way to
attempt to decompose a bivariate probability mass function (p.m.f.) into its
two univariate margins and a bivariate p.m.f.\ with uniform margins playing the
role of a discrete copula. From a practical perspective, the necessary
$I$-projections on Fr\'echet classes can be carried out using the iterative
proportional fitting procedure (IPFP), also known as Sinkhorn's algorithm or
matrix scaling in the literature. After providing conditions under which a
bivariate p.m.f.\ can be decomposed in the aforementioned sense, we
investigate, for starting bivariate p.m.f.s with rectangular supports,
nonparametric and parametric estimation procedures as well as goodness-of-fit
tests for the underlying discrete copula. Related asymptotic results are
provided and build upon a differentiability result for $I$-projections on
Fr\'echet classes which can be of independent interest. Theoretical results are
complemented by finite-sample experiments and a data example.
