三角分布应该可以执行您想要的任务，因为它将三个参数（最小值、众数、最大值）作为与您的标准相匹配的输入。您可以考虑其他分布，如标准分布、均匀分布等;但是，他们的输入参数都缺少或部分采用了您上面提到的三个输入参数中的一个，如果我处在您的位置，我会考虑三角分布，因为即使部分排除单个参数也会导致信息丢失。

import numpy as np
import matplotlib.pyplot as plt
h = plt.hist(np.random.triangular(-3, 0, 8, 100000), bins=200,
             density=True)
plt.show()

Numpy - Triangular Distribution

作为noted here：
与这些样本量一致的可能分布有无限多个。
但是，您可以引入其他假设来找到一些解决方案：

• 仅使用某些流行发行版的固定列表
• 对分布的参数添加约束

您可以将其视为优化问题：找到具有最佳拟合的分布及其参数（根据指定的最小值/最大值/平均值统计）。在伪代码中，解决方案如下所示：

candidates = []
for distribution in distributions:
    best_parameters, score = find_best_parameters(distribution, target_statistics)
    candidates.append((distribution, best_parameters, score))
best_distribution = sorted(candidates, key=lambda x: x[2])

使用此过程，您可以发现幂律分布可以产生与所需类似的统计信息：

s = stats.powerlaw(a=5.0909e-2, loc=1.00382, scale=4.122466e+2)
sample = s.rvs(size=100_000)
print(np.max(sample), np.min(sample), np.mean(sample))

最大值/最小值/平均值：

411.02946481216634 0.994030016 20.943683603008324

完整代码：

import numpy as np
from scipy import stats
import cma
from matplotlib import pyplot as plt

distributions_and_bounds = [
    (stats.cauchy, {'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.chi2, {'loc': [0, 1000], 'scale': [0, None]}),
    (stats.expon, {'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.exponpow, {'b': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.gamma, {'a': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.lognorm, {'s': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.norm, {'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.powerlaw, {'a': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.rayleigh, {'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.uniform, {'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.alpha, {'a': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.anglit, {'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.arcsine, {'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.burr, {'c': [0, None], 'd': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.argus, {'chi': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
    (stats.beta, {'a': [0, None], 'b': [0, None], 'loc': [-1000, 1000], 'scale': [0, None]}),
]

target_params = np.array([411, 1, 20.98])

candidates = []
for distribution, bounds in distributions_and_bounds:
    def objective(params):
        sample = distribution(*params).rvs(size=1_000)
        pred_params = np.array([np.max(sample), np.min(sample), np.mean(sample)])
        mse = (np.abs(target_params - pred_params)**2).mean()
        return mse

    x0 = np.ones(len(bounds))

    lower_bounds = [bound[0] for bound in bounds.values()]
    upper_bounds = [bound[1] for bound in bounds.values()]

    best_params, es = cma.fmin2(objective, x0, 1, {'bounds': [lower_bounds, upper_bounds]}, restarts=4)
    score = objective(best_params)
    candidates.append((score, distribution, best_params))

best_distribution = list(sorted(candidates, key=lambda x: x[0]))[0]
print(best_distribution)

这里，pycma包中的CMA-ES优化是为了简化。

让我们试试下面的函数：

import numpy as np
import random

def re_sample(min_v, max_v, mean_v, size):
    """
    Parameters
    ----------
    min_v  : Minimum value of the original population
    max_v  : Maximum value of the original population
    mean_v : Mean value of the original population
    size   : Number of observation we want to generate

    Returns
    -------

    sample : List of simulated values
    """

    s_min_to_mean=int(((max_v-mean_v)/(max_v-min_v))*size)
    sample_1=[random.uniform(min_v, mean_v) for i in range(s_min_to_mean)]
    sample_2=[random.uniform(mean_v, max_v) for i in range(size-s_min_to_mean)]

    sample=sample_1+sample_2

    sample=random.sample(sample, len(sample))

    sample=[round(x, 2) for x in sample] 

    return sample

当我测试这个功能时为：

sample = re_sample(1, 411, 20.98, 200)

print(np.mean(sample))
print(np.min(sample))
print(np.max(sample))
print(type(sample))
print(len(sample))
print(sample)

我得到以下输出：

>>> 19.8997
>>> 1.0
>>> 307.8
>>> <class 'list'>
>>> 200
>>> [20.55, 7.87, 3.48, 5.23, 18.54, 268.06, 1.71,....

快速编辑与细化（我后来意识到这一点）：您可以将平衡技巧应用于任何分布。
许多建议的解决方案的缺点是 * 使用浮点数达到MIN、MAX和AVERAGE的EXACT值的机会基本上为零 *。了解这一点意味着需要手动添加MIN和MAX值，但添加值会扰乱生成的分布。
一个简单的方法是，生成一个分布，将最小值和最大值相加，然后平衡它们以达到平均值：
1.设置最小值和最大值
1.计算平均值
1.添加点以补偿所需平均值的偏差（取决于MIN和MAX在所需平均值周围的不对称程度）
1.创建一个随机分布，在移动平均值后，该分布仍将在所需平均值和最近的边界条件之间拟合
1.将分布的均值转换为所需的真实均值
1.将生成的对称分布添加到点4之前的可用数据。
最初的3个步骤确保边界条件（最小值、最大值）不会扰乱平均值。步骤4-5创建一些数据，这些数据保证具有精确的所需平均值，并且将落在最小值和最大值之间。步骤6将数据组合到所需结果。

import math
import numpy as np

MAX, MIN, AVERAGE = [411, 3, 20.98]

data = [3, 411]

left = AVERAGE - MIN
right = MAX - AVERAGE
ratio = max(left, right)/min(left,right)

n = math.ceil(ratio) - 1
dx = math.ceil(ratio) - ratio  # this checks overcompensation due to working with integer numbers

data = data + [MIN]*(n) + [AVERAGE + left*dx]  # the second part compensates the overcompensation again :)

print(np.mean(data))
print(min(data))
print(max(data))

N = 1000

width = min(MAX-AVERAGE, AVERAGE-MIN)
print(width)

dist = np.random.normal(AVERAGE, width/3, N)

delta1 = np.mean(dist) - AVERAGE

dist = [x for x in dist if x > (MIN + delta1) and x < (MAX - delta1)]

delta2 = np.mean(dist) - AVERAGE
dist = [x - delta2 for x in dist]

full = data + dist

print(np.mean(full))
print(min(full))
print(max(full))

一个概率（函数）仅仅由它的最小值，平均值和最大值来定义是不够的。有（字面上）无限多的概率分布满足这些条件。
为了证明这一点，给出具有概率（最大-平均）/（最大-最小）的最小值和具有概率（平均-最小）/（最大-最小）的最大值的概率分布已经满足这些特征.
这一点很容易验证：

• 最小值和最大值是微不足道的。
• 平均值=最小值的概率 * 最小值+最大值的概率 * 最大值= {最小值 *（最大值-平均值）+最大值 （平均值-最小值）} /（最大值-最小值）=（-最小值 * 平均值+最大值 * 平均值）/（最大值-最小值）=（最大值-最小值） 平均值/（最大值-最小值）=平均值。

此外，正态分布是对称的 * 且 * 不限于观察值（例如：它没有最小值和最大值）。