[Sinkleberg]

Algorithm D compared the remaining population with n / 13. That kept it from switching to Algorithm A even when the remaining sample was dense. This uses the intended 13 * n threshold and passes the remaining population as a view.

This addresses the second suggestion in #998. The public algorithm selector is outside this change.

The boundary and handoff tests passed, along with the full package suite on Julia 1.10.13.

[devmotion]

Thanks! The threshold fix seems correct, I checked and it matches Vitter's pseudocode in both papers:

• Vitter (1984), "Faster Methods for Random Sampling", CACM 27(7):703–718 (PDF). p. 715 says alpha_inverse (≈ 1/α, with α ≈ 0.07) and that threshold is "always equal to alpha_inverse × n". p. 716 has threshold := alpha_inverse × n; while (n > 1) and (threshold < N) … threshold := threshold − alpha_inverse. seqsample_d! is a line-by-line port of that p. 716 version.
• Vitter (1987), "An Efficient Algorithm for Sequential Random Sampling", ACM TOMS 13(1):58–67 (PDF). In Appendix A2 (p. 66): negalphainv := −13; threshold := −negalphainv × n, with "we used α = 1/13". This is where the 13 comes from.

I also checked the sampler numerically with the PR's code. In the pure-D regime (N = 200, n = 3), the inclusion probabilities and the exact law of the first selected index match. In the mixed D→A regime, chi-square tests on inclusion counts are consistent with uniform. With dense samples from a Vector it's also faster, since A is used for the tail and the tail is no longer copied.

A few comments inline. The main one is that the existing statistical test for D no longer exercises D's rejection loop.

[devmotion]

With this change, the statistical test just above (seqsample_d!(3:12, …) with N = 10, n = 5) never enters D's rejection loop, because 13·5 ≥ 10. I instrumented it and got 0 loop iterations in 1000 runs, so it now only tests Algorithm A. The same goes for this test_rng_use call (N = 10, n = 6).

Could we add check_sample_norep cases where N > 13n, so D's loop and the D→A mix are covered? Something like:

# exercise D's rejection loop (N > 13n) and the mixed D -> A regime
for (N, k) in ((100, 3), (40, 3))
    a = zeros(Int, k, n)
    for j = 1:size(a,2)
        seqsample_d!(1:N, view(a,:,j))
    end
    check_sample_norep(a, (1, N), 5.0e-3; ordered=true)
end

Locally this gives an L2 deviation of about 0.002, well within the 5e-3 tolerance. Maybe also add a test_rng_use(seqsample_d!, 1:100, zeros(Int, 3)).

[devmotion]

Existing issue, but this PR makes it reachable. seqsample_d! has no return x, so the return value now depends on the path:

• after the handoff to A, it returns a SubArray of the tail of x;
• on the n == 1 path, it returns a scalar. For example, seqsample_d!(rng, 1:1000, zeros(Int, 3)) returns an Int.

seqsample_a! and seqsample_c! both return x. Could you add return x after this if/else?

While you're here: with an empty x (n = 0), this branch writes x[1] and throws a BoundsError. A and C guard against this with if k > 0. Fixing it is optional, but it's a one-liner.

[devmotion]

Nit: the comment mentions alpha, which is no longer defined. Maybe something like:

alpha_inv = 13 # 1/α with α = 1/13 as in Vitter (1987); switch to algorithm A once N <= n / α

Vitter uses an integer constant here, so 13 instead of 13.0 would also keep threshold an Int. Either is fine. It might also be worth mentioning the fallback to A and the 1987 TOMS reference in the docstring.

[devmotion]

The SamplingTailProbe type is more machinery than the rest of this file uses, and it's coupled to the exact @view(a[i+1:end]) call. Could we test the handoff behaviourally instead? When N <= 13n, D only draws vprime (one randexp) before calling A, so its output and RNG state should match randexp(rng); seqsample_a!(rng, a, x):

@testset "seqsample_d! hands off to seqsample_a! when N <= 13n" begin
    for N in (25, 26, 27), seed in 1:5
        rng_d = MersenneTwister(seed)
        x_d = seqsample_d!(rng_d, 1:N, zeros(Int, 2))
        rng_a = MersenneTwister(seed)
        randexp(rng_a)
        x_a = seqsample_a!(rng_a, 1:N, zeros(Int, 2))
        @test (x_d == x_a && rng_d == rng_a) == (N <= 26)
    end
end

This passes with the PR and fails on master for N = 25 and 26. Together with the check_sample_norep cases for N > 13n suggested above, I think that's enough, and the SamplingTailProbe testset could be dropped.