program define choi_lr_hyperg_prob, rclass
version 14

* This program calculates the probability mass function defined by
* equation 7 of Choi et al. 2015.
* This is the probability of observing y1 successes
* given n1 and n2 patients on treatments 1 and 2, respectively, and
* is conditioned on the marginal number of successes yplus = y1 + y2.

* Input:
*	n1 = patients on treatment 1
*	n2 = patients on treatment 2
*	y1 = successes on treatment 1
*	y2 = successes on treatment 2
*	psi = log odds ratio for success in Rx 1 vs Rx 2

* Output:
*	lnf = log(f), where f = pmf = equation 7 in Choi et al. 2015

* Equation 7 is of the form f = M/C where M and C can both overflow
* for large tables. We manipulate this formula to avoid overflow errors.
* Let k = log(M). Then log(f) = k - log(C) = -log(C*exp(-k)) where
* C*exp(-k) can be written in such a way to avoid overflow errors
* when n1 and/or n2 are large. For details see CalculatingEq7.pdf in
* \\biostat1158a\Users\dupontwd.VANDERBILT\Documents\Dropbox\OFFICE_desktop
* \ChoiFisherTest\Stata Program\ado\withZeros

    args n1 n2 y1 y2 psi
    local yplus = `y1' + `y2'
    local start = max(0,`yplus' - `n2')
    local end = min(`n1', `yplus')

* Calculate k = log(M), where f =M/C in equation 7

    local k = lnfactorial(`n1') - lnfactorial(`y1') - lnfactorial(`n1'-`y1') ///
            + lnfactorial(`n2') - lnfactorial(`yplus'-`y1')                  ///
            - lnfactorial(`n2' - `yplus' +`y1') + `psi'*`y1'

    local emkc = 0
    forvalues u = `start'/`end' {
	local emkc = `emkc' + exp(lnfactorial(`n1') -lnfactorial(`u')       ///
	           -lnfactorial(`n1'-`u') + lnfactorial(`n2')               ///
	           -lnfactorial(`yplus'-`u')-lnfactorial(`n2'-`yplus'+`u')  ///
	           + `psi'*`u' - `k')
	}
* emkc = C*exp(-k) where C is defined in equation 7

    local lnf = -log(`emkc') // = 0 -log(C*exp(-k)) = log(f)

    return scalar lnf = `lnf'

end // **************************************************************