如何在不更换的情况下逐步取样？

Question

Python必须my_sample = random.sample(range(100), 10)随机抽样而无需替换[0, 100).

假设我已经对n这些数字进行了采样,现在我想再进行一次采样而不进行替换(不包括之前的任何采样n),如何超高效地进行采样？

更新:从"合理有效"改为"超级有效"(但忽略常数因素)

Answer 1

如果你事先知道你会想要多个样本不重叠,最简单的是做random.shuffle()对list(range(100))(Python的3 -可以跳过list()在Python 2),然后根据需要剥离片.

s = list(range(100))
random.shuffle(s)
first_sample = s[-10:]
del s[-10:]
second_sample = s[-10:]
del s[-10:]
# etc

其他@ Chronial的答案相当有效.

Answer 2

简短的方法

如果采样的数量比人口少得多,只需采样,检查是否已经选择并重复采样.这可能听起来很愚蠢,但是你选择相同的数字会有一个指数衰减的可能性,所以它比O(n)你没有选择一小部分时快得多.

漫长的道路

Python使用一个梅森倍捻机为PRNG,这是好足够.我们可以完全使用其他东西来以可预测的方式生成非重叠数字.

这是秘密:

• 二次残基x² mod p在是2x < p和p是素数时是唯一的.

• 如果你"翻转"残留物,p - (x² % p)在这个时间也是如此p = 3 mod 4,结果将是剩余的空间.

• 这不是一个非常有说服力的数字传播,所以你可以增加功率,添加一些软糖常数,然后分布是相当不错的.

首先,我们需要生成素数:

from itertools import count
from math import ceil
from random import randrange

def modprime_at_least(number):
    if number <= 2:
        return 2

    number = (number // 4 * 4) + 3
    for number in count(number, 4):
        if all(number % factor for factor in range(3, ceil(number ** 0.5)+1, 2)):
            return number

您可能会担心生成素数的成本.对于10⁶元素,这需要十分之一毫秒.运行[None] * 10**6需要更长的时间,因为它只计算一次,这不是一个真正的问题.

此外,该算法不需要素数的精确值; 只需要一个最多是大于输入数的常数因子的东西.这可以通过保存值列表并搜索它们来实现.如果你进行线性扫描,那就是O(log number),如果你进行二分搜索,那就是O(log number of cached primes).事实上,如果你使用疾驰,你可以把它带到O(log log number),这基本上是常数(log log googol = 2).

然后我们实现了生成器

def sample_generator(up_to):
    prime = modprime_at_least(up_to+1)

    # Fudge to make it less predictable
    fudge_power = 2**randrange(7, 11)
    fudge_constant = randrange(prime//2, prime)
    fudge_factor = randrange(prime//2, prime)

    def permute(x):
        permuted = pow(x, fudge_power, prime) 
        return permuted if 2*x <= prime else prime - permuted

    for x in range(prime):
        res = (permute(x) + fudge_constant) % prime
        res = permute((res * fudge_factor) % prime)

        if res < up_to:
            yield res

并检查它是否有效:

set(sample_generator(10000)) ^ set(range(10000))
#>>> set()

现在,关于这个可爱的事情是,如果你忽视的首要考验,这大约是O(?n)这里n是元素的个数,该算法的时间复杂度O(k),这里k是样本sizeit的和O(1)内存使用情况!从技术上讲,这是O(?n + k),但实际上是O(k).

要求:

1. 您不需要经过验证的PRNG.这个PRNG比线性同余生成器(它很受欢迎; Java使用它)要好得多,但它不像Mersenne Twister那样经过验证.

2. 您不首先生成具有不同功能的任何项目.这可以通过数学避免重复,而不是检查.下一节我将展示如何删除此限制.

3. 简短的方法必须是不够的(k必须接近n).如果k只有一半n,请按照我原来的建议.

好处:

1. 极节省内存.这需要不断的记忆...甚至不是O(k)!

2. 生成下一个项目的恒定时间.这实际上是在不断的条款相当快,太:这不是因为快,内置梅森难题,但它的2倍之内.

3. 凉意.

要删除此要求:

您不首先生成具有不同功能的任何项目.这可以通过数学避免重复,而不是检查.

我在时间和空间复杂度方面做了最好的算法,这是我以前的生成器的简单扩展.

这是纲要(n是数字池的长度,k是"外来"键的数量):

初始化时间O(?n); O(log log n)所有合理的投入

由于O(?n)成本原因,这是我的算法中唯一在算法复杂性方面技术上并不完美的因素.实际上,这不会有问题,因为预先计算将其降低到O(log log n)不可估量地接近恒定时间.

如果您以任何固定百分比耗尽可迭代成本,则成本可以免费摊销.

这不是一个实际问题.

摊销O(1)密钥生成时间

显然这无法改进.

最坏的O(k)密钥生成时间

如果你有从外部生成的密钥,只要求它不能是这个生成器已经生成的密钥,那么这些密钥将被称为"外键".外键被假定为完全随机的.因此,任何能够从池中选择项目的功能都可以这样做.

因为可以有任意数量的外键并且它们可以是完全随机的,所以对于完美算法来说最坏的情况是O(k).

最坏的空间复杂性 O(k)

如果假定外键完全独立,则每个外键代表不同的信息项.因此必须存储所有密钥.算法恰好在看到密钥时丢弃密钥,因此内存成本将在生成器的生命周期内清除.

算法

嗯,这是我的算法.它实际上非常简单:

def sample_generator(up_to, previously_chosen=set(), *, prune=True):
    prime = modprime_at_least(up_to+1)

    # Fudge to make it less predictable
    fudge_power = 2**randrange(7, 11)
    fudge_constant = randrange(prime//2, prime)
    fudge_factor = randrange(prime//2, prime)

    def permute(x):
        permuted = pow(x, fudge_power, prime) 
        return permuted if 2*x <= prime else prime - permuted

    for x in range(prime):
        res = (permute(x) + fudge_constant) % prime
        res = permute((res * fudge_factor) % prime)

        if res in previously_chosen:
            if prune:
                previously_chosen.remove(res)

        elif res < up_to:
            yield res

更改就像添加:

if res in previously_chosen:
    previously_chosen.remove(res)

您可以previously_chosen随时通过添加到set您传入的内容来添加.实际上,您也可以从集合中删除以添加回潜在的池,尽管这只有在sample_generator尚未产生或跳过它时才有效.与prune=False.

所以有.很容易看出它满足所有要求,并且很容易看出要求是绝对的.请注意,如果您没有集合,它仍然可以通过将输入转换为集合来满足其最坏情况,尽管它会增加开销.

测试RNG的质量

从统计学的角度来说,我很好奇这个PRNG实际上有多好.

一些快速搜索引导我创建这三个测试,这些测试似乎都显示出良好的效果!

首先,一些随机数:

N = 1000000

my_gen = list(sample_generator(N))

target = list(range(N))
random.shuffle(target)

control = list(range(N))
random.shuffle(control)

这些被"洗牌"10⁶号码表,从0以10?-1使用Mersenne扭曲作为基线,一个使用我们的乐趣捏造PRNG,其他.第三是控制.

这是一个测试,它查看沿线两个随机数之间的平均距离.将差异与对照进行比较:

from collections import Counter

def birthdat_calc(randoms):
    return Counter(abs(r1-r2)//10000 for r1, r2 in zip(randoms, randoms[1:]))

def birthday_compare(randoms_1, randoms_2):
    birthday_1 = sorted(birthdat_calc(randoms_1).items())
    birthday_2 = sorted(birthdat_calc(randoms_2).items())

    return sum(abs(n1 - n2) for (i1, n1), (i2, n2) in zip(birthday_1, birthday_2))

print(birthday_compare(my_gen, target), birthday_compare(control, target))
#>>> 9514 10136

这小于每个的方差.

这是一个测试,它依次获取5个数字并查看元素的顺序.它们应该在所有120个可能的订单之间平均分配.

def permutations_calc(randoms):
    permutations = Counter()        

    for items in zip(*[iter(randoms)]*5):
        sorteditems = sorted(items)
        permutations[tuple(sorteditems.index(item) for item in items)] += 1

    return permutations

def permutations_compare(randoms_1, randoms_2):
    permutations_1 = permutations_calc(randoms_1)
    permutations_2 = permutations_calc(randoms_2)

    keys = sorted(permutations_1.keys() | permutations_2.keys())

    return sum(abs(permutations_1[key] - permutations_2[key]) for key in keys)

print(permutations_compare(my_gen, target), permutations_compare(control, target))
#>>> 5324 5368

这再次小于每个的方差.

这是一个测试,看看"运行"有多长,也就是说.连续增加或减少的部分.

def runs_calc(randoms):
    runs = Counter()

    run = 0
    for item in randoms:
        if run == 0:
            run = 1

        elif run == 1:
            run = 2
            increasing = item > last

        else:
            if (item > last) == increasing:
                run += 1

            else:
                runs[run] += 1
                run = 0

        last = item

    return runs

def runs_compare(randoms_1, randoms_2):
    runs_1 = runs_calc(randoms_1)
    runs_2 = runs_calc(randoms_2)

    keys = sorted(runs_1.keys() | runs_2.keys())

    return sum(abs(runs_1[key] - runs_2[key]) for key in keys)

print(runs_compare(my_gen, target), runs_compare(control, target))
#>>> 1270 975

这里的差异是非常大的,并且在几次执行中,我看起来两者都是均匀的.因此,该测试通过.

我提到了一个线性同余发生器,可能"更有成效".我已经做了一个糟糕的LCG我自己的LCG,看看这是否是一个准确的声明.

LCG,AFAICT,就像普通发电机一样,它们不是循环的.因此,我看过的大多数参考资料,又名.维基百科,仅涵盖了定义期间的内容,而不是如何在特定时期内制定强大的LCG.这可能会影响结果.

开始:

from operator import mul
from functools import reduce

# Credit http://stackoverflow.com/a/16996439/1763356
# Meta: Also Tobias Kienzler seems to have credit for my
#       edit to the post, what's up with that?
def factors(n):
    d = 2
    while d**2 <= n:
        while not n % d:
            yield d
            n //= d
        d += 1
    if n > 1:
       yield n

def sample_generator3(up_to):
    for modulier in count(up_to):
        modulier_factors = set(factors(modulier))
        multiplier = reduce(mul, modulier_factors)
        if not modulier % 4:
            multiplier *= 2

        if multiplier < modulier - 1:
            multiplier += 1
            break

    x = randrange(0, up_to)

    fudge_constant = random.randrange(0, modulier)
    for modfact in modulier_factors:
        while not fudge_constant % modfact:
            fudge_constant //= modfact

    for _ in range(modulier):
        if x < up_to:
            yield x

        x = (x * multiplier + fudge_constant) % modulier

我们不再检查质数,但我们确实需要做一些奇怪的事情.

• modulier ? up_to > multiplier, fudge_constant > 0
• a - 1必须被modulier...中的每个因素整除
• ......而fudge_constant必须互质与modulier

请注意,这些不是LCG的规则,而是具有完整周期的LCG,显然等于modulier.

我是这样做的:

• modulier至少尝试一下up_to,在条件满足时停止
  • 制定一系列因素,
  • Let multiplier的产品 with duplicates removed
  • If multiplier不小于modulier,继续下一个modulier
  • 让fudge_constant是数字小于modulier,随机选择的
  • 拆下的因素fudge_constant是在

This is not a very good way of generating it, but I don't see why it would ever impinge the quality of the numbers, aside from the fact that low fudge_constants并且multiplier比这些可能造成的完美发生器更常见的事实.

无论如何,结果令人震惊:

print(birthday_compare(lcg, target), birthday_compare(control, target))
#>>> 22532 10650

print(permutations_compare(lcg, target), permutations_compare(control, target))
#>>> 17968 5820

print(runs_compare(lcg, target), runs_compare(control, target))
#>>> 8320 662

总之,我的RNG是好的,线性同余生成器不是.考虑到Java使用线性同余生成器(尽管它只使用较低的位),我希望我的版本绰绰有余.

Answer 3

好的,我们走了.这应该是最快的非概率算法.它具有运行时O(k?log²(s) + f?log(f)) ? O(k?log²(f+k) + f?log(f)))和空间O(k+f).f是禁号的数量,s是禁止数字最长条纹的长度.对此的期望更加复杂,但显然受到约束f.如果你认为s^log?(s)比我的大f或只是不高兴的事实,s再次概率,您可以在日志部改变为分搜索forbidden[pos:]得到O(k?log(f+k) + f?log(f)).

这里的实际实现是O(k?(k+f)+f?log(f)),在列表中插入forbid是O(n).通过用blist sortedlist替换该列表,可以很容易地解决这个问题.

我还添加了一些评论,因为这个算法非常复杂.该lin部件与log部件相同,但需要s而不是log²(s)时间.

import bisect
import random

def sample(k, end, forbid):
    forbidden = sorted(forbid)
    out = []
    # remove the last block from forbidden if it touches end
    for end in reversed(xrange(end+1)):
        if len(forbidden) > 0 and forbidden[-1] == end:
            del forbidden[-1]
        else:
            break

    for i in xrange(k):
        v = random.randrange(end - len(forbidden) + 1)
        # increase v by the number of values < v
        pos = bisect.bisect(forbidden, v)
        v += pos
        # this number might also be already taken, find the
        # first free spot
        ##### linear
        #while pos < len(forbidden) and forbidden[pos] <=v:
        #    pos += 1
        #    v += 1
        ##### log
        while pos < len(forbidden) and forbidden[pos] <= v:
            step = 2
            # when this is finished, we know that:
            # • forbidden[pos + step/2] <= v + step/2
            # • forbidden[pos + step]   >  v + step
            # so repeat until (checked by outer loop):
            #   forbidden[pos + step/2] == v + step/2
            while (pos + step <= len(forbidden)) and \
                  (forbidden[pos + step - 1] <= v + step - 1):
                step = step << 1
            pos += step >> 1
            v += step >> 1

        if v == end:
            end -= 1
        else:
            bisect.insort(forbidden, v)
        out.append(v)
    return out

现在到比作"黑客"(在蟒蛇的默认实现),其Veedrac提出的,它具有空间O(f+k)和(n/(n-(f+k))是"猜测"的预计数)时间:

我只是策划了这对k=10和合理的大n=10000(它只是变得更极端的更大的n).而且我必须说:我只是实施了这个,因为它似乎是一个有趣的挑战,但即使我对这是多么极端感到惊讶:

让我们放大一下,看看发生了什么:

是的 - 对于您生成的第9998个数字,猜测更快.请注意,正如您在第一个图中所看到的那样,即使是我的单行也可能更快f/n(但仍然有相当可怕的大空间要求n).

为了推动这一点:你在这里花费时间的唯一事情是生成集合,因为这是fVeedrac方法中的因素.

所以我希望我的时间不会浪费,我设法说服你,Veedrac的方法就是这样.我可以理解为什么概率部分会给你带来麻烦,但也许可以想到这样的事实:hashmaps(= python dicts)和大量其他算法使用类似的方法,它们似乎做得很好.

你可能会害怕重复次数的变化.如上面所指出的,这遵循几何分布与p=n-f/n.所以标准偏差(=你应该"期望"的结果偏离预期的平均值)是

这与mean(?f?n < ?n² = n)基本相同.

****编辑**:
我刚才意识到这s也是n/(n-(f+k)).因此,我的算法的运行时更精确O(k?log²(n/(n-(f+k))) + f?log(f)).这是很好的,因为给出上面的图表,它证明了我的直觉,这比它快得多O(k?log(f+k) + f?log(f)).但请放心,这也不会改变上述结果,因为它f?log(f)是运行时中绝对主导的部分.

Answer 4

来自OP的读者注意:请考虑查看最初接受的答案以了解逻辑,然后了解这个答案.

Aaaaaand为了完整起见:这是死灵法师的答案的概念,但是进行了调整,因此需要一个禁止数字列表作为输入.这与我之前的答案中的代码相同,但我们forbid在生成数字之前构建了一个状态.

• 这是时间O(f+k)和记忆O(f+k).显然,这是最快的东西可能不朝的格式要求forbid(整理/套).我认为这使得它在某种程度上成为赢家^^.
• 如果forbid是一组,反复猜测的方法是更快O(k?n/(n-(f+k))),这是非常接近O(k)的f+k不是很接近n.
• 如果forbid排序,我的荒谬算法更快:
  

import random
def sample_gen(n, forbid):
    state = dict()
    track = dict()
    for (i, o) in enumerate(forbid):
        x = track.get(o, o)
        t = state.get(n-i-1, n-i-1)
        state[x] = t
        track[t] = x
        state.pop(n-i-1, None)
        track.pop(o, None)
    del track
    for remaining in xrange(n-len(forbid), 0, -1):
        i = random.randrange(remaining)
        yield state.get(i, i)
        state[i] = state.get(remaining - 1, remaining - 1)
        state.pop(remaining - 1, None)

用法:

gen = sample_gen(10, [1, 2, 4, 8])
print gen.next()
print gen.next()
print gen.next()
print gen.next()

Answer 5

这是一种不明确构建差异集的方法.但它确实使用了@ Veedrac的"接受/拒绝"逻辑形式.如果您不愿意随意改变基本序列,我担心这是不可避免的:

def sample(n, base, forbidden):
    # base is iterable, forbidden is a set.
    # Every element of forbidden must be in base.
    # forbidden is updated.
    from random import random
    nusable = len(base) - len(forbidden)
    assert nusable >= n
    result = []
    if n == 0:
        return result
    for elt in base:
        if elt in forbidden:
            continue
        if nusable * random() < n:
            result.append(elt)
            forbidden.add(elt)
            n -= 1
            if n == 0:
                return result
        nusable -= 1
    assert False, "oops!"

这是一个小司机:

base = list(range(100))
forbidden = set()
for i in range(10):
    print sample(10, base, forbidden)

Answer 6

好的,最后一次尝试;-)以改变基本序列为代价,这不需要额外的空间,并且需要与n每次sample(n)调用成比例的时间:

class Sampler(object):
    def __init__(self, base):
        self.base = base
        self.navail = len(base)
    def sample(self, n):
        from random import randrange
        if n < 0:
            raise ValueError("n must be >= 0")
        if n > self.navail:
            raise ValueError("fewer than %s unused remain" % n)
        base = self.base
        for _ in range(n):
            i = randrange(self.navail)
            self.navail -= 1
            base[i], base[self.navail] = base[self.navail], base[i]
        return base[self.navail : self.navail + n]

小司机:

s = Sampler(list(range(100)))
for i in range(9):
    print s.sample(10)
    print s.sample(1)
print s.sample(1)

实际上,这实现了在选择元素random.shuffle()之后的可恢复的暂停n. base没有被破坏,但是被置换了.

Answer 7

你可以实现一个洗牌生成器,基于维基百科的"Fisher - Yates shuffle#Modern method"

def shuffle_gen(src):
    """ yields random items from base without repetition. Clobbers `src`. """
    for remaining in xrange(len(src), 0, -1):
        i = random.randrange(remaining)
        yield src[i]
        src[i] = src[remaining - 1]

然后可以使用itertools.islice以下方法切片:

>>> import itertools
>>> sampler = shuffle_gen(range(100))
>>> sample1 = list(itertools.islice(sampler, 10))
>>> sample1
[37, 1, 51, 82, 83, 12, 31, 56, 15, 92]
>>> sample2 = list(itertools.islice(sampler, 80))
>>> sample2
[79, 66, 65, 23, 63, 14, 30, 38, 41, 3, 47, 42, 22, 11, 91, 16, 58, 20, 96, 32, 76, 55, 59, 53, 94, 88, 21, 9, 90, 75, 74, 29, 48, 28, 0, 89, 46, 70, 60, 73, 71, 72, 93, 24, 34, 26, 99, 97, 39, 17, 86, 52, 44, 40, 49, 77, 8, 61, 18, 87, 13, 78, 62, 25, 36, 7, 84, 2, 6, 81, 10, 80, 45, 57, 5, 64, 33, 95, 43, 68]
>>> sample3 = list(itertools.islice(sampler, 20))
>>> sample3
[85, 19, 54, 27, 35, 4, 98, 50, 67, 69]

Answer 8

这是我的Knuth shuffle版本,由蒂姆·彼得斯首先发布,由埃里克美化,然后由死灵法师进行了很好的空间优化.

这是基于Eric的版本,因为我确实发现他的代码非常漂亮:).

import random
def shuffle_gen(n):
    # this is used like a range(n) list, but we don’t store
    # those entries where state[i] = i.
    state = dict()
    for remaining in xrange(n, 0, -1):
        i = random.randrange(remaining)
        yield state.get(i,i)
        state[i] = state.get(remaining - 1,remaining - 1)
        # Cleanup – we don’t need this information anymore
        state.pop(remaining - 1, None)

用法:

out = []
gen = shuffle_gen(100)
for n in range(100):
    out.append(gen.next())
print out, len(set(out))

Answer 9

相当快的单线（O(n + m), n=range,m=old samplesize）：

next_sample = random.sample(set(range(100)).difference(my_sample), 10)

Answer 10

编辑:请参阅@TimPeters和@Chronial下面的更清洁版本.一个小编辑将其推到了顶部.

以下是我认为最有效的增量采样解决方案.由调用者维护的状态包括准备好由增量采样器使用的字典,以及该范围内剩余的数字计数,而不是先前采样的数字的列表.

以下是演示实现.与其他解决方案相比:

• 没有循环(没有标准的Python/Veedrac黑客; Python impl和Veedrac之间的共享信用)
• 时间的复杂性是 O(log(number_previously_sampled))
• 空间复杂性是 O(number_previously_sampled)

码:

import random

def remove (i, n, state):
  if i == n - 1:
    if i in state:
      t = state[i]
      del state[i]
      return t
    else:
      return i
  else:
    if i in state:
      t = state[i]
      if n - 1 in state:
        state[i] = state[n - 1]
        del state[n - 1]
      else:
        state[i] = n - 1
      return t
    else:
      if n - 1 in state:
        state[i] = state[n - 1]
        del state[n - 1]
      else:
        state[i] = n - 1
      return i

s = dict()
for n in range(100, 0, -1):
  print remove(random.randrange(n), n, s)

Answer 11

这是@ necromancer的酷解决方案的重写版本.将它包装在一个类中以使其更容易正确使用,并使用更多的dict方法来切割代码行.

from random import randrange

class Sampler:
    def __init__(self, n):
        self.n = n # number remaining from original range(n)
        # i is a key iff i < n and i already returned;
        # in that case, state[i] is a value to return
        # instead of i.
        self.state = dict()

    def get(self):
        n = self.n
        if n <= 0:
            raise ValueError("range exhausted")
        result = i = randrange(n)
        state = self.state
        # Most of the fiddling here is just to get
        # rid of state[n-1] (if it exists).  It's a
        # space optimization.
        if i == n - 1:
            if i in state:
                result = state.pop(i)
        elif i in state:
            result = state[i]
            if n - 1 in state:
                state[i] = state.pop(n - 1)
            else:
                state[i] = n - 1
        elif n - 1 in state:
            state[i] = state.pop(n - 1)
        else:
            state[i] = n - 1
        self.n = n-1
        return result

这是一个基本的驱动程序:

s = Sampler(100)
allx = [s.get() for _ in range(100)]
assert sorted(allx) == list(range(100))

from collections import Counter
c = Counter()
for i in range(6000):
    s = Sampler(3)
    one = tuple(s.get() for _ in range(3))
    c[one] += 1
for k, v in sorted(c.items()):
    print(k, v)

和样本输出:

(0, 1, 2) 1001
(0, 2, 1) 991
(1, 0, 2) 995
(1, 2, 0) 1044
(2, 0, 1) 950
(2, 1, 0) 1019

通过眼球,这种分布很好(如果你持怀疑态度,那就进行卡方检验).这里的一些解决方案没有给出每个置换具有相等的概率(即使它们以相等的概率返回n的每个k子集),因此random.sample()在这方面不同.