将伽马分布参数拟合到期望和概率

Question

我的问题如下。我测量了一堆不同的物理特性，并将方法和测量的不确定性一直传播到某种效率比。对于我所有的物理特性，正态分布似乎是一个不错的选择，并且对于前几次计算和相应的传播，我的不确定性相当低。我将所有不确定性作为 kp=2 的扩展不确定性转发，这意味着覆盖了所有可能值的 95.45%。

但是，对于我现在执行的计算，为了获得效率，我最终得到了 (16+/-31)% 这样的结果，这是不可能的，因为我的效率只能从 0 扩展到 1。我的假设是，我找到了效率的“正确”期望值，但我的概率分布应该是正偏伽马分布而不是正态分布。假设区间 [0; 0.16+0.31], [0; 0.16+0.31*3/2] 和 [0;1] 涵盖了所有可能值的 95.45%、99.73% 和 100%，我应该能够解析地计算 gamma 分布参数 alpha 和 beta。不幸的是，以下使用 sympy 的代码不起作用，因为 sympy 无法处理它。

有人知道如何解决我的问题吗？

import sympy as sy
# from sympy import init_printing
# from sympy import symbols
# from sympy import Eq
# from sympy import var
# from sympy import integrate
# from sympy import gamma
# from sympy import Pow
# from sympy import exp

R_std3 = float(31 * 3/2 * 1/100)
R = float(16 * 1/100)
R_plus = R + R_std3

alpha = sy.symbols('alpha', positive = True)
beta = sy.symbols('beta', positive = True)
r = sy.symbols('r', positive = True)

z = sy.Pow(beta, alpha) * sy.Pow(r, alpha-1) * sy.exp(-beta*r)
n = sy.gamma(alpha)
pdf = sy.Pow(n, -1) * z

system = [sy.Eq(R, sy.Pow(beta, -1) * alpha ), 
          sy.Eq(0.9973, sy.integrate(pdf, (r,0,R_plus) ))
          ]

sy.solve( system )

Answer 1

我不希望第二个（超越）方程有一个很好的解析解，除非有一些聪明的恒等式可以用来简化它。不过，您可以用数字解决这个问题：

In [46]: system
Out[46]: 
⎡       α           α⋅γ(α, 0.625⋅β)⎤
⎢0.16 = ─, 0.9973 = ───────────────⎥
⎣       β               Γ(α + 1)   ⎦

In [47]: nsolve(system, [alpha, beta], [1, 10])
Out[47]: 
⎡2.16365408317279⎤
⎢                ⎥
⎣ 13.52283801983 ⎦