3
$\begingroup$

This is my first time posting so I apologize in advance for any errors! I am struggling to understand the probability mass function (PMF) of the Fisher's noncentral hypergeometric distribution, which I am hoping to implement as a custom distribution in Stan. I believe I understand the theory behind the distribution and have read my fair share of urn examples at this point, but the notation in the PMF continues to confuse me. The PMF for the (central) hypergeometric distribution makes sense and I have successfully implemented it in R, obtaining the same results as the built-in function dhyper:

N = 1000  # total number of balls in urn
mx = 300  # number of red balls in urn
my = N - mx  # number of white balls in urn
n = 250  # total number of balls sampled
x = 150  # number of red balls sampled
y = n - x  # number of white balls sampled 

# R's hypergeometric distribution
dhyper(x, mx, my, n)

# my hypergeometric distribution
log_p <- lchoose(mx, x) + lchoose(my, y) - lchoose(N, n)
exp(log_p)

1 3.986613e-31

1 3.986613e-31

In the noncentral case, however, the denominator becomes a more complicated sum of binomial coefficients with indexing that I struggle to interpret. Additionally, superscripts start showing up on the noncentrality term, $\omega$, and it is not clear to me whether those are meant to be interpreted as exponents or are indices for the odds ratios with respect to different outcomes, or perhaps are just the individual odds for different outcomes, not ratios at all? From the Wikipedia entry:

enter image description here

I recognize that there are somewhat complex constraints to ensure everything stays positive, but I would love to find a clean expression of the (log)probability for a typical case in which the total number of samples is less than the population of either group (n < mx, n < my), if such a thing is possible.

$\endgroup$

1 Answer 1

2
$\begingroup$

The pmf is proportional to $${{m_X}\choose x} {{m_Y}\choose y} w_1^xw_2^y$$ where $w_x$ and $w_y$ are weights on red and black balls. The two binomial coefficients give the number of ways to get $(x,y)$ as the answer, $w_1^xw_2^y$ indicates how the non-centrality makes draws with counts $(x,y)$ more or less likely than draws with other proportions.

The $w$s only matter relatively, so it will be tidier to replace them by the ratio $w_1/w_2=\omega$ giving $${{m_X}\choose x} {{m_Y}\choose y} \omega^xw_2^n$$

Call that expression $Q(x,y)$. It's correct up to proportionality, so the actual pmf is

$$p(x,y) = \frac{Q(x,y)}{\sum_{\textrm{all the possible distinct $(x',y')$}} Q(x',y')}$$

One good thing happens here: the $w_2^n$ cancels out. We now have $$p(x,y) = \frac{{{m_X}\choose x} {{m_Y}\choose y} \omega^x}{\sum_{\textrm{all the possible distinct $(x',y')$}} {{m_X}\choose x'} {{m_Y}\choose y'} \omega^{x'}}$$

So the only thing we need is to decide what all possible distinct $(x',y')$ there are. The easy case should be large $m_X$ and $m_Y$, where we don't run out of red or black balls. In that scenario, $x'$ can be anything in $0,1,\dots,n$ and $y'=n-x'$ is then $n,n-1,\dots,0$

$$p(x,y) = \frac{{{m_X}\choose x} {{m_Y}\choose y} \omega^x}{\sum_{x'=0}^n {{m_X}\choose x'} {{m_Y}\choose {n-x'}} \omega^{x'}}$$

And in general it's the same, except you need to choose the starting and ending points for $x'$ so that no cells of the table end up negative.

$\endgroup$
2
  • $\begingroup$ Hi Thomas, thanks for the explanation. This helps a lot. Just to clarify, the $x$ and $y$ in $w_{1}^{x}w_{2}^{y}$ are exponents, not indices, right? Same with $\omega^{x}$ and $\omega^{x'}$? Also, assuming the simple case where $m_{X}$ and $m_{Y}$ are much larger than the number of samples $n$, is there a trick to calculating the sum in the denominator? $\endgroup$
    – dholcomb
    Commented Oct 1, 2020 at 16:49
  • 1
    $\begingroup$ Yes. the lowercase $x$ and $y$ are always the variables, and those are exponents. I don't know about algorithms that avoid the sum in the denominator but there at least should be a simple update formula to go from $(x',y')$ to $(x'+1,y'-1)$. And you could maybe start at the mean and work outwards, stopping when the summands get small, rather than starting at $(0,n)$ $\endgroup$ Commented Oct 1, 2020 at 20:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.