in rand/rand_distr/src/hypergeometric.rs [238:368]
fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> u64 {
use SamplingMethod::*;
let Hypergeometric { n1, n2, k, sign_x, offset_x, sampling_method } = *self;
let x = match sampling_method {
InverseTransform { initial_p: mut p, initial_x: mut x } => {
let mut u = rng.gen::<f64>();
while u > p && x < k as i64 { // the paper erroneously uses `until n < p`, which doesn't make any sense
u -= p;
p *= ((n1 as i64 - x as i64) * (k as i64 - x as i64)) as f64;
p /= ((x as i64 + 1) * (n2 as i64 - k as i64 + 1 + x as i64)) as f64;
x += 1;
}
x
},
RejectionAcceptance { m, a, lambda_l, lambda_r, x_l, x_r, p1, p2, p3 } => {
let distr_region_select = Uniform::new(0.0, p3);
loop {
let (y, v) = loop {
let u = distr_region_select.sample(rng);
let v = rng.gen::<f64>(); // for the accept/reject decision
if u <= p1 {
// Region 1, central bell
let y = (x_l + u).floor();
break (y, v);
} else if u <= p2 {
// Region 2, left exponential tail
let y = (x_l + v.ln() / lambda_l).floor();
if y as i64 >= i64::max(0, k as i64 - n2 as i64) {
let v = v * (u - p1) * lambda_l;
break (y, v);
}
} else {
// Region 3, right exponential tail
let y = (x_r - v.ln() / lambda_r).floor();
if y as u64 <= u64::min(n1, k) {
let v = v * (u - p2) * lambda_r;
break (y, v);
}
}
};
// Step 4: Acceptance/Rejection Comparison
if m < 100.0 || y <= 50.0 {
// Step 4.1: evaluate f(y) via recursive relationship
let mut f = 1.0;
if m < y {
for i in (m as u64 + 1)..=(y as u64) {
f *= (n1 - i + 1) as f64 * (k - i + 1) as f64;
f /= i as f64 * (n2 - k + i) as f64;
}
} else {
for i in (y as u64 + 1)..=(m as u64) {
f *= i as f64 * (n2 - k + i) as f64;
f /= (n1 - i) as f64 * (k - i) as f64;
}
}
if v <= f { break y as i64; }
} else {
// Step 4.2: Squeezing
let y1 = y + 1.0;
let ym = y - m;
let yn = n1 as f64 - y + 1.0;
let yk = k as f64 - y + 1.0;
let nk = n2 as f64 - k as f64 + y1;
let r = -ym / y1;
let s = ym / yn;
let t = ym / yk;
let e = -ym / nk;
let g = yn * yk / (y1 * nk) - 1.0;
let dg = if g < 0.0 {
1.0 + g
} else {
1.0
};
let gu = g * (1.0 + g * (-0.5 + g / 3.0));
let gl = gu - g.powi(4) / (4.0 * dg);
let xm = m + 0.5;
let xn = n1 as f64 - m + 0.5;
let xk = k as f64 - m + 0.5;
let nm = n2 as f64 - k as f64 + xm;
let ub = xm * r * (1.0 + r * (-0.5 + r / 3.0)) +
xn * s * (1.0 + s * (-0.5 + s / 3.0)) +
xk * t * (1.0 + t * (-0.5 + t / 3.0)) +
nm * e * (1.0 + e * (-0.5 + e / 3.0)) +
y * gu - m * gl + 0.0034;
let av = v.ln();
if av > ub { continue; }
let dr = if r < 0.0 {
xm * r.powi(4) / (1.0 + r)
} else {
xm * r.powi(4)
};
let ds = if s < 0.0 {
xn * s.powi(4) / (1.0 + s)
} else {
xn * s.powi(4)
};
let dt = if t < 0.0 {
xk * t.powi(4) / (1.0 + t)
} else {
xk * t.powi(4)
};
let de = if e < 0.0 {
nm * e.powi(4) / (1.0 + e)
} else {
nm * e.powi(4)
};
if av < ub - 0.25*(dr + ds + dt + de) + (y + m)*(gl - gu) - 0.0078 {
break y as i64;
}
// Step 4.3: Final Acceptance/Rejection Test
let av_critical = a -
ln_of_factorial(y) -
ln_of_factorial(n1 as f64 - y) -
ln_of_factorial(k as f64 - y) -
ln_of_factorial((n2 - k) as f64 + y);
if v.ln() <= av_critical {
break y as i64;
}
}
}
}
};
(offset_x + sign_x * x) as u64
}