# Hypergeometric: how do I construct a credibility interval around K (population successes) in R?

I have a problem for which I believe I should use the hypergeometric distribution, but I can't figure out how to do it in R.

Say I have a bag of marbles with known number ($N$) of marbles, but the number of successes (white marbles) in the bag ($K$) is unknown. I want to infer about K.

Given a sample from this population, where I see $n$ trials and $k$ successes, how do I infer something about the population $K$?

Ideally, I'd like to construct a prior distribution around $K$ and then use the sample to update it and get a Bayesian posterior credibility interval around $K$ (for a given credibility score $\alpha$), but I am struggling to complete this practically. I have read that the conjugate prior for the hypergeometric is the beta-binomial. I thought maybe there would be an R function or package that could take the prior parameters, and then update with a sample to give me a credibility interval, but haven't been able to find it.

If the Bayesian setting is difficult, perhaps a confidence interval would suffice... Can anyone either point me to some R functions or tutorials, or some resources that could help? Thanks.

EDIT: To add an example, I can do this in the case of the binomial distribution, to infer $p$, given a sample. I can construct a credibility interval with the binom package

k= 15
n= 25
library(binom)
binom.bayes(k, n, conf.level = .95, tol=.005, type="central")
# method  x  n shape1 shape2      mean     lower     upper  sig
#  bayes 15 25   15.5   10.5 0.5961538 0.4057793 0.7725105 0.05


To add a prior, because of the way updates work with the beta-binomial, I can just add counts to the $k$ and $n$ parameters according to the $a$ and $b$ parameters of the prior beta distribution.

In the beta-binomial example, $N$ is infinite, and I'm inferring $p$. What I want to do is take this exact situation and just extend it to the case where $N$ is finite (and known), and infer $K$ (which would be equivalent to inferring $p$). This changes the binomial to the hypergeometric.

I have also been interested in this question, and I initially came up with a solution similar to Phil's, where

$$P(K|k)=\frac{\frac{\binom{K}{k}*\binom{N-K}{n-k}}{\binom{N}{n}}}{\sum_{j=k}^{N-n+k} \frac{\binom{j}{k}*\binom{N-j}{n-k}}{\binom{N}{n}}}$$

Where the probability of K for a given k is the likelihood of K over the sum of the likelihoods of all Ks. This feels right to me, but I haven't been able to prove that it is precisely the correct answer. However, if it is, then the PDF simplifies (after some handy Wolfram Alpha) to:

$$f(K)=\frac{\binom{K}{k}*\binom{N-K}{n-k}}{\binom{N-k}{n-k}* {}_2 F_1(k+1, n-N; k-N, 1)}$$

However, since the combinations in the numerator and denominator can get rather large, I have found it helpful to use the following form for computation:

$$f(K)=\frac{\prod_{j=1}^{K-k} \frac{(K+j)*(N-K-n+k+j)}{j*(N-K+j)}}{{}_2 F_1(k+1, n-N; k-N, 1)}$$

The first and second moments are rather unwieldy, but do have a closed form:

$$E[K]=\frac{(N-n)*(k+1)* {}_2 F_1 (k+2, n-N+1; k-N+1; 1)}{(N-k)* {}_2 F_1 (k+1, n-N; k-N; 1)}$$ $$E[K^2]=\left((2k+1)* {}_2 F_1 (k+2, n-N+1, k-N+1, 1) + \frac{(k+2)*(n-N+1)* {}_2 F_1 (k+3,n-N+2,k-N+2,1)}{k - N + 1}\right)*\frac{(k+1)*(n-N)}{(k - N)* {}_2 F_1 (k+1, n-N, k-N, 1)} + k^2$$

Where $${}_2 F_1$$ is the hypergeometric function.

I have not been able to find a better method for calculating the CDF than simply summing up the terms of the PDF, but if you do that, you should have everything you need for a confidence interval.

Also note that the special cases where k=0 and k=1 are much easier to develp and do not require using the hypergeometric function.

You can solve your problem using the Maximum Likelihood method rather than the Bayesian method. Just calculate the hypergeometric probability of finding k for a large number of values of K using your known values of k, N and n. The value of K that provides the maximum for this probability is the expected value of K.