| 117 | */ |
| 118 | |
| 119 | static int64_t hypergeometric_hrua(bitgen_t *bitgen_state, |
| 120 | int64_t good, int64_t bad, int64_t sample) |
| 121 | { |
| 122 | int64_t mingoodbad, maxgoodbad, popsize; |
| 123 | int64_t computed_sample; |
| 124 | double p, q; |
| 125 | double mu, var; |
| 126 | double a, c, b, h, g; |
| 127 | int64_t m, K; |
| 128 | |
| 129 | popsize = good + bad; |
| 130 | computed_sample = MIN(sample, popsize - sample); |
| 131 | mingoodbad = MIN(good, bad); |
| 132 | maxgoodbad = MAX(good, bad); |
| 133 | |
| 134 | /* |
| 135 | * Variables that do not match Stadlober (1989) |
| 136 | * Here Stadlober |
| 137 | * ---------------- --------- |
| 138 | * mingoodbad M |
| 139 | * popsize N |
| 140 | * computed_sample n |
| 141 | */ |
| 142 | |
| 143 | p = ((double) mingoodbad) / popsize; |
| 144 | q = ((double) maxgoodbad) / popsize; |
| 145 | |
| 146 | // mu is the mean of the distribution. |
| 147 | mu = computed_sample * p; |
| 148 | |
| 149 | a = mu + 0.5; |
| 150 | |
| 151 | // var is the variance of the distribution. |
| 152 | var = ((double)(popsize - computed_sample) * |
| 153 | computed_sample * p * q / (popsize - 1)); |
| 154 | |
| 155 | c = sqrt(var + 0.5); |
| 156 | |
| 157 | /* |
| 158 | * h is 2*s_hat (See Stadlober's thesis (1989), Eq. (5.17); or |
| 159 | * Stadlober (1990), Eq. 8). s_hat is the scale of the "table mountain" |
| 160 | * function that dominates the scaled hypergeometric PMF ("scaled" means |
| 161 | * normalized to have a maximum value of 1). |
| 162 | */ |
| 163 | h = D1*c + D2; |
| 164 | |
| 165 | m = (int64_t) floor((double)(computed_sample + 1) * (mingoodbad + 1) / |
| 166 | (popsize + 2)); |
| 167 | |
| 168 | g = (logfactorial(m) + |
| 169 | logfactorial(mingoodbad - m) + |
| 170 | logfactorial(computed_sample - m) + |
| 171 | logfactorial(maxgoodbad - computed_sample + m)); |
| 172 | |
| 173 | /* |
| 174 | * b is the upper bound for random samples: |
| 175 | * ... min(computed_sample, mingoodbad) + 1 is the length of the support. |
| 176 | * ... floor(a + 16*c) is 16 standard deviations beyond the mean. |
no test coverage detected