Suppose you have a discrete random variable, $Y$, with a large number (say, $300$) of discrete (which happen to be nominal) possible outcomes. The mass function, $p(y)=P(Y=y)$, is unknown but a sample is taken from $Y$ from 20 populations $i=1,2,...,20$. Let $Y_{ik}$ be the $k$'th sample from the $i$'th population and $n_{i}$ be the sample size taken from population $i$. Then we have the empirical mass functions:
$$\hat{f}_{i}(y) = \frac{1}{n_i} \sum_{k=1}^{n_{i}} \mathcal{I}\{Y_{ik} = y\}$$
I'm trying to figure out a way of quantifying (and possibly testing) the difference in the distribution between populations. I've thought of one possible way to do this but I'd like to hear any other suggestions, even if they are completely different.
What I've done so far is to create a measure for discrepancy between $\hat{f}_{i}$ and $\hat{f}_{j}$. One could treat this has a vector in 300-dimensional space and quantify differences by some metric (e.g. Euclidean distance). I wasn't too comfortable with this idea so, I create a matrix
$$ D_{ij} = {\rm KS}(\hat{f}_{i}, \hat{f}_{j}) $$
where ${\rm KS}$ denote the Kolmogorov-Smirnov statistic. I then use $D$ as the distance matrix in multidimensional scaling, project the data down to 2 or 3 dimensions, and test heterogeneity between the points in this reduced dimensional space.
This was a very ad-hoc approach and I'm guessing this type of problem has been worked on before so any advice or references would be much appreciated.
Note: The sample sizes are small and vary between populations. The smallest populations may only have 5-10 samples (so there are MANY zero cells), while the larger ones may have more like 30.