在Python中基于P10、P50和P90值生成一个概率密度函数

Question

我想使用P10、P50和P90值作为输入： A） 生成一个概率密度函数（这感觉像一个Myerson分布，但我不知道如何在Python中实现这一点。Excel中有一个插件可以做到这一点；SIPMath） B） 运行模拟（蒙特卡洛？）在PDF上

我想模拟一下从a到B需要多长时间

P10 = 1 hour
P50 = 1.5 hours
P90 = 2.5 hours

也就是说，我在1小时或更短时间内从A运行到B的尝试次数的10%，在1.5小时或更短时间内从A运行到B的尝试次数的50%（即1.5是平均值）。10%的尝试我会花2.5个小时以上。在

谢谢你

Answer 1

事实证明，metalog是适合在这种情况下使用的发行版。一个非常灵活的分布，可以处理4种不同的情况：无边界、下限（具有最小值）、上限（最大值）和有界（最小值和最大值）。在

def metalog_multi(p10, p50, p90, numberofsamples, p0 = None, p100 = None):
    p10 = float(p10)
    p50 = float(p50)
    p90 = float(p90)
    if p0 != None:
        p0 = float(p0)
    if p100 != None:
        p100 = float(p100)


    samples = []
    for i in range(numberofsamples):
        x = random.random()
        if p0 == None and p100 == None:
            # unbound

            sample = p50 + 0.5 * (log((1 - 0.1) / 0.1)) ** (-1) * (p90 - p10) * log(x / (1 - x)) + ((1 - 2 * 0.1) * (log((1 - 0.1) / 0.1))) ** -1 * (1 - 2 * (p50 - p10) / (p90 - p10)) * (p90 - p10) * (x - 0.5) * log(x / (1 - x))

        elif p100 == None:
            # lower bound
            sample = p0 + e ** (log(p50 - p0) + 0.5 * (log((1 - 0.1) / 0.1)) ** -1 * log((p90 - p0) / (p10 - p0)) * log(x / (1 - x)) + ((1 - 2 * 0.1) * (log((1 - 0.1) / 0.1))) ** -1 * log(((p90 - p0) * (p10 - p0)) / (p50 - p0) ** 2) * (x - 0.5) * log(x / (1 - x)))
        elif p0 == None:
            # upper bound
            sample = p100 - e ** (-(-log(p100 - p50) - (0.5) * (log((1 - 0.1) / 0.1)) ** -1 * log((p100 - p90) / (p100 - p10)) * log(x / (1 - x)) - ((1 - 2 * 0.1) * (log((1 - 0.1) / 0.1))) ** -1 * log(((p100 - p90) * (p100 - p10)) / (p100 - p50) ** 2) * (x - 0.5) * log(x / (1 - x))))
        else:
            # bound
            sample = (p0 + p100 * e ** (log((p50 - p0) / (p100 - p50)) + (0.5) * (log((1 - 0.1) / 0.1)) ** -1 * log(((p90 - p0) / (p100 - p90)) / ((p10 - p0) / (p100 - p10))) * log(x / (1 - x)) + ((1 - 2 * 0.1) * (log((1 - 0.1) / 0.1))) ** -1 * log((((p90 - p0) / (p100 - p90)) * ((p10 - p0) / (p100 - p10))) / ((p50 - p0) / (p100 - p50)) ** 2) * (x - 0.5) * log(x / (1 - x)))) / (1 + e ** (log((p50 - p0) / (p100 - p50)) + (0.5) * (log((1 - 0.1) / 0.1)) ** -1 * log(((p90 - p0) / (p100 - p90)) / ((p10 - p0) / (p100 - p10))) * log(x / (1 - x)) + ((1 - 2 * 0.1) *(log((1 - 0.1) / 0.1))) ** -1 * log((((p90 - p0) / (p100 - p90)) * ((p10 - p0) / (p100 - p10))) / ((p50 - p0) / (p100 - p50)) ** 2) * (x - 0.5) * log(x / (1 - x))))
        samples.append(sample)
    return samples


p0_in = 10
p10_in = 20
p50_in = 40
p90_in = 80
p100_in = 250
numberofsamples = 10000
data = metalog_multi(p10_in, p50_in, p90_in, numberofsamples, p0 = p0_in)

p10_out = np.percentile(data,10)
p50_out = np.percentile(data,50)
p90_out = np.percentile(data,90)

Answer 2

这是我试图解决的问题。在b'=1的情况下，数据是对称的，我们应该将其视为正态分布。随着样本数量的增加，pX_out接近pX_in。我本希望能够设置上下两个障碍，但我还没有想出如何做到这一点。如有任何建议，我们将不胜感激。谢谢您。在

def myerson(p10, p50, p90, number_of_samples):
    b_mark = ((float(p90) - float(p50)) / (float(p50) - float(p10)))
    samples = []
    for i in range(number_of_samples):
        rand_numb = random.random()
        factor = norm.ppf(rand_numb, 0, 0.780304146072379)
        if 0.9999 < b_mark < 1.0001: 
            sample = p50 + (p90 - p50) * factor
        else:
            sample = p50 + (p90 - p50)*((b_mark**factor - 1)/(b_mark - 1))
        samples.append(sample)
    return samples

p10_in = 90
p50_in = 100
p90_in = 111
numberofsamples = 10000
data = myerson(p10_in, p50_in, p90_in, numberofsamples)

p10_out = np.percentile(data,10)
p50_out = np.percentile(data,50)
p90_out = np.percentile(data,90)

Answer 3

假设用Myerson分布来建模这个系统是合适的，那么根据Frontline Solvers，“[i]如果指定的百分位是等距的（由下面的参数b'来测量），那么Myerson分布就等价于正态分布。”你有一个简单的例子。在

当然，这不一定是真的，因为法线有无限的尾巴。您需要从左侧截断的正常总体中提取样本。在

你需要的（不可信的）正态分布的平均值是1.5小时，它的质量的40%在1小时和1.5小时的平均值之间。标准法线将其质量的40%置于-1.281551565544604和0之间。然后，给定一组标准正态随机偏差，z我们可以通过缩放它们来将它们转换为所需的（未加密的）偏差，0.5*(z+1.5)/1.28155，其中0.5是1小时到1.5小时之间的“距离”，1.28155是标准法线的相应“距离”。在

作为正态分布，一些小于零的随机变量可能产生。但是，使用scipy库我发现

>>> norm.cdf(0, loc=1.5, scale=0.5/1.28)
6.151715518325519e-05

我想说这是不太可能的，所以不值得费心去把它当作一个截断的正常值。在

因此为了获得您问题中定义的Myerson偏差样本，您可以这样做。在

^{pr2}$

loc和{}的值如我们所讨论的。size的值将是您需要的任何样本大小。在