蓄水池采样

Posted 2021-03-29 12:50:38 ‐ 9 min read

https://en.wikipedia.org/wiki/Reservoir_sampling

简单算法

ReservoirSample(S[1...n], R[1...k])
// fill the reservoir array
for i := 1 to k
R[i] := S[i]

for i := k + 1 to n
j := randomInteger(1, i)
if j <= k
R[j] := S[i]

一个优化算法

Algorithm L 优化了这个算法，通过计算下一个item进入蓄水池之前已经丢弃了多少项。关键点在于这个值满足几何分布，因此可以在常数时间计算

/* S has items to sample, R will contain the result */
ReservoirSample(S[1..n], R[1..k])
// fill the reservoir array
for i = 1 to k
R[i] := S[i]

/* random() generates a uniform (0,1) random number */
W := exp(log(random())/k)

while i <= n
i := i + floor(log(random())/log(1-W)) + 1
if i <= n
/* replace a random item of the reservoir with item i */
R[randomInteger(1,k)] := S[i]  // random index between 1 and k, inclusive
W := W * exp(log(random())/k)

/*
S is a stream of items to sample
S.Current returns current item in stream
S.Next advances stream to next position
min-priority-queue supports:
Count -> number of items in priority queue
Minimum -> returns minimum key value of all items
Extract-Min() -> Remove the item with minimum key
Insert(key, Item) -> Adds item with specified key
*/
ReservoirSample(S[1..?])
H := new min-priority-queue
while S has data
r := random()   // uniformly random between 0 and 1, exclusive
if H.Count < k
H.Insert(r, S.Current)
else
// keep k items with largest associated keys
if r > H.Minimum
H.Extract-Min()
H.Insert(r, S.Current)
S.Next
return items in H

权重蓄水池

Some applications require items' sampling probabilities to be according to weights associated with each item. For example, it might be required to sample queries in a search engine with weight as number of times they were performed so that the sample can be analyzed for overall impact on user experience. Let the weight of item i be , and the sum of all weights be W. There are two ways to interpret weights assigned to each item in the set:[4]

1. In each round, the probability of every unselected item to be selected in that round is proportional to its weight relative to the weights of all unselected items. If X is the current sample, then the probability of an item $i \notin X$ to be selected in the current round is $\frac{w_i}{W - \sum_{j \in X}w_j}$
2. The probability of each item to be included in the random sample is proportional to its relative weight, i.e. $\frac{w_i}{W}$. Note that this interpretation might not be achievable in some cases, e.g., $k = n$.

Algorithm A-Res

/*
S is a stream of items to sample
S.Current returns current item in stream
S.Weight  returns weight of current item in stream
S.Next advances stream to next position
The power operator is represented by ^
min-priority-queue supports:
Count -> number of items in priority queue
Minimum() -> returns minimum key value of all items
Extract-Min() -> Remove the item with minimum key
Insert(key, Item) -> Adds item with specified key
*/
ReservoirSample(S[1..?])
H := new min-priority-queue
while S has data
r := random() ^ (1/S.Weight)   // random() produces a uniformly random number in (0,1)
if H.Count < k
H.Insert(r, S.Current)
else
// keep k items with largest associated keys
if r > H.Minimum
H.Extract-Min()
H.Insert(r, S.Current)
S.Next
return items in H

Algorithms A-ExpJ

/*
S is a stream of items to sample
S.Current returns current item in stream
S.Weight  returns weight of current item in stream
S.Next advances stream to next position
The power operator is represented by ^
min-priority-queue supports:
Count -> number of items in the priority queue
Minimum -> minimum key of any item in the priority queue
Extract-Min() -> Remove the item with minimum key
Insert(Key, Item) -> Adds item with specified key
*/
ReservoirSampleWithJumps(S[1..?])
H := new min-priority-queue
while S has data and H.Count < k
r := random() ^ (1/S.Weight)   // random() produces a uniformly random number in (0,1)
H.Insert(r, S.Current)
S.Next
X := log(random()) / log(H.Minimum) // this is the amount of weight that needs to be jumped over
while S has data
X := X - S.Weight
if X <= 0
t := H.Minimum ^ S.Weight
r := random(t, 1) ^ (1/S.Weight) // random(x, y) produces a uniformly random number in (x, y)

H.Extract-Min()
H.Insert(r, S.Current)

X := log(random()) / log(H.Minimum)
S.Next
return items in H

Algorithm A-Chao

/*
S has items to sample, R will contain the result
S[i].Weight contains weight for each item
*/
WeightedReservoir-Chao(S[1..n], R[1..k])
WSum := 0
// fill the reservoir array
for i := 1 to k
R[i] := S[i]
WSum := WSum + S[i].Weight
for i := k+1 to n
WSum := WSum + S[i].Weight
p := S[i].Weight / WSum // probability for this item
j := random();          // uniformly random between 0 and 1
if j <= p               // select item according to probability
R[randomInteger(1,k)] := S[i]  //uniform selection in reservoir for replacement

Relation to Fisher-Yates shuffle

Fihser-Yates shuffle

Shuffle(S[1...n], R[1...n])
R[1] := S[1]
for i from 2 to n do
j := randomInteger(1, i)
R[i] := R[j]
R[j] := S[i]

ReservoirSample(S[1..n], R[1..k])
R[1] := S[1]
for i from 2 to k do
j := randomInteger(1, i)  // inclusive range
R[i] := R[j]
R[j] := S[i]
for i from k + 1 to n do
j := randomInteger(1, i)  // inclusive range
if (j <= k)
R[j] := S[i]