A geometric combination estimator is proposed for $d$-dimensional ordinal contingency tables. The proposed estimator is nonnegative. It is shown that, assuming sufficient smoothness and boundary conditions for the underlying probabilities, the rate of convergence of mean summed squared error (MSSE) of this estimator is $O(K^{-1}N^{-8/(d+8)})$ for $d$-dimensional tables $(d \leq 4)$ with $K$ cells and sample size $N$. This rate is optimal under the smoothness assumptions, and is faster than that attained by nonnegative kernel estimates. Boundary kernels for multidimensional tables are also developed for the proposed estimator to relax restrictive boundary conditions, resulting in summed squared error (SSE) being of order $O_p(K^{-1}N^{-8/(d+8)})$ for all $d \geq 1$. The behavior of the new estimator is investigated through simulations and applications to real data. It is shown that even for relatively small tables, these estimators are superior to nonnegative kernel estimators, in sharp contrast to the relatively unimpressive performance of such estimators for continuous data.