【问题标题】:Generate random numbers in a specific range with a standard deviation?生成具有标准偏差的特定范围内的随机数?
【发布时间】:2017-01-01 20:54:06
【问题描述】:

我已经知道如何在一个范围内生成随机数。我可以通过使用来做到这一点

rand.nextInt((max - min) + 1) + min;

问题是我还想为这些数字设置标准偏差。数字也需要是正数,并且它们不在 0 和 1 之间

编辑我删除了 ThreadLocalRandom 类,因为我无法在该类中设置种子,并且这些随机数应该可以在不同的系统中重现。

【问题讨论】:

  • @statboy 你必须指定你想要什么样的分布。您的 nextInt 示例生成了一个 discrete uniform distribution,其方差 (StdDev) 已由 minmax 参数确定(请参阅维基百科链接)。
  • 根据您的最新编辑,我删除了我的答案,因为它显然没有解决您编辑的要求(事后)数字也必须是正数,并且它们不在 0 和1。您将需要通过具有这些限制的“标准偏差”来进一步说明您的意思。
  • @StefanZobel 缩放和翻译的beta distribution 提供了很大的灵活性。
  • @StefanZobel 如果 OP 想要将 Beta 迁移到 1 到 c 之间的范围内,请将其缩放 c-1 并加 1。正如您所描述的那样,方差会受到影响——加法常数不会'不要改变它,你用(c-1)^2缩放原始Beta的方差,其中插入符号表示幂。

标签: java random


【解决方案1】:

选择有界分布的标准差(或方差)只能受到取决于所选分布和区间的边界(min, max) 的约束。某些分布可能允许您使方差任意小(例如Beta distribution),其他分布(例如Uniform distribution)一旦设置了界限(min, max),就不允许任何灵活性。在任何情况下,您都永远无法使方差任意大——边界确实可以防止这种情况发生(它们总是会输入分布方差的表达式)。

我将通过一个非常简单的示例来说明这一点,该示例无需任何第三方库即可实现。假设您想要区间 (min, max) 上的对称分布,对称意味着分布的均值 E(X) 位于区间的中间:E(X) = (min + max)/2

x = a + (b - a) * rnd.nextDouble() 中使用Random 的nextDouble 将为您提供区间a <= x < b 中的均匀分布随机变量,该变量具有固定方差Var(X) = (b - a)^2 / 12(不是我们想要的)。

OTH,在相同的间隔 (a, b) 上模拟对称的 triangular distribution 会给我们一个随机变量,其均值相同但只有一半的方差:Var(X) = (b - a)^2 / 24(也是固定的,所以也不是我们想要的)。

带有参数(a < b < c < d) 的对称trapezoidal distribution 位于区间(a, d) 上的均匀分布和三角形分布的中间某处。对称条件意味着d - c = b - a,下面我将距离b - a称为x或“位移”(这个名字是我编的,不是专业术语)。

如果让x 从上方接近 0.0,梯形将开始看起来非常类似于均匀分布,并且其方差将趋于最大可能值(d - a)^2 / 12。如果让x 从下方逼近最大可能值(d - a)/2,梯形将看起来非常类似于对称三角形分布,其方差将逼近(d - a)^2 / 24) 的最小可能值(但请注意,我们应该远离为了不破坏方差公式或我们的梯形算法,这些极值很少)。

因此,考虑到您的目标标准差必须位于(0.2041(d - a), 0.2886(d - a)) 给出的开放范围(大致)内.为方便起见,我们假设 a = min = 2.0d = max = 10.0 为我们提供了这个可能的 stddev 范围:(1.6328, 2.3088)。让我们进一步假设我们想要构建一个 stddev 为 2.0 的分布(当然,它必须在允许的范围内)。

解决这个问题需要 3 个步骤:

1) 我们需要给定的方差公式min, max 和位移的允许值x

2) 我们需要以某种方式“反转”这个表达式,以便为我们的目标方差提供 x 的值

3) 一旦我们知道x 的值,我们必须构造一个具有对称梯形分布的随机变量,参数为(min, max, x)

第 1 步

/**
 * Variance of a symmetric trapezoidal distribution with parameters
 * {@code a < b < c < d} and the length of {@code d - c = b - a}
 * (by symmetry) identified by {@code x}.
 * 
 * @param a support lower bound
 * @param d support upper bound
 * @param x length of {@code d - c = b - a}, constrained to lie in the open
 *          interval {@code (0, (d-a)/2)}
 * @return variance of the symmetric trapezoidal distribution defined by
 *         the triple {@code (a, d, x)}
 */
static double varSymTrapezoid(double a, double d, double x) {
    if (a <= 0.0 || d <= 0.0 || a >= d) {
        throw new IllegalArgumentException();
    }
    if (x <= 0.0 || x >= (d - a) / 2) {
        throw new IllegalArgumentException();
    }
    double b = a + x;
    double c = d - x;
    double b3 = pow(b, 3);
    double c3 = pow(c, 3);
    double ex2p1 = pow(b, 4) / 4 - a * b3 / 3 + pow(a, 4) / 12;
    double ex2p2 = (c3 / 3 - b3 / 3) * (d - c);
    double ex2p3 = pow(c, 4) / 4 - d * c3 / 3 + pow(d, 4) / 12;
    double ex2 = (ex2p1 + ex2p2 + ex2p3) / ((d - b) * (d - c));
    return ex2 - pow((a + d) / 2, 2);
}

请注意,此公式仅对对称梯形分布有效。例如,如果您使用 2.5 (varSymTrapezoid(2.0, 10.0, 2.5)) 的位移调用此方法,它会返回大约为 3.0416 的方差,这太低了(我们需要 4.0),这意味着 2.5 的位移太很多(较高的位移产生较低的方差)。

方差表达式是x 中的四阶多项式,我不想通过解析求解。然而,对于允许范围内的目标x,这个表达式是单调递减的,所以我们可以构造一个函数,使我们的目标方差过零,并通过简单的bisection 解决这个问题。这是

第 2 步

/**
 * Find the displacement {@code x} for the given {@code stddev} by simple
 * bisection.
 * @param min support lower bound
 * @param max support upper bound
 * @param stddev the standard deviation we want
 * @return the length {@code x} of {@code d - c = b - a} that yields a
 * standard deviation roughly equal to {@code stddev}
 */
static double bisect(double min, double max, double stddev) {
    final double eps = 1e-4;
    final double var = pow(stddev, 2);
    int iters = 0;
    double a = eps;
    double b = (max - min) / 2 - eps;
    double x = eps;
    double dx = b - a;

    while (abs(dx) > eps && iters < 150 && eval(min, max, x, var) != 0.0) {
        x = ((a + b) / 2);
        if ((eval(min, max, a, var) * eval(min, max, x, var)) < 0.0) {
            b = x;
            dx = b - a;
        } else {
            a = x;
            dx = b - a;
        }
        iters++;
    }
    if (abs(eval(min, max, x, var)) > eps) {
        throw new RuntimeException("failed to find solution");
    }
    return x;
}

/**
 * Function whose root we want to find.
 */
static double eval(double min, double max, double x, double var) {
    return varSymTrapezoid(min, max, x) - var;
}

调用bisect 方法并使用标准偏差的期望值2.0 (bisect(2.0, 10.0, 2.0)) 为我们提供所需的位移:~ 1.1716。既然知道x 的值,我们要做的最后一件事就是构造一个适当分布的随机变量,它是

第 3 步

概率论的一个众所周知的事实是,两个独立的均匀分布随机变量X1 ~ U[a1, b1]X2 ~ U[a2, b2]之和是一个在区间[a1 + a2, b1 + b2]提供的对称梯形分布随机变量a1 + b2 &lt; a2 + b1(案例 1)或 a2 + b1 &lt; a1 + b2(案例 2)。我们必须避免 a2 + b1 = a1 + b2 的情况(情况 3),因为这样总和具有我们不想要的对称三角形分布。

我们将选择案例 1 (a1 + b2 &lt; a2 + b1)。在这种情况下,b2 - a2 的长度将等于“位移”x

所以,我们所要做的就是选择区间边界 a1、a2、b1 和 b2 使得 a1 + a2 = minb1 + b2 = maxb2 - a2 = x 并且满足上述不等式:

/**
 * Return a pseudorandom double for the symmetric trapezoidal distribution
 * defined by the triple {@code (min, max, x)}
 * @param min support lower bound
 * @param max support upper bound
 * @param x length of {@code max - c = b - min}, constrained to lie in the
 *          open interval {@code (0, (max-min)/2)}
 */
public static double symTrapezoidRandom(double min, double max, double x) {
    final double a1 = 0.5 * min;
    final double a2 = a1;

    final double b1 = max - a2 - x;
    final double b2 = a2 + x;

    if ((a1 + b2) >= (a2 + b1)) {
        throw new IllegalArgumentException();
    }

    double u = a1 + (b1 - a1) * rnd.nextDouble();
    double v = a2 + (b2 - a2) * rnd.nextDouble();
    return u + v;
}

反复调用symTrapezoidRandom(2.0, 10.0, 1.1716) 会为您提供具有所需分布的随机变量。

您可以使用其他更复杂的发行版(例如 Beta)做非常相似的事情。这将为您提供其他(通常更灵活)的可接受差异界限,但您需要像 commons.math 这样的第 3 方库。

代码中的abspowsqrt指的是静态导入的java.lang.Math方法,rnd是java.util.Random的一个实例。

【讨论】:

  • 你能告诉我你从哪里得到对称梯形的方差公式吗?我在任何地方都找不到。
  • @apophis 我用笔和纸手动计算了它。有趣的是,here 给出的公式是错误的,尽管它看起来很合理。也许他们只是交换了变量。 This 更糟糕的是,甚至 E(X) 都不正确。我找不到一个使用(a, b, c, d) 表示的正确方程,所以我不得不自己推导它。
猜你喜欢
  • 1970-01-01
  • 2020-06-17
  • 2010-10-01
  • 2016-05-04
  • 1970-01-01
  • 2015-04-08
  • 1970-01-01
相关资源
最近更新 更多