4

仔细阅读了上一个问题 Random numbers that add to 100: Matlab

我正在努力解决一个类似但稍微复杂的问题。

我想创建一个总和为 1 的 n 个元素的数组，但是我想要一个额外的约束，即每个元素的最小增量（或者如果您喜欢有效数字的数量）是固定的。

例如，如果我想要 10 个总和为 1 的数字而没有任何约束，则以下工作完美：

 num_stocks=10;
 num_simulations=100000;
 temp = [zeros(num_simulations,1),sort(rand(num_simulations,num_stocks-1),2),ones(num_simulations,1)];
 weights = diff(temp,[],2);

我愚蠢地认为通过缩放这个我可以添加如下约束

 num_stocks=10;
 min_increment=0.001;
 num_simulations=100000;
 scaling=1/min_increment;

 temp2 = [zeros(num_simulations,1),sort(round(rand(num_simulations,num_stocks-1)*scaling)/scaling,2),ones(num_simulations,1)];
 weights2 = diff(temp2,[],2);

然而，尽管这适用于小值 n 和小增量值，如果例如 n=1,000 并且增量为 0.1%，那么在大量试验中，第一个和最后一个数字的平均值始终低于 0.1%。

我确信对此有一个合乎逻辑的解释/解决方案，但我一直在努力寻找并找到它，并且想知道是否有人会这么好心地指出我正确的方向。将问题置于上下文中，创建随机股票投资组合（因此总和为 1）。

提前致谢

感谢您到目前为止的回复，只是为了澄清（因为我认为我最初的问题可能措辞不当），权重具有 0.1% 的固定增量，因此 0%、0.1%、0.2% 等。

我最初确实尝试过使用整数

 num_stocks=1000;
 min_increment=0.001;
 num_simulations=100000;
 scaling=1/min_increment;

 temp = [zeros(num_simulations,1),sort(randi([0 scaling],num_simulations,num_stocks-1),2),ones(num_simulations,1)*scaling];
 weights = (diff(temp,[],2)/scaling);
 test=mean(weights);

但更糟糕的是，第一个和最后一个权重的平均值远低于 0.1% .....

编辑以反映弗洛里斯的出色回答并澄清

我用来解决这个问题的原始代码（在找到这个论坛之前）是

function x = monkey_weights_original(simulations,stocks)
stockmatrix=1:stocks;
base_weight=1/stocks;
r=randi(stocks,stocks,simulations);
x=histc(r,stockmatrix)*base_weight;
end

这运行得非常快，考虑到我想运行总共 10,000,000 次模拟，这很重要，对 1,000 只股票进行 10,000 次模拟只需 2 秒多一点，而我正在使用并行工具箱在 8 核机器上运行整个代码。

它还准确地给出了我在均值方面寻找的分布，我认为获得 1 只股票 100% 的投资组合的可能性与获得每只股票 0.1% 的投资组合的可能性一样（尽管我很高兴得到纠正）。

我的问题是，虽然它适用于 1,000 只股票和 0.1% 的增量，我猜它适用于 100 只股票和 1% 的增量，但随着股票数量的减少，每个选择变成一个非常大的百分比（在极端情况下）拥有 2 只股票，您将始终获得 50/50 的投资组合）。

实际上，我认为这个解决方案就像 Floris 建议的二项式解决方案（但更有限）

但是我的问题出现了，因为我想让我的方法更灵活，并且有可能说 3 只股票和 1% 的增量，我当前的代码无法正确处理，因此我是如何偶然发现关于 stackoverflow 的原始问题的

Floris 的递归方法将得到正确的答案，但考虑到问题的规模，速度将是一个主要问题。

原始研究的一个例子是here

http://www.huffingtonpost.com/2013/04/05/monkeys-stocks-study_n_3021285.html

我目前正在努力扩展它，使其在投资组合权重和指数中的股票数量方面具有更大的灵活性，但我的编程和概率论能力似乎是一个限制因素......

4

3 回答 3

2

我可以看到的一个问题是您的公式允许数字为零 - 当舍入操作导致两个连续数字在排序后相同时。不确定您是否认为这是一个问题 - 但我建议您考虑一下（这意味着您的模型投资组合中的股票少于 N 支，因为其中一只股票的贡献为零）。

另一件要注意的是，在你的分布中得到极值的概率是你希望它们的一半：如果你有从 0 到 1000 的均匀分布的数字，你round他们，四舍五入的数字0在间隔[0 0.5>；那些圆形1来自[0.5 1.5>- 两倍大。最后一个数字（四舍五入到1000）再次来自较小的区间：[999.5 1000]。因此，您不会像您想象的那样经常获得第一个和最后一个数字。如果不是round你使用floor我认为你会得到你期望的答案。

编辑

我对此进行了更多考虑，并想出了一个缓慢但（我认为）准确的方法来做到这一点。基本思想是这样的：

  1. 用整数来思考；不是以 0.001 的步长划分区间 0 - 1，而是以整数步长划分区间 0 - 1000
  2. 如果我们尝试将 N 分成 m 个区间，那么一个 step 的平均大小应该是 N/m；但作为整数，我们希望区间是二项式分布的
  3. 这提出了一种算法，在该算法中，我们选择第一个区间作为具有均值的二项式分布变量(N/m)——调用第一个值v1；然后将剩余的间隔N - v1分成m-1步骤；我们可以递归地这样做。

以下代码实现了这一点：

% random integers adding up to a definite sum
function r = randomInt(n, limit)
% returns an array of n random integers
% whose sum is limit
% calls itself recursively; slow but accurate
if n>1
    v = binomialRandom(limit, 1 / n);
    r = [v randomInt(n-1, limit - v)];
else
    r = limit;
end

function b = binomialRandom(N, p)
b = sum(rand(1,N)<p); % slow but direct

要获得 10000 个实例，请按如下方式运行：

tic
portfolio = zeros(10000, 10);
for ii = 1:10000
  portfolio(ii,:) = randomInt(10, 1000);
end
toc

这在一台普通机器（单线程）上运行了 3.8 秒——当然，获得二项式分布随机变量的方法是减慢它的速度；有功能更有效的统计工具箱，但我没有。如果您增加粒度（例如，通过设置limit=10000），它会减慢更多，因为您增加了生成的随机数样本的数量；上面的limit = 10000循环需要 13.3 秒才能完成。

作为测试，我发现mean(portfolio)'如下std(portfolio)'（带limit=1000）：

100.20  9.446
 99.90  9.547
100.09  9.456
100.00  9.548
100.01  9.356
100.00  9.484
 99.69  9.639
100.06  9.493
 99.94  9.599
100.11  9.453

对我来说，这看起来像是一个非常有说服力的“扁平”分布。我们希望这些数字呈二项式分布，均值为 100，标准差为sqrt(p*(1-p)*n). 在这种情况下，p=0.1我们期望s = 9.4868. 我实际得到的值再次非常接近。

我意识到这对于较大的 值是低效的limit，并且我没有尝试提高效率。我发现当你开发新东西时，清晰度胜过速度。但是例如，您可以预先计算 的累积二项分布p=1./(1:10)，然后进行随机查找；但如果你只做一次，100,000 个实例，它会在一分钟内运行；除非你打算做很多次，否则我不会打扰。但是，如果有人想改进此代码，我会很高兴收到他们的来信。

于 2013-07-23T15:18:52.177 回答
2

最终我解决了这个问题！

我在约翰霍普金斯大学找到了 2 位学者的论文“从单位单纯形中均匀采样” http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf

在论文中，他们概述了朴素算法是如何不起作用的，其方式非常类似于木片对添加到 100 个问题的随机数的回答。然后，他们继续表明 David Schwartz 建议的方法也可能略有偏差，并提出了一种似乎可行的修改算法。

如果您想要总和为 y 的 x 个数字

  • 从 1 到 x+y-1 范围内均匀采样 x-1 个随机数而不进行替换
  • 对它们进行排序
  • 在开头添加一个零，在末尾添加 x+y
  • 区分它们并从每个值中减去 1
  • 如果你想像我一样缩放它们，然后除以 y

我花了一段时间才意识到为什么当原始方法不起作用时这有效，并且归结为获得零权重的概率（正如弗洛里斯在他的回答中强调的那样）。要在原始版本中为除第一个或最后一个权重之外的所有权重获得零权重，您的随机数必须有 2 个相同的值，但对于第一个和最后一个，则随机数为零或最大数将导致权重为零这更有可能。在修改后的算法中，零和最大数字不在随机选择集合中，并且仅当您选择两个连续数字时才会出现零权重，这对于每个位置都是同样可能的。

我在 Matlab 中编码如下

function weights = unbiased_monkey_weights(num_simulations,num_stocks,min_increment)

scaling=1/min_increment;
sample=NaN(num_simulations,num_stocks-1);

    for i=1:num_simulations
      allcomb=randperm(scaling+num_stocks-1);
      sample(i,:)=allcomb(1:num_stocks-1);
    end

temp = [zeros(num_simulations,1),sort(sample,2),ones(num_simulations,1)*(scaling+num_stocks)];
weights = (diff(temp,[],2)-1)/scaling;

end

显然循环有点笨拙，因为我使用的是 2009 版本，所以 randperm 函数只允许您生成整个集合的排列，但是尽管如此，我可以在我笨重的笔记本电脑上在 5 秒内为 1,000 个数字运行 10,000 次模拟，这是足够快。

平均重量现在是正确的，作为快速测试，我复制了生成 3 个数字的木片，这些数字总和为 1，最小增量为 0.01%，而且看起来也正确

在此处输入图像描述

谢谢大家的帮助，我希望这个解决方案在未来对其他人有用

于 2013-07-23T20:03:15.087 回答
1

简单的答案是使用在没有最小增量的情况下运行良好的方案，然后转换问题。一如既往，要小心。有些方法不会产生统一的数字集。

因此，假设我想要 11 个总和为 100 的数字，最小增量为 5。我首先会找到 11 个总和为 45 的数字，样本没有下限（除了零）。我可以使用来自文件交换的工具。最简单的是在区间 [0,45] 中简单地采样 10 个数字。对它们进行排序，然后找出不同之处。

X = diff([0,sort(rand(1,10)),1]*45);

向量 X 是总和为 45 的数字样本。但向量 Y 的总和为 100，最小值为 5。

Y = X + 5;

当然，如果您希望找到具有给定约束的多组数字，这很容易矢量化。

于 2013-07-23T15:32:30.243 回答