c++ - 根据其概率选择一个矩阵单元

Question

我有一个正实值的2D矩阵，存储如下:

vector<vector<double>> matrix;

每个单元格可以具有等于或大于0的值，并且该值表示选择该单元格的可能性。特别是，例如，值等于3的单元格与值为1的单元格相比，被选择的概率是其三倍。

我需要随机选择矩阵的 N单元格(0 <= N <=单元格总数)，但要根据其被选择的概率。

我怎样才能做到这一点？

该算法应尽可能快。

Answer 1

我描述了两种方法，A和B。

A的工作时间约为N * number of cells，并使用空间O(log number of cells)。 N很小的时候很好。

B在大约(number of cells + N) * O(log number of cells)的时间上工作，并使用空间O(number of cells)。因此，当N很大(甚至是“medium”)但使用更多的内存时，这是很好的，实际上，由于某些原因，在某些情况下它可能会变慢。

方法A:

您需要做的第一件事是规范化条目。 (对我来说，尚不清楚您是否假设它们已被规范化。)这意味着，将所有条目相加并除以总和。 (这部分可能很慢，因此最好假设或要求已发生。)

然后您像这样采样:

选择矩阵的随机[i,j]条目(通过从整数i,j到0的范围内均匀地随机选择n-1)。

在p范围内选择一个均匀随机的实数[0, 1]。

检查是否matrix[i][j] > p。如果是这样，则返回对[i][j]。如果不是，请返回步骤1。

为什么这样做？我们在第3步以任何特定输出结束的概率等于选择 [i][j]的概率(每个条目都相同)乘以 p数足够小的概率。这与 matrix[i][j]值成正比，因此抽样正在选择具有正确比例的每个条目。在第3步中，我们也有可能回到起点，这是否会使事情产生偏差？基本上没有原因是，假设我们随意选择一个数字 k，然后考虑算法的分布，条件是在 k回合之后精确停止。在假设我们停止在 k的第一个回合的前提下，无论我们选择什么值 k，采样的分布都必须由上述参数完全正确。因为如果我们消除 p太小的情况，那么其他可能性的比例都正确。由于分布对于我们可能要限制的 k的每个值都是完美的，并且总体分布(不是对 k限制的条件)是 k每个值的分布的平均值，因此总体分布也很理想。

如果您想严格地分析通常需要的回合数，则可以通过分析任何特定回合我们实际上在步骤3处停止的概率来做到。由于各回合是独立的，因此每个回合都是相同的，并且从统计上来说，这意味着算法的运行时间是泊松分布的。这意味着它紧密地集中在均值附近，我们可以通过知道该概率来确定均值。

考虑到我们选择了任何特定的条目 [i][j]，可以通过考虑我们在步骤3停止的条件概率来确定在步骤3停止的概率。通过条件期望的公式，您可以得出

Pr[ stop at step 3 ] = sum_{i,j} ( 1/(n^2) * Matrix[i,j] )

由于我们假设矩阵是规范化的，因此总和减少为 1/n^2。因此，无论矩阵中的条目是什么，预期的回合数都约为 n^2(即 n^2直到一个恒定因子)。您不能希望做得比我想的要好得多-读取矩阵的所有条目所花费的时间几乎相同，而且很难从分布中进行采样，甚至无法读取所有。

注意:我所描述的是一种正确采样单个元素的方法-要从一个矩阵中获取 N元素，您只需重复 N次即可。

方法B:

基本上，您只想计算一个直方图并从中进行逆采样，以便您知道完全正确的分布。计算直方图是昂贵的，但是一旦有了它，获取样本既便宜又容易。

在C++中，它可能看起来像这样:

// Make histogram
typedef unsigned int uint;
typedef std::pair<uint, uint> upair;
typedef std::map<double, upair> histogram_type;
histogram_type histogram;
double cumulative = 0.0f;
for (uint i = 0; i < Matrix.size(); ++i) {
  for (uint j = 0; j < Matrix[i].size(); ++j) {
    cumulative += Matrix[i][j];
    histogram[cumulative] = std::make_pair(i,j);
  }
}

std::vector<upair> result;
for (uint k = 0; k < N; ++k) {
  // Do a sample (this should never repeat... if it does not find a lower bound you could also assert false quite reasonably since it means something is wrong with rand() implementation)
  while(1) {
    double p = cumulative * rand(); // Or, for best results use std::mt19937 or boost::mt19937 and sample a real in the range [0,1] here.
    histogram_type::iterator it = histogram::lower_bound(p);
    if (it != histogram.end()) {
      result.push_back(it->second);
      break;
    }
  }
}
return result;

在这里，制作直方图的时间类似于 number of cells * O(log number of cells)，因为插入 map 需要时间 O(log n)。您需要一个有序的数据结构，以便以后进行重复采样时获得便宜的查询 N * O(log number of cells)。可能您可以选择一种更专业的数据结构来加快速度，但是我认为改进的空间有限。

编辑:正如@Bob__在注释中指出的那样，在方法(B)中，如果矩阵很大，即使在此行上使用 double类型，也可能由于浮点舍入而出现一些错误:

cumulative += Matrix[i][j];

问题是，如果 cumulative比 Matrix[i][j]大得多，超出了浮点精度可以处理的范围，那么每次执行此语句时，您可能会观察到明显的错误，这些错误会不断累积，从而带来严重的不准确性。

正如他建议的那样，如果发生这种情况，最简单的解决方法是先对 Matrix[i][j]值进行排序。为了安全起见，您甚至可以在常规实现中执行此操作-渐渐地对这些家伙进行排序不会花费比您已经拥有的更多的时间。