C - 在相对于均值的区间内生成随机数
C - generate random numbers within an interval with respect to a mean
我需要在一个区间内生成一组随机数，该区间也恰好有一个平均值。例如，最小值 = 1000，最大值 = 10000，平均值为 7000。我知道如何在一个范围内创建数字，但我正在为求平均值而苦苦挣扎。有我可以使用的功能吗？
使用所谓的接受拒绝方法可以最轻松地完成您正在寻找的内容。
将间隔分成更小的间隔。
指定一个概率密度函数 (PDF)，也可以是一个非常简单的函数，例如阶跃函数。对于高斯分布，您的左右步长将低于中间步长，即（请参见下图，该图像具有更一般的分布）。
在整个区间内生成一个随机数。如果此时生成的数字大于您的 PDF 的值，则拒绝生成的数字。
重复这些步骤，直到获得所需的分数
编辑 1
高斯 PDF 的概念证明。
好的，基本思路如图(a)所示。
- Define/Pick 您的概率密度函数 (PDF)。从统计学上讲，PDF 是随机变量的函数，它描述了在 measurement/experiment 中找到值 x 的概率。函数可以是随机变量
x 的 PDF，如果它满足：1) f(x) >= 0 和 2) 它被归一化（意味着它求和或积分，直到值 1）。
- 获取 PDF 的最大值 (
max) 和 "zero points" (z1 < z2)。一些 PDF 的零点可能在无穷大。在这种情况下，确定截止点 (z1, z2)，其中 PDF(z1>x>z2) < eta 您自己选择 eta。基本上意味着，设置一些小值 eta 然后说你的零点是 PDF(x) 的值小于 eta. 的那些值
- 定义随机生成器的间隔
Ch(z1, z2, max)。这是您生成随机变量的时间间隔。
- 生成一个随机变量 x 使得
z1<x<z2。
- 在
(0, max) 范围内生成第二个不相关的随机变量 y。如果 y 的值小于 PDF(x) 拒绝两个随机生成的值 (x,y) 并返回到步骤 4。如果生成的值 y 大于 PDF(x) 接受值 x 作为分布上随机生成的点并 return 它。
这是为 Gaussian PDF 重现类似行为的代码。
#include "Random.h"
#include <fstream>
using namespace std;
double gaus(double a, double b, double c, double x)
{
return a*exp( -((x-b)*(x-b)/(2*c*c) ));
}
double* random_on_a_gaus_distribution(double inter_a, double inter_b)
{
double res [2];
double a = 1.0; //currently parameters for the Gaussian
double b = 2.0; //are defined here to avoid having
double c = 3.0; //a long function declaration line.
double x = kiss::Ran(inter_a, inter_b);
double y = kiss::Ran(0.0, 1.0);
while (y>gaus(a,b,c,x)) //keep creating values until step 5. is satisfied.
{
x = kiss::Ran(inter_a, inter_b); //this is interval (z1, z2)
y = kiss::Ran(0.0, 1.0); //this is the interval (0, max)
}
res[0] = x;
res[1] = y;
return res; //I return (x,y) for plot reasons, only x is the randomly
} //generated value you're looking for.
void main()
{
double* x;
ofstream f;
f.open("test.txt");
for(int i=0; i<100000; i++)
{
//see bellow how I got -5 and 10 to be my interval (z1, z2)
x = random_on_a_gaus_distribution(-5.0, 10.0);
f << x[0]<<","<<x[1]<<endl;
}
f.close();
}
步骤 1
因此，首先我们在名为 gaus 的函数中定义高斯 PDF 的一般外观。简单的。
然后我们定义一个函数random_on_a_gaus_distribution，它使用定义良好的高斯函数。在 experiment\measurement 中，我们将通过拟合函数得到系数 a, b, c。我为这个例子选择了一些随机的（1、2、3），你可以选择满足你的 HW 分配的那些（即：使高斯均值为 7000 的系数）。
步骤 2 和 3
我用 wolfram mathematica 绘制高斯。使用参数 1、2、3 也可以查看 max 和 (z1, z2) 最合适的值。你可以see the graph yourself。函数的最大值是 1.0，通过称为 eyeballin 的古老科学方法，我估计截止点是 -5.0 和 10.0。
为了使 random_on_a_gaus_distribution 更通用，您可以更严格地执行步骤 2) 并定义 eta，然后在连续的点中计算您的函数，直到 PDF 小于 eta。这样做的危险是你的截止点可能相距很远，对于非常单调的函数来说这可能需要很长时间。此外，您必须自己找到最大值。这通常很棘手，但是一个更简单的问题是最小化函数的负数。对于一般情况，这也可能很棘手，但 "undoable" 并非如此。最简单的方法是像我一样作弊，只对几个函数进行硬编码。
步骤 4 和 5
然后你bash走了。只要不断创造新的和新的点，直到你达到满意的命中率。 请注意 returned 号码 x 是 一个随机数。您将无法在两个连续创建的 x 值或第一个创建的 x 和第百万个值之间找到逻辑 link。
然而，在我们分布的 x_max 区间内接受的 x 值的数量大于在 PDF(x) < PDF(x_max) 区间内创建的 x 值的数量].
这只是意味着您的随机数将在所选区间内以这样的方式加权，即随机变量的较大 PDF 值 x 将对应于在该值附近的小区间内接受的更多随机点比 xi 的任何其他值 PDF(xi)<PDF(x).
我 return 编辑了 x 和 y 以便能够绘制下面的图表，但是您要 return 实际上只是 x。我用 matplotlib 画了图。
最好只显示随机创建的分布变量的直方图。这表明 PDF 函数平均值附近的 x 值最有可能被接受，因此将创建更多具有这些近似值的随机创建变量。
此外，我假设您会对 kiss 随机数生成器的实现感兴趣。 拥有一台非常好的发电机非常重要。我敢说在某种程度上 kiss 可能不会削减它（经常使用 mersene twister）。
Random.h
#pragma once
#include <stdlib.h>
const unsigned RNG_MAX=4294967295;
namespace kiss{
// unsigned int kiss_z, kiss_w, kiss_jsr, kiss_jcong;
unsigned int RanUns();
void RunGen();
double Ran0(int upper_border);
double Ran(double bottom_border, double upper_border);
}
namespace Crand{
double Ran0(int upper_border);
double Ran(double bottom_border, double upper_border);
}
Kiss.cpp
#include "Random.h"
unsigned int kiss_z = 123456789; //od 1 do milijardu
unsigned int kiss_w = 378295763; //od 1 do milijardu
unsigned int kiss_jsr = 294827495; //od 1 do RNG_MAX
unsigned int kiss_jcong = 495749385; //od 0 do RNG_MAX
//KISS99*
//Autor: George Marsaglia
unsigned int kiss::RanUns()
{
kiss_z=36969*(kiss_z&65535)+(kiss_z>>16);
kiss_w=18000*(kiss_w&65535)+(kiss_w>>16);
kiss_jsr^=(kiss_jsr<<13);
kiss_jsr^=(kiss_jsr>>17);
kiss_jsr^=(kiss_jsr<<5);
kiss_jcong=69069*kiss_jcong+1234567;
return (((kiss_z<<16)+kiss_w)^kiss_jcong)+kiss_jsr;
}
void kiss::RunGen()
{
for (int i=0; i<2000; i++)
kiss::RanUns();
}
double kiss::Ran0(int upper_border)
{
unsigned velicinaIntervala = RNG_MAX / upper_border;
unsigned granicaIzbora= velicinaIntervala*upper_border;
unsigned slucajniBroj = kiss::RanUns();
while(slucajniBroj>=granicaIzbora)
slucajniBroj = kiss::RanUns();
return slucajniBroj/velicinaIntervala;
}
double kiss::Ran (double bottom_border, double upper_border)
{
return bottom_border+(upper_border-bottom_border)*kiss::Ran0(100000)/(100001.0);
}
此外还有标准的 C 随机生成器：
CRands.cpp
#include "Random.h"
//standardni pseudo random generatori iz C-a
double Crand::Ran0(int upper_border)
{
return rand()%upper_border;
}
double Crand::Ran (double bottom_border, double upper_border)
{
return (upper_border-bottom_border)*rand()/((double)RAND_MAX+1);
}
也值得对上面的 (b) 图进行评论。当你有一个非常糟糕的 PDF 时，PDF(x) 会在大数字和非常小的数字之间有很大差异。
问题在于区间区域 Ch(x) 将很好地匹配 PDF 的极值，但是由于我们为 PDF(x) 的小值创建了一个随机变量 y 作为出色地;接受该值的机会很小！更有可能的是，生成的 y 值在该点将始终大于 PDF(x)。这意味着您将花费大量周期来创建不会被选中的数字，并且您选择的所有随机数都将在本地绑定到 PDF 的 max。
这就是为什么不在所有地方使用相同的 Ch(x) 间隔，而是定义一组参数化的间隔通常很有用。然而，这给代码增加了相当多的复杂性。
你在哪里设置你的限制？如何处理边缘情况？何时以及如何确定您确实需要突然使用这种方法？计算 max 现在可能并不那么简单，具体取决于您最初设想的执行此操作的方法。
此外，现在您必须纠正这样一个事实，即在 Ch(x) 框高度较低的区域更容易接受更多数字，这会扭曲原始 PDF。
这可以通过权衡在较低边界中创建的数字与较高和较低边界的高度比来纠正，基本上你再重复一次y步骤。从 0 到 1 创建一个随机数 z 并将其与比率 lower_height/higher_height 进行比较，保证 <1。如果 z 小于比率：接受 x，如果大于则拒绝。
也可以通过编写一个接受对象指针的函数来概括所呈现的代码。通过定义您自己的 class 即 function，它通常会描述函数，在某个点有一个 eval 方法，能够存储您的参数，计算并存储它自己的 max/min 值和 zero/cutoff 点，你不必像我一样通过或在函数中定义它们。
祝你玩得开心！
tl;dr：提高均匀的 0 到 1 分布的 (1 - m) / m 次方，其中 m 是所需的均值（在 0 和 1 之间） . Shift/scale 根据需要。
我很好奇如何实现这一点。我认为梯形是最简单的方法，但是你会受到限制，因为你可以获得的最极端的平均值是三角形，这并不是那么极端。数学开始变得困难，所以我恢复了一种似乎很有效的纯经验方法。
无论如何，对于分布，如何从均匀的 [0, 1) 分布开始并将值提高到某个任意幂。将它们平方，分布向右移动。对它们进行平方根，然后它们向左移动。你可以走到任何你想要的极端，并尽可能地推动分布。
def randompow(p):
return random.random() ** p
（所有内容都写在Python中，但应该很容易翻译。如果有不清楚的地方，尽管问。random.random() returns 从0到1浮动）
那么，我们如何调整这个功率呢？那么，均值似乎如何随着不同的幂而变化？
看起来像某种 S 形曲线。有很多 sigmoid functions，但双曲正切似乎工作得很好。
那里不是 100%，让我们尝试在 X 方向缩放它...
# x are the values from -3 to 3 (log transformed from the powers used)
# y are the empirically-determined means given all those powers
def fitter(tanscale):
xsc = tanscale * x
sigtan = np.tanh(xsc)
sigtan = (1 - sigtan) / 2
resid = sigtan - y
return sum(resid**2)
fit = scipy.optimize.minimize(fitter, 1)
钳工说最好的比例因子是1.1514088816214016。残差实际上很低，所以听起来不错。
实现我没有谈到的所有数学的逆看起来像：
def distpow(mean):
p = 1 - (mean * 2)
p = np.arctanh(p) / 1.1514088816214016
return 10**p
这使我们能够在第一个函数中使用以获得分布的任何均值。工厂函数可以 return 一种从具有所需均值的分布中生成一堆数字的方法
def randommean(mean):
p = distpow(mean)
def f():
return random.random() ** p
return f
怎么样？相当不错的小数点后 3-4 位：
for x in [0.01, 0.1, 0.2, 0.4, 0.5, 0.6, 0.8, 0.9, 0.99]:
f = randommean(x)
# sample the distribution 10 million times
mean = np.mean([f() for _ in range(10000000)])
print('Target mean: {:0.6f}, actual: {:0.6f}'.format(x, mean))
Target mean: 0.010000, actual: 0.010030
Target mean: 0.100000, actual: 0.100122
Target mean: 0.200000, actual: 0.199990
Target mean: 0.400000, actual: 0.400051
Target mean: 0.500000, actual: 0.499905
Target mean: 0.600000, actual: 0.599997
Target mean: 0.800000, actual: 0.799999
Target mean: 0.900000, actual: 0.899972
Target mean: 0.990000, actual: 0.989996
一个更简洁的函数，它只给你一个给定平均值的值（不是工厂函数）：
def randommean(m):
p = np.arctanh(1 - (2 * m)) / 1.1514088816214016
return random.random() ** (10 ** p)
编辑： 拟合平均值的自然对数而不是 log10 得到的残差接近 0.5。做一些数学运算来简化 arctanh 给出：
def randommean(m):
'''Return a value from the distribution 0 to 1 with average *m*'''
return random.random() ** ((1 - m) / m)
从这里开始，移动、重新缩放和四舍五入分布应该相当容易。截断为整数可能最终会将平均值移动 1（或半个单位？），因此这是一个未解决的问题（如果重要的话）。
您只需定义在 [1000, 7000] 中运行的 2 个分布 dist1 和在 [7000, 10000] 中运行的 dist2。
我们称 m1 为 dist1 的平均值，m2 为 dist2 的平均值。
您正在寻找 dist1 和 dist2 之间的混合，其平均值为 7000。
您必须调整权重 (w1, w2 = 1-w1) 例如：
7000 = w1 * m1 + w2 * m2
这导致：
w1 = (m2 - 7000) / (m2 - m1)
使用 OpenTURNS 库，代码如下所示：
import openturns as ot
dist1 = ot.Uniform(1000, 7000)
dist2 = ot.Uniform(7000, 10000)
m1 = dist1.getMean()[0]
m2 = dist2.getMean()[0]
w = (m2 - 7000) / (m2 - m1)
dist = ot.Mixture([dist1, dist2], [w, 1 - w])
print ("Mean of dist = ", dist.getMean())
>>> Mean of dist = [7000]
现在调用dist.getSample(N)即可抽取大小为N的样本。例如：
print(dist.getSample(10))
>>> [ X0 ]
0 : [ 3019.97 ]
1 : [ 7682.17 ]
2 : [ 9035.1 ]
3 : [ 8873.59 ]
4 : [ 5217.08 ]
5 : [ 6329.67 ]
6 : [ 9791.22 ]
7 : [ 7786.76 ]
8 : [ 7046.59 ]
9 : [ 7088.48 ]
我需要在一个区间内生成一组随机数，该区间也恰好有一个平均值。例如，最小值 = 1000，最大值 = 10000，平均值为 7000。我知道如何在一个范围内创建数字，但我正在为求平均值而苦苦挣扎。有我可以使用的功能吗？
使用所谓的接受拒绝方法可以最轻松地完成您正在寻找的内容。
将间隔分成更小的间隔。 指定一个概率密度函数 (PDF)，也可以是一个非常简单的函数，例如阶跃函数。对于高斯分布，您的左右步长将低于中间步长，即（请参见下图，该图像具有更一般的分布）。
在整个区间内生成一个随机数。如果此时生成的数字大于您的 PDF 的值，则拒绝生成的数字。
重复这些步骤，直到获得所需的分数
编辑 1
高斯 PDF 的概念证明。
好的，基本思路如图(a)所示。
- Define/Pick 您的概率密度函数 (PDF)。从统计学上讲，PDF 是随机变量的函数，它描述了在 measurement/experiment 中找到值 x 的概率。函数可以是随机变量
x的 PDF，如果它满足：1)f(x) >= 0和 2) 它被归一化（意味着它求和或积分，直到值 1）。 - 获取 PDF 的最大值 (
max) 和 "zero points" (z1 < z2)。一些 PDF 的零点可能在无穷大。在这种情况下，确定截止点(z1, z2)，其中PDF(z1>x>z2) < eta您自己选择eta。基本上意味着，设置一些小值eta然后说你的零点是PDF(x)的值小于 eta. 的那些值
- 定义随机生成器的间隔
Ch(z1, z2, max)。这是您生成随机变量的时间间隔。 - 生成一个随机变量 x 使得
z1<x<z2。 - 在
(0, max)范围内生成第二个不相关的随机变量y。如果y的值小于PDF(x)拒绝两个随机生成的值(x,y)并返回到步骤 4。如果生成的值y大于PDF(x)接受值x作为分布上随机生成的点并return它。
这是为 Gaussian PDF 重现类似行为的代码。
#include "Random.h"
#include <fstream>
using namespace std;
double gaus(double a, double b, double c, double x)
{
return a*exp( -((x-b)*(x-b)/(2*c*c) ));
}
double* random_on_a_gaus_distribution(double inter_a, double inter_b)
{
double res [2];
double a = 1.0; //currently parameters for the Gaussian
double b = 2.0; //are defined here to avoid having
double c = 3.0; //a long function declaration line.
double x = kiss::Ran(inter_a, inter_b);
double y = kiss::Ran(0.0, 1.0);
while (y>gaus(a,b,c,x)) //keep creating values until step 5. is satisfied.
{
x = kiss::Ran(inter_a, inter_b); //this is interval (z1, z2)
y = kiss::Ran(0.0, 1.0); //this is the interval (0, max)
}
res[0] = x;
res[1] = y;
return res; //I return (x,y) for plot reasons, only x is the randomly
} //generated value you're looking for.
void main()
{
double* x;
ofstream f;
f.open("test.txt");
for(int i=0; i<100000; i++)
{
//see bellow how I got -5 and 10 to be my interval (z1, z2)
x = random_on_a_gaus_distribution(-5.0, 10.0);
f << x[0]<<","<<x[1]<<endl;
}
f.close();
}
步骤 1
因此，首先我们在名为 gaus 的函数中定义高斯 PDF 的一般外观。简单的。
然后我们定义一个函数random_on_a_gaus_distribution，它使用定义良好的高斯函数。在 experiment\measurement 中，我们将通过拟合函数得到系数 a, b, c。我为这个例子选择了一些随机的（1、2、3），你可以选择满足你的 HW 分配的那些（即：使高斯均值为 7000 的系数）。
步骤 2 和 3
我用 wolfram mathematica 绘制高斯。使用参数 1、2、3 也可以查看 max 和 (z1, z2) 最合适的值。你可以see the graph yourself。函数的最大值是 1.0，通过称为 eyeballin 的古老科学方法，我估计截止点是 -5.0 和 10.0。
为了使 random_on_a_gaus_distribution 更通用，您可以更严格地执行步骤 2) 并定义 eta，然后在连续的点中计算您的函数，直到 PDF 小于 eta。这样做的危险是你的截止点可能相距很远，对于非常单调的函数来说这可能需要很长时间。此外，您必须自己找到最大值。这通常很棘手，但是一个更简单的问题是最小化函数的负数。对于一般情况，这也可能很棘手，但 "undoable" 并非如此。最简单的方法是像我一样作弊，只对几个函数进行硬编码。
步骤 4 和 5
然后你bash走了。只要不断创造新的和新的点，直到你达到满意的命中率。 请注意 returned 号码 x 是 一个随机数。您将无法在两个连续创建的 x 值或第一个创建的 x 和第百万个值之间找到逻辑 link。
然而，在我们分布的 x_max 区间内接受的 x 值的数量大于在 PDF(x) < PDF(x_max) 区间内创建的 x 值的数量].
这只是意味着您的随机数将在所选区间内以这样的方式加权，即随机变量的较大 PDF 值 x 将对应于在该值附近的小区间内接受的更多随机点比 xi 的任何其他值 PDF(xi)<PDF(x).
我 return 编辑了 x 和 y 以便能够绘制下面的图表，但是您要 return 实际上只是 x。我用 matplotlib 画了图。
最好只显示随机创建的分布变量的直方图。这表明 PDF 函数平均值附近的 x 值最有可能被接受，因此将创建更多具有这些近似值的随机创建变量。
此外，我假设您会对 kiss 随机数生成器的实现感兴趣。 拥有一台非常好的发电机非常重要。我敢说在某种程度上 kiss 可能不会削减它（经常使用 mersene twister）。
Random.h
#pragma once
#include <stdlib.h>
const unsigned RNG_MAX=4294967295;
namespace kiss{
// unsigned int kiss_z, kiss_w, kiss_jsr, kiss_jcong;
unsigned int RanUns();
void RunGen();
double Ran0(int upper_border);
double Ran(double bottom_border, double upper_border);
}
namespace Crand{
double Ran0(int upper_border);
double Ran(double bottom_border, double upper_border);
}
Kiss.cpp
#include "Random.h"
unsigned int kiss_z = 123456789; //od 1 do milijardu
unsigned int kiss_w = 378295763; //od 1 do milijardu
unsigned int kiss_jsr = 294827495; //od 1 do RNG_MAX
unsigned int kiss_jcong = 495749385; //od 0 do RNG_MAX
//KISS99*
//Autor: George Marsaglia
unsigned int kiss::RanUns()
{
kiss_z=36969*(kiss_z&65535)+(kiss_z>>16);
kiss_w=18000*(kiss_w&65535)+(kiss_w>>16);
kiss_jsr^=(kiss_jsr<<13);
kiss_jsr^=(kiss_jsr>>17);
kiss_jsr^=(kiss_jsr<<5);
kiss_jcong=69069*kiss_jcong+1234567;
return (((kiss_z<<16)+kiss_w)^kiss_jcong)+kiss_jsr;
}
void kiss::RunGen()
{
for (int i=0; i<2000; i++)
kiss::RanUns();
}
double kiss::Ran0(int upper_border)
{
unsigned velicinaIntervala = RNG_MAX / upper_border;
unsigned granicaIzbora= velicinaIntervala*upper_border;
unsigned slucajniBroj = kiss::RanUns();
while(slucajniBroj>=granicaIzbora)
slucajniBroj = kiss::RanUns();
return slucajniBroj/velicinaIntervala;
}
double kiss::Ran (double bottom_border, double upper_border)
{
return bottom_border+(upper_border-bottom_border)*kiss::Ran0(100000)/(100001.0);
}
此外还有标准的 C 随机生成器： CRands.cpp
#include "Random.h"
//standardni pseudo random generatori iz C-a
double Crand::Ran0(int upper_border)
{
return rand()%upper_border;
}
double Crand::Ran (double bottom_border, double upper_border)
{
return (upper_border-bottom_border)*rand()/((double)RAND_MAX+1);
}
也值得对上面的 (b) 图进行评论。当你有一个非常糟糕的 PDF 时，PDF(x) 会在大数字和非常小的数字之间有很大差异。
问题在于区间区域 Ch(x) 将很好地匹配 PDF 的极值，但是由于我们为 PDF(x) 的小值创建了一个随机变量 y 作为出色地;接受该值的机会很小！更有可能的是，生成的 y 值在该点将始终大于 PDF(x)。这意味着您将花费大量周期来创建不会被选中的数字，并且您选择的所有随机数都将在本地绑定到 PDF 的 max。
这就是为什么不在所有地方使用相同的 Ch(x) 间隔，而是定义一组参数化的间隔通常很有用。然而，这给代码增加了相当多的复杂性。
你在哪里设置你的限制？如何处理边缘情况？何时以及如何确定您确实需要突然使用这种方法？计算 max 现在可能并不那么简单，具体取决于您最初设想的执行此操作的方法。
此外，现在您必须纠正这样一个事实，即在 Ch(x) 框高度较低的区域更容易接受更多数字，这会扭曲原始 PDF。
这可以通过权衡在较低边界中创建的数字与较高和较低边界的高度比来纠正，基本上你再重复一次y步骤。从 0 到 1 创建一个随机数 z 并将其与比率 lower_height/higher_height 进行比较，保证 <1。如果 z 小于比率：接受 x，如果大于则拒绝。
也可以通过编写一个接受对象指针的函数来概括所呈现的代码。通过定义您自己的 class 即 function，它通常会描述函数，在某个点有一个 eval 方法，能够存储您的参数，计算并存储它自己的 max/min 值和 zero/cutoff 点，你不必像我一样通过或在函数中定义它们。
祝你玩得开心！
tl;dr：提高均匀的 0 到 1 分布的 (1 - m) / m 次方，其中 m 是所需的均值（在 0 和 1 之间） . Shift/scale 根据需要。
我很好奇如何实现这一点。我认为梯形是最简单的方法，但是你会受到限制，因为你可以获得的最极端的平均值是三角形，这并不是那么极端。数学开始变得困难，所以我恢复了一种似乎很有效的纯经验方法。
无论如何，对于分布，如何从均匀的 [0, 1) 分布开始并将值提高到某个任意幂。将它们平方，分布向右移动。对它们进行平方根，然后它们向左移动。你可以走到任何你想要的极端，并尽可能地推动分布。
def randompow(p):
return random.random() ** p
（所有内容都写在Python中，但应该很容易翻译。如果有不清楚的地方，尽管问。random.random() returns 从0到1浮动）
那么，我们如何调整这个功率呢？那么，均值似乎如何随着不同的幂而变化？
看起来像某种 S 形曲线。有很多 sigmoid functions，但双曲正切似乎工作得很好。
那里不是 100%，让我们尝试在 X 方向缩放它...
# x are the values from -3 to 3 (log transformed from the powers used)
# y are the empirically-determined means given all those powers
def fitter(tanscale):
xsc = tanscale * x
sigtan = np.tanh(xsc)
sigtan = (1 - sigtan) / 2
resid = sigtan - y
return sum(resid**2)
fit = scipy.optimize.minimize(fitter, 1)
钳工说最好的比例因子是1.1514088816214016。残差实际上很低，所以听起来不错。
实现我没有谈到的所有数学的逆看起来像：
def distpow(mean):
p = 1 - (mean * 2)
p = np.arctanh(p) / 1.1514088816214016
return 10**p
这使我们能够在第一个函数中使用以获得分布的任何均值。工厂函数可以 return 一种从具有所需均值的分布中生成一堆数字的方法
def randommean(mean):
p = distpow(mean)
def f():
return random.random() ** p
return f
怎么样？相当不错的小数点后 3-4 位：
for x in [0.01, 0.1, 0.2, 0.4, 0.5, 0.6, 0.8, 0.9, 0.99]:
f = randommean(x)
# sample the distribution 10 million times
mean = np.mean([f() for _ in range(10000000)])
print('Target mean: {:0.6f}, actual: {:0.6f}'.format(x, mean))
Target mean: 0.010000, actual: 0.010030
Target mean: 0.100000, actual: 0.100122
Target mean: 0.200000, actual: 0.199990
Target mean: 0.400000, actual: 0.400051
Target mean: 0.500000, actual: 0.499905
Target mean: 0.600000, actual: 0.599997
Target mean: 0.800000, actual: 0.799999
Target mean: 0.900000, actual: 0.899972
Target mean: 0.990000, actual: 0.989996
一个更简洁的函数，它只给你一个给定平均值的值（不是工厂函数）：
def randommean(m):
p = np.arctanh(1 - (2 * m)) / 1.1514088816214016
return random.random() ** (10 ** p)
编辑： 拟合平均值的自然对数而不是 log10 得到的残差接近 0.5。做一些数学运算来简化 arctanh 给出：
def randommean(m):
'''Return a value from the distribution 0 to 1 with average *m*'''
return random.random() ** ((1 - m) / m)
从这里开始，移动、重新缩放和四舍五入分布应该相当容易。截断为整数可能最终会将平均值移动 1（或半个单位？），因此这是一个未解决的问题（如果重要的话）。
您只需定义在 [1000, 7000] 中运行的 2 个分布 dist1 和在 [7000, 10000] 中运行的 dist2。
我们称 m1 为 dist1 的平均值，m2 为 dist2 的平均值。
您正在寻找 dist1 和 dist2 之间的混合，其平均值为 7000。
您必须调整权重 (w1, w2 = 1-w1) 例如：
7000 = w1 * m1 + w2 * m2
这导致：
w1 = (m2 - 7000) / (m2 - m1)
使用 OpenTURNS 库，代码如下所示：
import openturns as ot
dist1 = ot.Uniform(1000, 7000)
dist2 = ot.Uniform(7000, 10000)
m1 = dist1.getMean()[0]
m2 = dist2.getMean()[0]
w = (m2 - 7000) / (m2 - m1)
dist = ot.Mixture([dist1, dist2], [w, 1 - w])
print ("Mean of dist = ", dist.getMean())
>>> Mean of dist = [7000]
现在调用dist.getSample(N)即可抽取大小为N的样本。例如：
print(dist.getSample(10))
>>> [ X0 ]
0 : [ 3019.97 ]
1 : [ 7682.17 ]
2 : [ 9035.1 ]
3 : [ 8873.59 ]
4 : [ 5217.08 ]
5 : [ 6329.67 ]
6 : [ 9791.22 ]
7 : [ 7786.76 ]
8 : [ 7046.59 ]
9 : [ 7088.48 ]