My problem is the following. I measured a bunch of different physical properties and propagated the methodic and measurement uncertainties all the way to some kind of efficiency ratio. For all my physical properties, a normal distribution seemed to be a good choice, and for the first few calculations and corresponding propagations I had fairly low uncertainties. I forwarded all uncertainties as an extended uncertainty of kp=2 meaning coverage of 95.45% of all possible values.
However, for the calculations I am performing right now in order to get the efficiency I end up with results like (16+/-31)%, which is not possible since my efficiency can only spread from 0 to 1. My assumption would be, that I found the "correct" expectation value for the efficiency, but my probability distribution should be a positive skewed gamma-distribution rather than a normal-distribution. With the assumption that the intervals [0; 0.16+0.31], [0; 0.16+0.31*3/2] and [0;1] cover 95.45%, 99.73% and 100% of all possible values I should be able to calculate the gamma-distribution parameters alpha and beta analytically. Unfortunately, the following code using sympy doesn't work because sympy can't handle it.
Does someone have an idea how to solve my problem?
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 )
I wouldn't expect the second (transcendental) equation to have a nice analytic solution unless there's some clever identity that you can use to simplify it. You can solve this numerically though:
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 ⎦