matlab - 以最小增量 : Matlab 加 1 的随机数

Question

仔细阅读了上一个问题 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%.....

编辑以反射(reflect) Floris 的出色回答并澄清

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

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% 的增量，我当前的代码无法正确处理，因此我是如何偶然发现原始问题的计算器

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

这里是原始研究的一个例子

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

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

Answer 1

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

另一件需要注意的事情是，在你的分布中获得极值的概率是你想要的概率的一半:如果你有从 0 到 1000 的均匀分布的数字，并且你round 他们，舍入到 0 的数字在区间 [0 0.5> 内；舍入到 1 的那些来自 [0.5 1.5> - 两倍大。最后一个数字(四舍五入到 1000)也是来自一个较小的区间:[999.5 1000]。因此，您不会像您想的那样经常获得第一个和最后一个数字。如果你使用 floor 而不是 round 我想你会得到你期望的答案。

编辑

我进一步考虑了这个问题，并想出了一个缓慢但(我认为)准确的方法来做到这一点。基本思想是这样的:

1. 用整数来思考；不要以 0.001 的步长划分区间 0 - 1，而是以整数步长划分区间 0 - 1000
2. 如果我们尝试将 N 分成 m 个间隔，则步长的平均大小应为 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 个实例，它会在一分钟内运行；除非你打算多次这样做，否则我不会打扰。但如果有人想改进这段代码，我很乐意听取他们的意见。