11

我正在研究一个模拟系统。我很快就会有几个模拟输入值的真实世界分布的实验数据（直方图）。

当模拟运行时，我希望能够产生与测量分布相匹配的随机值。我宁愿这样做而不存储原始直方图。有什么好的方法

  1. 将直方图映射到一组代表分布的参数？
  2. 在运行时生成基于这些参数的值？

编辑：输入数据是几种不同类型事件的事件持续时间。我期望不同的类型会有不同的分布函数。

4

6 回答 6

19

至少有两种选择：

  1. 对直方图进行积分并进行数值反转。
  2. 拒绝

数值积分

来自William R. Gibbs的《现代物理学计算》 ：

人们总是可以对 [函数] 进行数值积分并反转 [ cdf ]，但这通常不是很令人满意，尤其是在pdf快速变化的情况下。

您实际上构建了一个表格，将范围[0-1)转换为目标分布中的适当范围。然后扔出你通常的（高质量）PRNG 并与表格一起翻译。它很麻烦，但清晰、可行且完全通用。

拒绝：

归一化目标直方图，然后

  1. 掷骰子以x随机选择范围内的位置（ ）。
  2. 再次抛出，如果新的随机数小于这个bin中的归一化直方图，就选择这个点。否则转到 (1)。

再次，头脑简单但清晰且有效。分布可能很慢，概率非常低（长尾峰）。

使用这两种方法，如果不需要阶跃函数直方图，您可以使用分段多项式拟合或样条曲线来近似数据以生成平滑曲线——但将其留到以后，因为它可能是过早的优化。

对于特殊情况，可能存在更好的方法。

所有这些都是非常标准的，如果需要更多细节，应该出现在任何数值分析教科书中。

于 2009-01-08T02:27:52.870 回答
2

有关该问题的更多信息将很有用。例如，直方图是什么类型的值？它们是分类的（例如颜色、字母）还是连续的（例如高度、时间）？

如果直方图超过分类数据，我认为可能难以对分布进行参数化，除非类别之间存在许多相关性。

如果直方图超过连续数据，您可能会尝试使用高斯混合来拟合分布。也就是说，尝试使用 $\sum_{i=1}^n w_i N(m_i,v_i)$ 拟合直方图，其中 m_i 和 v_i 是均值和方差。然后，当您想要生成数据时，您首先从 1..n 中以与权重 w_i 成比例的概率对 i 进行采样，然后像从任何高斯中一样对 x ~ n(m_i,v_i) 进行采样。

无论哪种方式，您都可能想了解更多关于混合模型的信息。

于 2009-01-08T03:09:03.413 回答
1

因此，正如@dmckee 所说 ，为了生成给定的概率分布，我想要的是一个分位数函数，它是 累积分布函数的倒数。

问题变成了：生成和存储描述给定连续直方图的分位数函数的最佳方法是什么？我有一种感觉，答案将很大程度上取决于输入的形状——如果它遵循任何类型的模式，应该对最一般的情况进行简化。我会在我去的时候在这里更新。

编辑：

本周我的一次谈话让我想起了这个问题。如果我放弃将直方图描述为等式，而只存储表格，我可以在 O(1) 时间内进行选择吗？事实证明，您可以在不损失任何精度的情况下，以 O(N lgN) 的构建时间为代价。

创建一个包含 N 个项目的数组。对数组进行均匀随机选择将找到概率为 1/N 的项目。对于每个项目，存储实际应选择该项目的命中分数，以及如果该项目未选择则将选择的另一项目的索引。

加权随机采样，C 实现： 

//data structure
typedef struct wrs_data {
  double share; 
  int pair;
  int idx;
} wrs_t;


//sort helper
int wrs_sharecmp(const void* a, const void* b) {
  double delta = ((wrs_t*)a)->share - ((wrs_t*)b)->share;
  return (delta<0) ? -1 : (delta>0);
}


//Initialize the data structure
wrs_t* wrs_create(int* weights, size_t N) {
  wrs_t* data = malloc(sizeof(wrs_t));
  double sum = 0;
  int i;
  for (i=0;i<N;i++) { sum+=weights[i]; }
  for (i=0;i<N;i++) {
    //what percent of the ideal distribution is in this bucket?
    data[i].share = weights[i]/(sum/N); 
    data[i].pair = N;
    data[i].idx = i;
  }
  //sort ascending by size
  qsort(data,N, sizeof(wrs_t),wrs_sharecmp);

  int j=N-1; //the biggest bucket
  for (i=0;i<j;i++) {
    int check = i;
    double excess = 1.0 - data[check].share;
    while (excess>0 && i<j) {
      //If this bucket has less samples than a flat distribution,
      //it will be hit more frequently than it should be.  
      //So send excess hits to a bucket which has too many samples.
      data[check].pair=j; 
      // Account for the fact that the paired bucket will be hit more often,
      data[j].share -= excess;  
      excess = 1.0 - data[j].share;
      // If paired bucket now has excess hits, send to new largest bucket at j-1
      if (excess >= 0) { check=j--;} 
    }
  }
  return data;
}


int wrs_pick(wrs_t* collection, size_t N)
//O(1) weighted random sampling (after preparing the collection).
//Randomly select a bucket, and a percentage.
//If the percentage is greater than that bucket's share of hits, 
// use it's paired bucket.
{
  int idx = rand_in_range(0,N);
  double pct = rand_percent();
  if (pct > collection[idx].share) { idx = collection[idx].pair; }
  return collection[idx].idx;
}

编辑 2：经过一番研究，我发现甚至可以在 O(N) 时间内完成构建。通过仔细跟踪，您无需对数组进行排序即可找到大小箱。 在这里更新了实现

于 2009-01-08T18:17:44.747 回答
0

如果您需要提取大量具有离散点加权分布的样本，请查看类似问题的答案。

但是，如果您需要使用直方图逼近一些连续随机函数，那么您最好的选择可能是 dmckee 的数字积分答案。或者，您可以使用锯齿，并将点存储在左侧，并在两点之间选择一个统一的数字。

于 2009-01-08T02:52:29.110 回答
0

要从直方图（原始或缩减）中进行选择， Walker 的别名方法 既快速又简单。

于 2009-10-20T12:08:17.747 回答
-3

对于正态分布，以下可能会有所帮助：

http://en.wikipedia.org/wiki/Normal_distribution#Generating_values_for_normal_random_variables

于 2009-01-08T02:18:55.167 回答