R Project Sprint 2023 - Enhancing sample.int for unequal probability sampling

[giscus[bot]]

R Project Sprint 2023 - Enhancing sample.int for unequal probability sampling

https://contributor.r-project.org/r-project-sprint-2023/projects/enhance-sample/

[tslumley]

In some very limited experiments, sampling::UPmaxentropy seems to be numerically unstable for population sizes as low as 10,000 (getting NaN values). Less critically, it also takes Nn memory to sample n observations out of N. I'm currently leaning towards Tillé's elimination algorithm (sampling::UPtille, perhaps translated into C) but I need to look at it more.

I think the documentation should indicate that Poisson sampling (ie, random sample size but correct probabilities) and unequal probability sampling with replacement are fast and will be suitable for some applications where without-replacement sampling is slow. It might even be worth having something like method=c("sequential","elimination","poisson") to make Poisson sampling an explicit option.

Finally (for now): unequal probability without replacement differs from all the other methods in that the probabilities must sum to the desired sample size. What should we do if they don't sum to (nearly) an integer? Should there be some sort of rescaling option to allow the user to specify a sample size that isn't the sum of the probabilities (probably not)?

[tslumley]

Profiling sampling::UPtille shows it spends most of its time in inclusionprobabilities. Converting this to C speeds up the function quite a bit: https://gist.github.com/tslumley/d3882e8cc50c610f45052431e73b783a

[mmaechler]

Finally (for now): unequal probability without replacement differs from all the other methods in that the probabilities must sum to the desired sample size. ...

My take on this was that this is their (sampling package) convention.
I thought we should rather keep R's sample() / sample.int() "rule" and arguments (n, size, prob)
Their convention allows them to work with 2 (or even 1!) arguments instead of 3 , and while that may be convenient in their context,
I strongly believe R should keep the meaning of prob to mean that sum(prob) == 1.
From their pik we would extract sum(pik) and length(pik) and translate to our (n, size, prob)
... or is there a flaw in my thinking?

[tslumley]

No, I think sum(prob) really does have to equal size, because you can't just rescale probabilities like that in without-replacement sampling -- that, in a sense, is the whole reason it's hard.

If prob is c(0.7, 0.2, 0.1) and you want two elements, you can't just double the probabilities and use pik=c(1.4, 0.4, 0.2) because those aren't probabilities. Less obviously, if you set pik[i]=2*prob[i] that doesn't ever get you double the chance of sampling i that you'd get with pik[i]=prob[i]. Straightforward rescaling does work with-replacement and it does work for Poisson sampling, but it doesn't work here -- and if it did, the current sequential algorithm would be ok.

[tslumley]

We could of course ask the user to do prob=pik/sum(pik) and then do pik=size*prob internally, but I think that is actually more confusing, since we don't require sum(prob)=1 even in the current interface.

[dickoa]

It's true that in its current implementation, the sample parameter prob is just a measure of size. The only condition they need to satisfy is to be non-negative and not all zero. An easier way to easily integrates pik for the new algorithms is to run something similar to sampling::inclusionprobabilities(prob, size), we can have a helper function similar to the one you build @tslumley for sample_tille.

[dickoa]

I was exploring existing algorithms with good enough properties, and I think the Brewer procedure is a good candidate. Here's a C translation of sampling::UPbrewer

https://gist.github.com/dickoa/eed351acfcdf7383ee900491efcd1ff6

It's 8~9 times faster than UPbrewer on average on some benchmarks, and sampling::UPbrewer was already scaling better than the majority of the UP algorithms.

[tslumley]

As an additional note: the sampling package functions returns a vector of inclusion indicators, but sample.int returns a vector of selected values. We need to select the actual values and permute them randomly -- ie the moral equivalent of

sample(seq_len(n) %in% UPbrewer(prob))

[tslumley]

Actually, it's not at all clear what the sampling order even means for proper PPS sampling (or Poisson sampling)

[tslumley]

R code and Rd in https://gist.github.com/tslumley/16bdc33a5ccfc05d0b8e5d690cb397fd