c# - 获取两个值之间的 n 个不同的随机数，其总和等于给定数

Question

我想在总和为给定数字的范围内找到不同的随机数。

注意:我在 stackoverflow 中发现了类似的问题，但是它们并没有完全解决这个问题(即它们不考虑范围的负 lowerLimit)。

如果我希望随机数的总和等于 1，我只需生成所需的随机数，计算总和并将每个随机数除以总和即可；但是在这里我需要一些不同的东西；我将需要我的随机数加起来不同于 1，但我的随机数仍然必须在给定范围内。

示例:我需要 30 个介于 -50 和 50 之间的不同随机数(非整数)，其中 30 个生成的数字之和必须等于 300；我写了下面的代码，但是当 n 远大于范围 (upperLimit - lowerLimit) 时它不会工作，该函数可能返回范围 [lowerLimit - upperLimit] 之外的数字。对改进当前解决方案有何帮助？

static void Main(string[] args)
{
    var listWeights = GetRandomNumbersWithConstraints(30, 50, -50, 300);
}

private static List<double> GetRandomNumbersWithConstraints(int n, int upperLimit, int lowerLimit, int sum)
{
    if (upperLimit <= lowerLimit || n < 1)
        throw new ArgumentOutOfRangeException();

    Random rand = new Random(Guid.NewGuid().GetHashCode());
    List<double> weight = new List<double>();

    for (int k = 0; k < n; k++)
    {
        //multiply by rand.NextDouble() to avoid duplicates
        double temp = (double)rand.Next(lowerLimit, upperLimit) * rand.NextDouble();

        if (weight.Contains(temp))
            k--;
        else
            weight.Add(temp);
    }

    //divide each element by the sum
    weight = weight.ConvertAll<double>(x => x / weight.Sum());  //here the sum of my weight will be 1 

    return weight.ConvertAll<double>(x => x * sum);
}

编辑 - 澄清

运行当前代码将生成以下 30 个数字，它们相加为 300。但是这些数字不在 -50 和 50 之间

-4.425315699
67.70219958
82.08592061
46.54014109
71.20352208
-9.554070146
37.65032717
-75.77280868
24.68786878
30.89874589
142.0796933
-1.964407284
9.831226893
-15.21652248
6.479463312
49.61283063
118.1853036
-28.35462683
49.82661159
-65.82706541
-29.6865969
-54.5134262
-56.04708803
-84.63783048
-3.18402453
-13.97935982
-44.54265204
112.774348
-2.911427266
-58.94098071

Answer 1

好的，这里是如何完成的

我们将使用 Dirichlet Distribution ，这是 [0...1] 范围内的随机数 x_i 的分布，使得

求和_i x_i = 1

因此，在线性重新缩放后，总和的条件将自动满足。 Dirichlet 分布由 α_i 参数化，但我们假设所有 RN 都来自相同的边缘分布，因此每个索引只有一个参数 α。

对于合理大的 α 值，采样的随机数的均值 =1/n，方差 ~1/(n * α)，因此较大的 α 导致随机值更接近均值。

好的，现在回到重新缩放，

v_i = A + B*x_i

我们必须得到A 和B。正如@HansKesting 正确指出的那样，只有两个自由参数，我们只能满足两个约束，但你有三个。所以我们会严格满足下限约束、和值约束，但偶尔会违反上限约束。在这种情况下，我们只需丢弃整个样本，然后再做一个。

同样，我们有一个旋钮可以转动，α 变大意味着我们接近平均值并且不太可能达到上限。当 α = 1 时，我很少得到好的样本，但当 α = 10 时，我得到接近 40% 的好样本。当 α = 16 时，我得到接近 80% 的好样本。

Dirichlet 采样是通过 Gamma 分布完成的，使用来自 MathDotNet 的代码.

代码，使用 .NET Core 2.1 测试

using System;

using MathNet.Numerics.Distributions;
using MathNet.Numerics.Random;

class Program
{
    static void SampleDirichlet(double alpha, double[] rn)
    {
        if (rn == null)
            throw new ArgumentException("SampleDirichlet:: Results placeholder is null");

        if (alpha <= 0.0)
            throw new ArgumentException($"SampleDirichlet:: alpha {alpha} is non-positive");

        int n = rn.Length;
        if (n == 0)
            throw new ArgumentException("SampleDirichlet:: Results placeholder is of zero size");

        var gamma = new Gamma(alpha, 1.0);

        double sum = 0.0;
        for(int k = 0; k != n; ++k) {
            double v = gamma.Sample();
            sum  += v;
            rn[k] = v;
        }

        if (sum <= 0.0)
            throw new ApplicationException($"SampleDirichlet:: sum {sum} is non-positive");

        // normalize
        sum = 1.0 / sum;
        for(int k = 0; k != n; ++k) {
            rn[k] *= sum;
        }
    }

    static bool SampleBoundedDirichlet(double alpha, double sum, double lo, double hi, double[] rn)
    {
        if (rn == null)
            throw new ArgumentException("SampleDirichlet:: Results placeholder is null");

        if (alpha <= 0.0)
            throw new ArgumentException($"SampleDirichlet:: alpha {alpha} is non-positive");

        if (lo >= hi)
            throw new ArgumentException($"SampleDirichlet:: low {lo} is larger than high {hi}");

        int n = rn.Length;
        if (n == 0)
            throw new ArgumentException("SampleDirichlet:: Results placeholder is of zero size");

        double mean = sum / (double)n;
        if (mean < lo || mean > hi)
            throw new ArgumentException($"SampleDirichlet:: mean value {mean} is not within [{lo}...{hi}] range");

        SampleDirichlet(alpha, rn);

        bool rc = true;
        for(int k = 0; k != n; ++k) {
            double v = lo + (mean - lo)*(double)n * rn[k];
            if (v > hi)
                rc = false;
            rn[k] = v;
        }
        return rc;
    }

    static void Main(string[] args)
    {
        double[] rn = new double [30];

        double lo = -50.0;
        double hi =  50.0;

        double alpha = 10.0;

        double sum = 300.0;

        for(int k = 0; k != 1_000; ++k) {
            var q = SampleBoundedDirichlet(alpha, sum, lo, hi, rn);
            Console.WriteLine($"Rng(BD), v = {q}");
            double s = 0.0;
            foreach(var r in rn) {
                Console.WriteLine($"Rng(BD),     r = {r}");
                s += r;
            }
            Console.WriteLine($"Rng(BD),    summa = {s}");
        }
    }
}

更新

通常，当人们问这样的问题时，有一个隐含的假设/要求——所有的随机数都应该以相同的方式分布。这意味着如果我从采样数组中为索引为 0 的项目绘制边际概率密度函数 (PDF)，我将获得与为数组中的最后一个项目绘制边际概率密度函数相同的分布。人们通常对随机数组进行采样，将其传递给其他例程来做一些有趣的事情。如果项目 0 的边际 PDF 与最后一个索引项目的边际 PDF 不同，那么仅恢复数组将产生与使用此类随机值的代码截然不同的结果。

在这里，我使用我的抽样例程绘制了原始条件 ([-50...50] sum=300) 下第 0 项和最后一项 (#29) 的随机数分布。看起来很相似，不是吗？

好的，这是您的采样例程的图片，相同的原始条件([-50...50] sum=300)，相同的样本数

更新二

用户应该检查采样例程的返回值，并在(且仅当)返回值为真时接受和使用采样数组。这是接受/拒绝方法。作为说明，下面是用于直方图样本的代码:

        int[] hh = new int[100]; // histogram allocated

        var s = 1.0; // step size
        int k = 0;   // good samples counter
        for( ;; ) {
            var q = SampleBoundedDirichlet(alpha, sum, lo, hi, rn);
            if (q) // good sample, accept it
            {
                var v = rn[0]; // any index, 0 or 29 or ....
                var i = (int)((v - lo) / s);
                i = System.Math.Max(i, 0);
                i = System.Math.Min(i, hh.Length-1);
                hh[i] += 1;

                ++k;
                if (k == 100000) // required number of good samples reached
                    break;
            }
        }
        for(k = 0; k != hh.Length; ++k)
        {
            var x = lo + (double)k * s + 0.5*s;
            var v = hh[k];
            Console.WriteLine($"{x}     {v}");
        }