Tweedie分布是一种泊松分布和伽马分布的复合分布，它有三个参数： 一个是p，当p=1，Tweedie就是泊松分布，当p=2，Tweedie就是伽马分布。 第二个参数是μ，是Tweedie分布的期望。 第三个参数是ϕ，控制Tweedie分布的方Tweedie分布很适合用来描述某些场景的随机变量，比如商店中某样商品的日出货量；使用该分布来模拟随机变量的话，最大特点是能以一定的概率生成数值为0的样本 。 假设目前已知某样商品的日出货量可以用特定参数的Tweedie分布进行拟合，商店想进行库存优化，希望仓库储备的商品仅需满足特定概率P（比如，90%）情况下的安全库存，请编写一个函数，能够在给定p、μ、ϕ等分布参数的情况下，给出任意P的安全库存。

根据定义，Tweedie分布的概率密度函数为：

$$
f(y;\mu,\phi,p) = \frac{1}{Z(\mu,\phi,p)} \exp\left\{\frac{y\mu^{1-p}}{\phi(1-p)}\right\} \cdot \exp\left\{-\frac{\left(y\mu^{1-p}\right)^{2-p}}{\phi(2-p)}\right\}
$$

其中，$Z(\mu,\phi,p)$是归一化系数。为了求解安全库存，我们需要找到一个$y_0$，使得：

$$
P(Y \leq y_0) = \int_0^{y_0} f(y;\mu,\phi,p)dy = P
$$

这个积分一般是无法解析求解的，我们可以使用数值积分方法，比如Simpson法则或者Monte Carlo方法进行求解。下面是一个使用Simpson法则的求解代码实现：

import numpy as np
from scipy.integrate import quad

def inv_tweedie_cdf(p, mu, phi, p_val):
    def integrand(y, mu, phi, p):
        return np.exp((y * mu**(1-p)) / (phi * (1-p))) * np.exp(-(y * mu**(1-p))**(2-p) / (phi * (2-p)))
    
    def simpson_integral(f, a, b, n):
        h = (b-a) / n
        x = np.linspace(a, b, n+1)
        y = f(x, mu, phi, p)
        res = h/3 * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-1:2]))
        return res
    
    def f(y, mu, phi, p, p_val):
        return simpson_integral(integrand, 0, y, 1000) / simpson_integral(integrand, 0, np.inf, 1000) - p_val
    
    y0 = quad(f, 0, np.inf, args=(mu, phi, p, p_val))[0]
    return y0

使用时，只需要指定分布参数$p, \mu, \phi$以及所需概率$P$即可，函数将返回对应的安全库存$y_0$，例如：

p = 1.5
mu = 10
phi = 2
P = 0.9

y0 = inv_tweedie_cdf(p, mu, phi, P)
print(y0)

这个函数虽然使用了数值积分方法，但是由于Simpson法则的收敛速度非常快，所以在实际使用中效率还是非常高的。