根据 Python 中的 P10、P50 和 P90 值生成概率密度函数

Question

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

上的模拟 (Monte Carlo?)

示例：我想模拟从 A 到 B 运行需要多长时间。

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

意味着 10% 的尝试 I 运行在 1 小时或更短时间内从 A 到 B，50% 的尝试 I 运行在 1.5 小时或更短时间内从 A 到 B（即 1.5是均值）。并且 10% 的尝试我将花费超过 2.5 小时。

谢谢

Answer 1

假设用 Myerson 分布对这个系统建模是合适的，然后，根据 Frontline Solvers，“[i]f 指定的百分位数是等距的（由下面的参数 b' 测量），那么 Myerson分布等同于正态分布。”你很幸运，遇到了一个简单的案例。

当然这不可能完全正确，因为法线有无限尾巴。您需要从左侧被截断的正常人群中抽取样本。

您需要的（未截断的）正态分布的平均值为 1.5 小时，并将其质量的 40% 放在 1 小时和 1.5 小时的平均值之间。标准正态将其质量的 40% 置于 -1.2815515655446004 和 0 之间。然后，给定标准正态随机偏差的供应，z 我们可以通过缩放它们 [=13] 将它们转换为所需的（未截断的）偏差=]，其中 0.5 是 1 小时到 1.5 小时之间的 'distance'，1.28155 是标准正态对应的 'distance'。

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

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

我会说这不太可能，因此不值得将其视为截断的正常值。

因此，要获取问题中定义的 Myerson 偏差样本，您可以这样做。

>>> from scipy.stats import norm
>>> sample = norm.rvs(loc=1.5, scale=0.5/1.28, size=100)

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

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

事实证明，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)

根据 Python 中的 P10、P50 和 P90 值生成概率密度函数

Generate a Probability Density Function based on P10, P50 and P90 values in Python

python

montecarlo

python-2.7

probability-density