如何生成满足python中特定均值和中位数的随机数？
How to generate random numbers to satisfy a specific mean and median in python?
我想生成 n 个随机数，例如 n=200，其中可能值的范围在 2 到 40 之间，平均值为 12，中位数为 6.5。
我到处搜索，找不到解决办法。我尝试了以下脚本，它适用于小数字，例如 20，对于大数字，它需要很长时间并返回结果。
n=200
x = np.random.randint(0,1,size=n) # initalisation only
while True:
if x.mean() == 12 and np.median(x) == 6.5:
break
else:
x=np.random.randint(2,40,size=n)
谁能帮助我改进它，即使在 n=5000 左右时也能快速获得结果？
获得真正接近您想要的结果的一种方法是生成两个长度为 100 的独立随机范围，满足您的中位数约束并包括所有期望的数字范围。然后通过连接数组，平均值将在 12 左右但不完全等于 12。但是因为它只是意味着你正在处理你可以通过调整这些数组之一来简单地生成你的预期结果。
In [162]: arr1 = np.random.randint(2, 7, 100)
In [163]: arr2 = np.random.randint(7, 40, 100)
In [164]: np.mean(np.concatenate((arr1, arr2)))
Out[164]: 12.22
In [166]: np.median(np.concatenate((arr1, arr2)))
Out[166]: 6.5
以下是一个矢量化且非常优化的解决方案，针对任何其他使用 for 循环或 python 级代码通过约束随机序列创建的解决方案：
import numpy as np
import math
def gen_random():
arr1 = np.random.randint(2, 7, 99)
arr2 = np.random.randint(7, 40, 99)
mid = [6, 7]
i = ((np.sum(arr1 + arr2) + 13) - (12 * 200)) / 40
decm, intg = math.modf(i)
args = np.argsort(arr2)
arr2[args[-41:-1]] -= int(intg)
arr2[args[-1]] -= int(np.round(decm * 40))
return np.concatenate((arr1, mid, arr2))
演示：
arr = gen_random()
print(np.median(arr))
print(arr.mean())
6.5
12.0
函数背后的逻辑：
为了让我们有一个符合该标准的随机数组，我们可以将 3 个数组连接在一起 arr1、mid 和 arr2。 arr1 和 arr2 各有 99 个项目，mid 有 2 个项目 6 和 7，因此最终结果为 6.5 作为中位数。现在我们创建两个随机数组，每个数组的长度都是 99。要使结果的均值为 12，我们需要做的就是找出当前总和与 12 * 200 之间的差值，然后从 N 个最大的数字中减去结果在这种情况下，我们可以从 arr2 中选择它们并使用 N=50.
编辑：
如果结果中的浮点数不是问题，您实际上可以缩短函数如下：
import numpy as np
import math
def gen_random():
arr1 = np.random.randint(2, 7, 99).astype(np.float)
arr2 = np.random.randint(7, 40, 99).astype(np.float)
mid = [6, 7]
i = ((np.sum(arr1 + arr2) + 13) - (12 * 200)) / 40
args = np.argsort(arr2)
arr2[args[-40:]] -= i
return np.concatenate((arr1, mid, arr2))
如果您有一堆具有正确中值和均值的较小数组，您可以将它们组合起来生成一个更大的数组。
所以...您可以像您目前所做的那样预先生成较小的数组，然后将它们随机组合以获得更大的 n。当然，这会导致有偏差的随机样本，但听起来您只是想要近似随机的东西。
这是工作 (py3) 代码，它生成一个大小为 5000 的样本，具有您想要的属性，它是从大小为 4、6、8、10、...、18 的较小样本构建的。
请注意，我更改了较小的随机样本的构建方式：如果中位数为 6.5，一半的数字必须 <= 6，另一半必须 >= 7，因此我们独立生成这两半。这大大加快了速度。
import collections
import numpy as np
import random
rs = collections.defaultdict(list)
for i in range(50):
n = random.randrange(4, 20, 2)
while True:
x=np.append(np.random.randint(2, 7, size=n//2), np.random.randint(7, 41, size=n//2))
if x.mean() == 12 and np.median(x) == 6.5:
break
rs[len(x)].append(x)
def random_range(n):
if n % 2:
raise AssertionError("%d must be even" % n)
r = []
while n:
i = random.randrange(4, min(20, n+1), 2)
# Don't be left with only 2 slots left.
if n - i == 2: continue
xs = random.choice(rs[i])
r.extend(xs)
n -= i
random.shuffle(r)
return r
xs = np.array(random_range(5000))
print([(i, list(xs).count(i)) for i in range(2, 41)])
print(len(xs))
print(xs.mean())
print(np.median(xs))
输出：
[(2, 620), (3, 525), (4, 440), (5, 512), (6, 403), (7, 345), (8, 126), (9, 111), (10, 78), (11, 25), (12, 48), (13, 61), (14, 117), (15, 61), (16, 62), (17, 116), (18, 49), (19, 73), (20, 88), (21, 48), (22, 68), (23, 46), (24, 75), (25, 77), (26, 49), (27, 83), (28, 61), (29, 28), (30, 59), (31, 73), (32, 51), (33, 113), (34, 72), (35, 33), (36, 51), (37, 44), (38, 25), (39, 38), (40, 46)]
5000
12.0
6.5
输出的第一行显示最终数组中有 620 个 2、52 个 3、440 个 4 等。
在这里，您想要一个小于平均值的中值。这意味着均匀分布是不合适的：您需要很多小值而较少大值。
具体来说，您希望小于或等于 6 的值与大于或等于 7 的值的数量一样多。
确保中位数为 6.5 的一种简单方法是在 [2 - 6] 范围内使用与 [7 - 40] 范围内相同数量的值。如果您在两个范围内都选择均匀分布，则理论平均值为 13.75，与所需的 12 相差不远。
权重的微小变化可以使理论平均值更接近：如果我们使用 [5, 4, 3, 2, 1, 1, ..., 1 ] 作为 random.choices 的 [ 7, 8, ..., 40 ] 范围，我们发现该范围的理论平均值为 19.98，足够接近预期的 20.
示例代码：
>>> pop1 = list(range(2, 7))
>>> pop2 = list(range(7, 41))
>>> w2 = [ 5, 4, 3, 2 ] + ( [1] * 30)
>>> r1 = random.choices(pop1, k=2500)
>>> r2 = random.choices(pop2, w2, k=2500)
>>> r = r1 + r2
>>> random.shuffle(r)
>>> statistics.mean(r)
12.0358
>>> statistics.median(r)
6.5
>>>
所以我们现在有一个 5000 个值的分布，其中位数恰好为 6.5，平均值为 12.0358（这个 是 随机的，另一个测试会给出稍微不同的结果价值）。如果我们想要 12 的精确平均值，我们只需要调整一些值。这里 sum(r) 是 60179，而它应该是 60000，所以我们必须减少 175 个既不是 2（会超出范围）也不是 7（会改变中位数）的值。
最后，一个可能的生成器函数可能是：
def gendistrib(n):
if n % 2 != 0 :
raise ValueError("gendistrib needs an even parameter")
n2 = n//2 # n / 2 in Python 2
pop1 = list(range(2, 7)) # lower range
pop2 = list(range(7, 41)) # upper range
w2 = [ 5, 4, 3, 2 ] + ( [1] * 30) # weights for upper range
r1 = random.choices(pop1, k=n2) # lower part of the distrib.
r2 = random.choices(pop2, w2, k=n2) # upper part
r = r1 + r2
random.shuffle(r) # randomize order
# time to force an exact mean
tot = sum(r)
expected = 12 * n
if tot > expected: # too high: decrease some values
for i, val in enumerate(r):
if val != 2 and val != 7:
r[i] = val - 1
tot -= 1
if tot == expected:
random.shuffle(r) # shuffle again the decreased values
break
elif tot < expected: # too low: increase some values
for i, val in enumerate(r):
if val != 6 and val != 40:
r[i] = val + 1
tot += 1
if tot == expected:
random.shuffle(r) # shuffle again the increased values
break
return r
它真的很快：我可以 timeit gendistrib(10000) 不到 0.02 秒。但它不应该用于小分布（小于 1000）
好的，您正在查看具有不少于 4 个参数的分布 - 其中两个定义范围，另外两个负责所需的均值和中位数。
我可以从脑海中想到两种可能性：
截断正态分布，看here for details. You have already range defined, and have to recover μ and σ from mean and median. It will require solving couple of nonlinear equation, but quite doable in python. Sampling could be done using https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.truncnorm.html
4 参数 Beta 分布，参见 here for details. Again, recovering α and β in Beta distribution from mean and median will require solving couple of non-linear equations. Knowing them sampling would be easy via https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.beta.html
更新
在这里你可以如何从均值到 mu 截断正态：Truncated normal with a given mean
虽然这个 post 已经有一个公认的答案，但我想贡献一个通用的非整数方法。它不需要循环或测试。这个想法是采用紧凑支持的PDF。采纳 Kasrâmvd 已接受答案的想法，在左右间隔中进行两次分布。选择形状参数，使平均值下降到给定值。这里有趣的机会是可以创建一个连续的 PDF，即在间隔连接处没有跳跃。
例如，我选择了 beta 发行版。为了在边界处有有限的 non-zero 值，我选择左边的 beta =1 和右边的 alpha = 1。
查看 PDF 的定义和均值的要求，连续性给出两个方程：
4.5 / alpha = 33.5 / beta2 + 6.5 * alpha / ( alpha + 1 ) + 6.5 + 33.5 * 1 / ( 1 + beta ) = 24
这是一个比较容易求解的二次方程。刚刚使用 scipy.stat.beta like
from scipy.stats import beta
import matplotlib.pyplot as plt
import numpy as np
x1 = np.linspace(2, 6.5, 200 )
x2 = np.linspace(6.5, 40, 200 )
# i use s and t not alpha and beta
s = 1./737 *(np.sqrt(294118) - 418 )
t = 1./99 *(np.sqrt(294118) - 418 )
data1 = beta.rvs(s, 1, loc=2, scale=4.5, size=20000)
data2 = beta.rvs(1, t, loc=6.5, scale=33.5, size=20000)
data = np.concatenate( ( data1, data2 ) )
print np.mean( data1 ), 2 + 4.5 * s/(1.+s)
print np.mean( data2 ), 6.5 + 33.5/(1.+t)
print np.mean( data )
print np.median( data )
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.hist(data1, bins=13, density=True )
ax.hist(data2, bins=67, density=True )
ax.plot( x1, beta.pdf( x1, s, 1, loc=2, scale=4.5 ) )
ax.plot( x2, beta.pdf( x2, 1, t, loc=6.5, scale=33.5 ) )
ax.set_yscale( 'log' )
plt.show()
提供
>> 2.661366939244768 2.6495436216856976
>> 21.297348804473618 21.3504563783143
>> 11.979357871859191
>> 6.5006779033245135
所以结果符合要求，看起来像：
我想生成 n 个随机数，例如 n=200，其中可能值的范围在 2 到 40 之间，平均值为 12，中位数为 6.5。
我到处搜索，找不到解决办法。我尝试了以下脚本，它适用于小数字，例如 20，对于大数字，它需要很长时间并返回结果。
n=200
x = np.random.randint(0,1,size=n) # initalisation only
while True:
if x.mean() == 12 and np.median(x) == 6.5:
break
else:
x=np.random.randint(2,40,size=n)
谁能帮助我改进它，即使在 n=5000 左右时也能快速获得结果？
获得真正接近您想要的结果的一种方法是生成两个长度为 100 的独立随机范围，满足您的中位数约束并包括所有期望的数字范围。然后通过连接数组，平均值将在 12 左右但不完全等于 12。但是因为它只是意味着你正在处理你可以通过调整这些数组之一来简单地生成你的预期结果。
In [162]: arr1 = np.random.randint(2, 7, 100)
In [163]: arr2 = np.random.randint(7, 40, 100)
In [164]: np.mean(np.concatenate((arr1, arr2)))
Out[164]: 12.22
In [166]: np.median(np.concatenate((arr1, arr2)))
Out[166]: 6.5
以下是一个矢量化且非常优化的解决方案，针对任何其他使用 for 循环或 python 级代码通过约束随机序列创建的解决方案：
import numpy as np
import math
def gen_random():
arr1 = np.random.randint(2, 7, 99)
arr2 = np.random.randint(7, 40, 99)
mid = [6, 7]
i = ((np.sum(arr1 + arr2) + 13) - (12 * 200)) / 40
decm, intg = math.modf(i)
args = np.argsort(arr2)
arr2[args[-41:-1]] -= int(intg)
arr2[args[-1]] -= int(np.round(decm * 40))
return np.concatenate((arr1, mid, arr2))
演示：
arr = gen_random()
print(np.median(arr))
print(arr.mean())
6.5
12.0
函数背后的逻辑：
为了让我们有一个符合该标准的随机数组，我们可以将 3 个数组连接在一起 arr1、mid 和 arr2。 arr1 和 arr2 各有 99 个项目，mid 有 2 个项目 6 和 7，因此最终结果为 6.5 作为中位数。现在我们创建两个随机数组，每个数组的长度都是 99。要使结果的均值为 12，我们需要做的就是找出当前总和与 12 * 200 之间的差值，然后从 N 个最大的数字中减去结果在这种情况下，我们可以从 arr2 中选择它们并使用 N=50.
编辑：
如果结果中的浮点数不是问题，您实际上可以缩短函数如下：
import numpy as np
import math
def gen_random():
arr1 = np.random.randint(2, 7, 99).astype(np.float)
arr2 = np.random.randint(7, 40, 99).astype(np.float)
mid = [6, 7]
i = ((np.sum(arr1 + arr2) + 13) - (12 * 200)) / 40
args = np.argsort(arr2)
arr2[args[-40:]] -= i
return np.concatenate((arr1, mid, arr2))
如果您有一堆具有正确中值和均值的较小数组，您可以将它们组合起来生成一个更大的数组。
所以...您可以像您目前所做的那样预先生成较小的数组，然后将它们随机组合以获得更大的 n。当然，这会导致有偏差的随机样本，但听起来您只是想要近似随机的东西。
这是工作 (py3) 代码，它生成一个大小为 5000 的样本，具有您想要的属性，它是从大小为 4、6、8、10、...、18 的较小样本构建的。
请注意，我更改了较小的随机样本的构建方式：如果中位数为 6.5，一半的数字必须 <= 6，另一半必须 >= 7，因此我们独立生成这两半。这大大加快了速度。
import collections
import numpy as np
import random
rs = collections.defaultdict(list)
for i in range(50):
n = random.randrange(4, 20, 2)
while True:
x=np.append(np.random.randint(2, 7, size=n//2), np.random.randint(7, 41, size=n//2))
if x.mean() == 12 and np.median(x) == 6.5:
break
rs[len(x)].append(x)
def random_range(n):
if n % 2:
raise AssertionError("%d must be even" % n)
r = []
while n:
i = random.randrange(4, min(20, n+1), 2)
# Don't be left with only 2 slots left.
if n - i == 2: continue
xs = random.choice(rs[i])
r.extend(xs)
n -= i
random.shuffle(r)
return r
xs = np.array(random_range(5000))
print([(i, list(xs).count(i)) for i in range(2, 41)])
print(len(xs))
print(xs.mean())
print(np.median(xs))
输出：
[(2, 620), (3, 525), (4, 440), (5, 512), (6, 403), (7, 345), (8, 126), (9, 111), (10, 78), (11, 25), (12, 48), (13, 61), (14, 117), (15, 61), (16, 62), (17, 116), (18, 49), (19, 73), (20, 88), (21, 48), (22, 68), (23, 46), (24, 75), (25, 77), (26, 49), (27, 83), (28, 61), (29, 28), (30, 59), (31, 73), (32, 51), (33, 113), (34, 72), (35, 33), (36, 51), (37, 44), (38, 25), (39, 38), (40, 46)]
5000
12.0
6.5
输出的第一行显示最终数组中有 620 个 2、52 个 3、440 个 4 等。
在这里，您想要一个小于平均值的中值。这意味着均匀分布是不合适的：您需要很多小值而较少大值。
具体来说，您希望小于或等于 6 的值与大于或等于 7 的值的数量一样多。
确保中位数为 6.5 的一种简单方法是在 [2 - 6] 范围内使用与 [7 - 40] 范围内相同数量的值。如果您在两个范围内都选择均匀分布，则理论平均值为 13.75，与所需的 12 相差不远。
权重的微小变化可以使理论平均值更接近：如果我们使用 [5, 4, 3, 2, 1, 1, ..., 1 ] 作为 random.choices 的 [ 7, 8, ..., 40 ] 范围，我们发现该范围的理论平均值为 19.98，足够接近预期的 20.
示例代码：
>>> pop1 = list(range(2, 7))
>>> pop2 = list(range(7, 41))
>>> w2 = [ 5, 4, 3, 2 ] + ( [1] * 30)
>>> r1 = random.choices(pop1, k=2500)
>>> r2 = random.choices(pop2, w2, k=2500)
>>> r = r1 + r2
>>> random.shuffle(r)
>>> statistics.mean(r)
12.0358
>>> statistics.median(r)
6.5
>>>
所以我们现在有一个 5000 个值的分布，其中位数恰好为 6.5，平均值为 12.0358（这个 是 随机的，另一个测试会给出稍微不同的结果价值）。如果我们想要 12 的精确平均值，我们只需要调整一些值。这里 sum(r) 是 60179，而它应该是 60000，所以我们必须减少 175 个既不是 2（会超出范围）也不是 7（会改变中位数）的值。
最后，一个可能的生成器函数可能是：
def gendistrib(n):
if n % 2 != 0 :
raise ValueError("gendistrib needs an even parameter")
n2 = n//2 # n / 2 in Python 2
pop1 = list(range(2, 7)) # lower range
pop2 = list(range(7, 41)) # upper range
w2 = [ 5, 4, 3, 2 ] + ( [1] * 30) # weights for upper range
r1 = random.choices(pop1, k=n2) # lower part of the distrib.
r2 = random.choices(pop2, w2, k=n2) # upper part
r = r1 + r2
random.shuffle(r) # randomize order
# time to force an exact mean
tot = sum(r)
expected = 12 * n
if tot > expected: # too high: decrease some values
for i, val in enumerate(r):
if val != 2 and val != 7:
r[i] = val - 1
tot -= 1
if tot == expected:
random.shuffle(r) # shuffle again the decreased values
break
elif tot < expected: # too low: increase some values
for i, val in enumerate(r):
if val != 6 and val != 40:
r[i] = val + 1
tot += 1
if tot == expected:
random.shuffle(r) # shuffle again the increased values
break
return r
它真的很快：我可以 timeit gendistrib(10000) 不到 0.02 秒。但它不应该用于小分布（小于 1000）
好的，您正在查看具有不少于 4 个参数的分布 - 其中两个定义范围，另外两个负责所需的均值和中位数。
我可以从脑海中想到两种可能性：
截断正态分布，看here for details. You have already range defined, and have to recover μ and σ from mean and median. It will require solving couple of nonlinear equation, but quite doable in python. Sampling could be done using https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.truncnorm.html
4 参数 Beta 分布，参见 here for details. Again, recovering α and β in Beta distribution from mean and median will require solving couple of non-linear equations. Knowing them sampling would be easy via https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.beta.html
更新
在这里你可以如何从均值到 mu 截断正态：Truncated normal with a given mean
虽然这个 post 已经有一个公认的答案，但我想贡献一个通用的非整数方法。它不需要循环或测试。这个想法是采用紧凑支持的PDF。采纳 Kasrâmvd 已接受答案的想法，在左右间隔中进行两次分布。选择形状参数，使平均值下降到给定值。这里有趣的机会是可以创建一个连续的 PDF，即在间隔连接处没有跳跃。
例如，我选择了 beta 发行版。为了在边界处有有限的 non-zero 值，我选择左边的 beta =1 和右边的 alpha = 1。 查看 PDF 的定义和均值的要求，连续性给出两个方程：
4.5 / alpha = 33.5 / beta2 + 6.5 * alpha / ( alpha + 1 ) + 6.5 + 33.5 * 1 / ( 1 + beta ) = 24
这是一个比较容易求解的二次方程。刚刚使用 scipy.stat.beta like
from scipy.stats import beta
import matplotlib.pyplot as plt
import numpy as np
x1 = np.linspace(2, 6.5, 200 )
x2 = np.linspace(6.5, 40, 200 )
# i use s and t not alpha and beta
s = 1./737 *(np.sqrt(294118) - 418 )
t = 1./99 *(np.sqrt(294118) - 418 )
data1 = beta.rvs(s, 1, loc=2, scale=4.5, size=20000)
data2 = beta.rvs(1, t, loc=6.5, scale=33.5, size=20000)
data = np.concatenate( ( data1, data2 ) )
print np.mean( data1 ), 2 + 4.5 * s/(1.+s)
print np.mean( data2 ), 6.5 + 33.5/(1.+t)
print np.mean( data )
print np.median( data )
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.hist(data1, bins=13, density=True )
ax.hist(data2, bins=67, density=True )
ax.plot( x1, beta.pdf( x1, s, 1, loc=2, scale=4.5 ) )
ax.plot( x2, beta.pdf( x2, 1, t, loc=6.5, scale=33.5 ) )
ax.set_yscale( 'log' )
plt.show()
提供
>> 2.661366939244768 2.6495436216856976
>> 21.297348804473618 21.3504563783143
>> 11.979357871859191
>> 6.5006779033245135
所以结果符合要求，看起来像：