  When the target of inference is a real-valued function of probability
parameters in the k-sample multinomial problem, variance estimation may be
challenging. In small samples, methods like the nonparametric bootstrap or
delta method may perform poorly. We propose a novel general method in this
setting for computing exact p-values and confidence intervals which means that
type I error rates are correctly bounded and confidence intervals have at least
nominal coverage at all sample sizes. Our method is applicable to any
real-valued function of multinomial probabilities, accommodating an arbitrary
number of samples with varying category counts. We describe the method and
provide an implementation of it in R, with some computational optimization to
ensure broad applicability. Simulations demonstrate our method's ability to
maintain correct coverage rates in settings where the nonparametric bootstrap
fails.
