measure the fit of two distributions The population of a county is 150,000, ages and sexes are known. Somebody selected a sample of 2,000 people, the selection rules are not quite clear but the age and sex distributions are given. My objective is to find out what selection rules were used. I built a model which does this, and have my own sample of about the same size. The age and gender profiles are close to the original sample, but the fit is less than 100% perfect. I'd like to assess the quality of fit, for example, estimate the probability that a randomly selected sample would be as close to the original one as mine. The trouble is I do not know how to do this. Could you help?

Let's just look at age groups. Suppose there are five groups and that
the population proportions are given by vector pr=(.2, .2, .3, .2, .1).
Take two samples of size $2000,$ very much smaller than the population.
set.seed(902)
pr = c(.2,.2,.3,.2,.1)
x1 = sample(1:5, 2000, rep=T, p = pr)
x2 = sample(1:5, 2000, rep=T, p = pr)

There are relatively small differences in category counts between the tow samples.
However, treating categories as nominal, a chi-squared test finds no
significant difference between the two samples at the 5% level.
The P-value $0.55 > 0.05 = 5\%$ of the test is probability of a worse fit (according to the chi-squared
criterion, comparing observed and expected counts.).
t1 = tabulate(t1);  t2 = tabulate(t2)
TBL = rbind(t1,t2);  TBL
t1 = tabulate(x1);  t2 = tabulate(x2)
TBL = rbind(t1,t2);  TBL
   [,1] [,2] [,3] [,4] [,5]
t1  408  390  624  371  207
t2  412  366  637  399  186

par(mfrow=c(1,2))
 barplot(t1, names.arg=1:5, xlab="Sample 1", col="skyblue2")
 barplot(t2, names.arg=1:5, xlab="Sample 2", col="skyblue2")
par(mfrow=c(1,1))


chisq.test(TBL)

    Pearson's Chi-squared test

data:  TBL
X-squared = 3.0558, df = 4, p-value = 0.5485

In the figure below, the vertical dotted line is at the chi-squared statistic. The area under the density curve to the right of that line
is the P-value of the test.


Treating the age categories as ordinal, a Wilcoxon rank sum
test finds no significant difference in the locations of the two
samples.
wilcox.test(x1,x2)

        Wilcoxon rank sum test with continuity correction

data:  x1 and x2
W = 1994000, p-value = 0.8668
alternative hypothesis: true location shift is not equal to 0

boxplot(x1, x2, col="skyblue2", horizontal=T)


However, I know of no test that looks at your desired "probability that a randomly selected sample would be as close to the original one as mine."
I don't know your purpose in comparing the two samples. Perhaps a better question is whether for practical purposes there is an important difference
in counts between the two. That is a judgment call---IMHO not a statistical question.
Notes: (1) You could test whether your Sample 2 is a bad match for the
population proportion. Of course, your sample is very unlikely to be
an exact match to the population proportion, but as long as you use
random sampling, the P-value will tend to be above $0.05,$ sometimes just barely, sometimes greatly, and (of course about $5\%$ of the time) below $0.05.$ [Because the chi-squared statistic, which depends on counts, is discrete, the distribution of the P-value under $H_0$ will not be exacty uniform in $(0,1).]$
chisq.test(t2, p = pr)

        Chi-squared test for given probabilities

data:  t2
X-squared = 6.5142, df = 4, p-value = 0.1639


    Chi-squared test for given probabilities

data:  t2
X-squared = 6.5142, df = 4, p-value = 0.1639

(2) While some chi-squared tests are often called "goodness-of-fit" tests, but they reject 'fit' when the chi-squared statistic is large
and we reject when the chi-squared statistic is large (and P-value is small). So it might be better to call these chi-squared statistics "badness-of-fit" statistics--the larger the statistic, the worse the fit.
