如何高效地生成一组具有预定义分布的唯一随机数集?

3

我有一个包含某些概率分布的项目地图:

Map<SingleObjectiveItem, Double> itemsDistribution;

给定一个特定的m,我需要从上述分布中生成m个元素的Set

目前为止,我使用的是朴素的方法:

while(mySet.size < m)
   mySet.add(getNextSample(itemsDistribution));
getNextSample(...) 方法根据其概率从分布中获取一个对象。现在,随着 m 的增加,性能严重下降。对于 m = 500itemsDistribution.size() = 1000 个元素,存在太多 thrashing,并且该函数在 while 循环中停留的时间过长。生成 1000 个这样的集合,您将拥有一个爬行应用程序。
是否有更有效的方法来生成具有“预定义”分布的唯一随机数集?大多数集合洗牌技术和类似技术都是均匀随机的。如何解决这个问题?
更新:循环将“至少”调用 getNextSample(...) “1 + 2 + 3 + ... + m = m(m+1)/2” 次。也就是说,在第一次运行时,我们肯定会为集合获得一个样本。第二次迭代,它可能被调用至少两次,以此类推。如果 getNextSample 是顺序的,即通过整个累积分布查找样本,则循环的运行时间复杂度至少为:n*m(m+1)/2,“n”是分布中元素的数量。如果 m = cn; 0<c<=1,则循环至少为 Sigma(n^3)。那也是下限!
如果我们用二分搜索替换顺序搜索,则复杂度至少为 Sigma(log n * n^2)。效率高,但提升可能不会很大。
此外,由于我调用上述循环 k 次,以生成 k 个这样的集合,因此无法从分布中删除元素。这些集合是项目随机“计划”的一部分。因此是一个“项目”集合。

一个元素可以被多次选择吗?如果不行,那么地图中的值的确切形式意义是什么?它不能仅仅是选择元素的概率,因为当我们已经选择了一些元素并且不能再次选择它们时,这些值失去了某些概率属性。最明显的是,它们不再总和为1。此外,选择项目的顺序可能会干扰选择集的总体概率。例如,从{1,2,3}中,先选择1再选择2的概率可能与先选择2再选择1的概率不同 - 你可能希望在这方面保持一致性。 - Gassa
8个回答

3
首先,在二维平面上生成一些随机点。
然后应用你的分布函数。
现在找到所有在分布函数内的点,并选择它们的x坐标,这样你就得到了按照所需分布的随机数,如下图所示。

不确定这会有多快:假设输入元素服从均匀分布,接受一个点的概率是1/n,也就是说,平均而言,我们需要采样m*n个点才能得到大小为m的集合。那可是相当多的。 - meriton
好的 - 使用正态分布并覆盖3个标准差,大约有1/3的点会在图形下方。 - Ebbe M. Pedersen
但其他发行版比那还要糟糕得多。我只是想指出,拒绝取样需要一个相对较小的分布边界框,而并非所有分布都存在这样的边界框。也就是说,您的答案只对某些分布有效。 - meriton

1
如果您不太关心随机性质,那么我会这样做:
  1. 创建伪随机数缓冲区

    double buff[MAX]; // [edit1] 双精度浮点型伪随机数

    • MAX 是缓冲区大小,应足够大...例如1024*128
    • 类型可以是任何类型(float,int,DWORD...)
  2. 填充缓冲区

    你有一组数字范围 x = < x0,x1 > 和概率函数 probability(x) ,由你的概率分布定义,所以做这个:

    for (i=0,x=x0;x<=x1;x+=stepx)
     for (j=0,n=probability(x)*MAX,q=0.1*stepx/n;j<n;j++,i++) // [edit1] 独特的伪随机数
      buff[i]=x+(double(i)*q);                                // [edit1] ...
    

    其中 stepx 是你的项目精度(对于整数类型 = 1) 现在 buff[] 数组具有与你所需相同的分布,但它不是伪随机的。此外,您应添加检查以确保 j 不是 >= MAX 以避免数组溢出,并且在最后,buff[] 的实际大小是 j (由于舍入可能小于 MAX)

  3. 打乱 buff[]

    只需交换几个循环中的 buff[i]buff[j],其中 i 是循环变量,j 是伪随机的 <0-MAX)

  4. 编写你的伪随机函数

    它只需从缓冲区返回数字。第一次调用返回 buff[0],第二次返回 buff[1],依此类推...对于标准生成器,当你到达 buff[] 的末尾时,再次打乱 buff[] 并从 buff[0] 开始。但是,由于您需要唯一的数字,因此不能到达缓冲区的末尾,因此将 MAX 设置足够大以满足您的任务,否则无法保证唯一性。

[笔记]

MAX 应该足够大,以存储您想要的整个分布。如果不够大,则可能会完全缺少低概率项。

[编辑1] - 稍微调整了答案以符合问题需求(由 meriton 指出感谢)

PS。 初始化的复杂度为 O(N),获取数字的复杂度为 O(1)。


这似乎是不正确的:您没有确保返回的数字是唯一的。此外,权重可能不是整数(问题使用双精度表示它们)。 - meriton
@meriton 哦,我之前忽略了唯一的部分...但是数字可以像我写的那样是双倍的。缓冲填充只需要进行一些小的调整即可获得唯一的数字(在填充时添加一些小的精度漂移<stepx/probability(x)),但当然你不能两次使用缓冲区,所以你只能使用 rnd 最多 MAX 次。 - Spektre

1

问题不太可能是您展示的循环:

设 n 为分布的大小,I 为 getNextSample 调用次数。我们有 I = sum_i(C_i),其中 C_i 是在集合大小为 i 时 getNextSample 的调用次数。要找到 E[C_i],观察到 C_i 是具有 λ = 1 - i / n 的 泊松过程 的到达间隔时间,因此具有 指数分布 λ。因此,E[C_i] = 1/λ = 因此 E[C_i] = 1 / (1 - i / n) <= 1 / (1 - m / n)。因此,E[I] < m / (1 - m / n)。

即,抽取一个大小为m = n/2的集合将平均只需要不到2m = n次调用getNextSample。如果这很“慢”且“缓慢”,那可能是因为getNextSample很慢。这实际上并不奇怪,因为分布以不适当的方式传递给方法(因为方法必须必然遍历整个分布以找到随机元素)。
如果m < 0.8n,则以下应该更快速。
class Distribution<T> {
    private double[] cummulativeWeight;
    private T[] item;
    private double totalWeight;

    Distribution(Map<T, Double> probabilityMap) {
        int i = 0;

        cummulativeWeight = new double[probabilityMap.size()];
        item = (T[]) new Object[probabilityMap.size()];

        for (Map.Entry<T, Double> entry : probabilityMap.entrySet()) {
            item[i] = entry.getKey();
            totalWeight += entry.getValue();
            cummulativeWeight[i] = totalWeight;
            i++;
        }
    }

    T randomItem() {
        double weight = Math.random() * totalWeight;
        int index = Arrays.binarySearch(cummulativeWeight, weight);
        if (index < 0) {
            index = -index - 1;
        }
        return item[index];
    }

    Set<T> randomSubset(int size) {
        Set<T> set = new HashSet<>();
        while(set.size() < size) {
            set.add(randomItem());
        }
        return set;
    }
}



public class Test {

    public static void main(String[] args) {
        int max = 1_000_000;
        HashMap<Integer, Double> probabilities = new HashMap<>();
        for (int i = 0; i < max; i++) {
            probabilities.put(i, (double) i);
        }

        Distribution<Integer> d = new Distribution<>(probabilities);
        Set<Integer> set = d.randomSubset(max / 2);
        //System.out.println(set);
    }
}

预期运行时间为 O(m / (1 - m / n) * log n)。在我的电脑上,对于一个大小为 1,000,000 的集合的子集大小为 500,000,计算大约需要 3 秒钟。
正如我们所看到的,当 m 接近 n 时,预期运行时间趋近于无穷大。如果这是一个问题(即 m > 0.9n),则以下更复杂的方法应该效果更好:
Set<T> randomSubset(int size) {
    Set<T> set = new HashSet<>();
    while(set.size() < size) {
        T randomItem = randomItem();
            remove(randomItem); // removes the item from the distribution
            set.add(randomItem);
    }
    return set;
}

为了高效地实现 remove,需要使用不同的分布表示,例如二叉树,每个节点存储其根节点下子树的总权重。
但这相当复杂,因此如果 m 明显小于 n,则不建议采用该方法。

它为什么是1/(1-c)?如果getNextSample(...)运行时间为O(n)(顺序执行,不幸的是),那么循环预计将运行1 + 2 + 3 + ... + m = m(m+1)/2。如果m = cn,那就容易得到O(n^2)。二分查找肯定会使它更有效率。但我认为它不会有太大的差距。因为循环的期望值仍然是m(m+1)/2。我错过了什么吗? - PhD
我为运行时分析添加了一个证明草图。此外,n表示我们从中取样的集合的大小。因此,E [算法执行时间] = O(E [I] * n)= O(m /(1-c)* n)。话虽如此,你认为用二分查找替换线性查找不会显著提高执行时间,这让我感到困惑。 - meriton
哦,它会的。我在质疑加速是否真的显著。我正在尝试将其更改为二分查找,以查看其效果如何。一个快速问题:你为什么假设指数分布?从概念上讲... - PhD
好的,如果n = 1000,n / log(n) ~= 100 ... 我会编辑关于为什么分布是指数的原因。 - meriton
你是对的。问题不仅仅在于循环,还在于数据捕获的方式。我采纳了你的建议并创建了一个类似的类来解决这个问题。它确实大大加快了速度。只是将循环改为二分查找并没有像我猜测的那样有太大帮助。 - PhD

0

你的表现取决于你的getNextSample函数的工作方式。如果在选择下一个项目时必须遍历所有概率,那么速度可能会很慢。

从列表中选择几个唯一的随机项的好方法是先对列表进行洗牌,然后弹出列表中的项。您可以使用给定的分布一次对列表进行洗牌。从那时起,选择您的m项只需弹出列表即可。

这是一个概率洗牌的实现:

List<Item> prob_shuffle(Map<Item, int> dist)
{
    int n = dist.length;
    List<Item> a = dist.keys();
    int psum = 0;
    int i, j;

    for (i in dist) psum += dist[i];

    for (i = 0; i < n; i++) {
        int ip = rand(psum);    // 0 <= ip < psum
        int jp = 0;

        for (j = i; j < n; j++) {
            jp += dist[a[j]];
            if (ip < jp) break;
        }

        psum -= dist[a[j]];

        Item tmp = a[i];
        a[i] = a[j];
        a[j] = tmp;
    }
    return a;
}

这不是Java,而是在C实现后的伪代码,请谨慎参考。其思想是通过不断从未洗牌区域中选择项目来将项目附加到洗牌区域。

在这里,我使用整数概率。(概率不必添加到特定值,只要“越大越好”即可。)您可以使用浮点数,但由于不准确性,当选择一个项目时,您可能会超出数组范围。然后应该使用项目n-1。如果添加了这个安全网,甚至可以有概率为零的项目,总是最后被选中。

可能有一种方法可以加速选择循环,但我真的看不出来。交换使任何预计算都无用。


0
你应该实现自己的随机数生成器(使用 MonteCarlo 方法或任何良好的均匀生成器,如 Mersen Twister),并基于反演方法(这里)。
例如:指数律:生成 [0,1] 中的均匀随机数 u,然后你的指数律随机变量将是:ln(1-u)/(-lambda) 其中 lambda 是指数律参数,ln 是自然对数。
希望对你有所帮助 ;).

0

我认为你有两个问题:

  1. 你的 itemDistribution 不知道你需要一个集合,所以当你正在构建的集合变得很大时,你会选择很多已经在集合中的元素。如果你从集合中开始并删除元素,对于非常小的集合,你将遇到相同的问题。

    你为什么不在选取后从 itemDistribution 中删除元素呢?那么你就不会重复选择同一个元素了吧?

  2. itemDistribution 的数据结构选择让我感到可疑。你希望 getNextSample 操作快速。从值到概率的映射是否强制你为每个 getNextSample 遍历大部分映射?我不擅长统计学,但你能否用另一种方式表示 itemDistribution,比如从概率到映射,或者所有较小概率的总和+元素集合的概率?


0

在表格中累积您的概率

               Probability
Item       Actual  Accumulated
Item1       0.10      0.10
Item2       0.30      0.40
Item3       0.15      0.55
Item4       0.20      0.75
Item5       0.25      1.00

生成一个0.0到1.0之间的随机数,并使用二分查找找到第一个总和大于生成数的项目。这个项目将以所需的概率被选择。

0

Ebbe的方法被称为拒绝抽样

我有时使用一种简单的方法,使用一个反累积分布函数,它是将0到1之间的数字X映射到Y轴上的函数。 然后你只需要生成一个在0到1之间均匀分布的随机数,并将该函数应用于它。 这个函数也被称为“分位函数”。

例如,假设你想生成一个正态分布的随机数。 它的累积分布函数被称为Phi。 它的反函数被称为probit。 有许多方法可以生成正态变量,这只是一个例子。

你可以很容易地构造任何你喜欢的单变量分布的近似累积分布函数,以表格的形式呈现。 然后你可以通过表格查找和插值来反转它。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接