_cdf方法本质上与泛型实现相同，因此其余部分假设它被删除。我还假设_pdf实现是正确的，我还没有看到任何其他建议。当mu s为0且sigma s为1时，它简化为柯西分布，即使参数改变，它似乎也会积分为1，PDF似乎与随机变量一致（见下文）。
使用quad集成PDF以获得CDF肯定会非常慢，但有一件事会让它更快，那就是在PDF中使用scipy.special.ndtr代替scipy.stats.norm.cdf。SciPy分布有很多开销，但特殊函数直接跳转到评估正态分布的CDF。

from scipy import stats, special
x = 1.5
%timeit special.ndtr(x)
# 1.43 µs ± 726 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit stats.norm.cdf(x)
# 107 µs ± 21.5 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

字符串
这将使您的分发的CDF加速10倍至20倍。
提高RVS的速度也很容易。_rvs的通用实现将CDF（数字）反转以获得PPF，并传入均匀分布的随机数以获得每个观察值（逆变换采样）。我们重写_rvs以仅生成两个正态随机变量并取比值。

def _rvs(self, mu_x, mu_y, sigma_x, sigma_y, size=None, random_state=None):
      x = random_state.normal(mu_x, sigma_x, size)
      y = random_state.normal(mu_y, sigma_y, size)
      return x/y

型
要检查此RVS实现与PDF之间的一致性，请执行以下操作：

rng = np.random.default_rng(8235834852452842)
nr_wiki_inst = normal_ratio_wiki(name='normal_ratio_wiki', seed=rng)

mu_x = 1
mu_y = 2
sigma_x = 3
sigma_y = 4
dist = nr_wiki_inst(mu_x, mu_y, sigma_x, sigma_y)

rvs = dist.rvs(size = 100_000)
plt.hist(rvs, density=True, bins=np.linspace(-10, 10, 300))

x = np.linspace(-10, 10, 300)
pdf = dist.pdf(x)
plt.plot(x, pdf)

型
x1c 0d1x的数据
这从分布中生成了100_000个随机变量（以8 ms为单位），使我们能够直观地将直方图与PDF进行比较。一致性很好，直方图的明显过量可以解释为直方图在范围[-10，10]上归一化，而不是整个真实的线。
你可以使用scipy.stats.sampling.NumericalInversePolynomial得到一个快速（但近似）的PPF方法。更容易描述NumericalInverseHermite的作用：在多个点上计算CDF，并使用这些数据插值PPF。NumericalInversePolynomial也可以使用插值，但更高级，它只需要PDF。

dist2 = sampling.NumericalInversePolynomial(dist, center=0, random_state=rng)
# show that `dist2`'s `ppf` method is the inverse of the original `cdf`
dist2.ppf(dist.cdf(1.5))  # 1.5000000001244032

型
我认为fit有问题，部分原因是它不知道sigma_x和sigma_y必须为正。为此，您需要定义_argcheck：

def _argcheck(self, mu_x, mu_y, sigma_x, sigma_y,):
        return (sigma_x > 0.) & (sigma_y > 0.)

型
但是还有其他问题，我想我现在不会调试它，相反，我们可以定义_shape_info，这样我们就可以在这个发行版中使用scipy.stats.fit函数，而不是发行版的fit方法。

# _ShapeInfo is private/not documented
# so don't use this for production code
from scipy.stats._distn_infrastructure import _ShapeInfo
...
    def _shape_info(self):
        # parameter name, discrete (True or False), domain, domain includes endpoints
        return [_ShapeInfo("mu_x", False, (-np.inf, np.inf), (False, False)),
                _ShapeInfo("mu_y", False, (-np.inf, np.inf), (False, False)),
                _ShapeInfo("sigma_x", False, (0, np.inf), (False, False)),
                _ShapeInfo("sigma_y", False, (0, np.inf), (False, False))]

型
现在可以使用fit函数：

# generate random variates to fit
rvs = dist.rvs(size=200)

# assign rough bounds to parameters
bounds = dict(mu_x=(-10, 10), mu_y=(-10, 10), sigma_x=(1e-1, 10), sigma_y=(1e-1, 10))
res = stats.fit(nr_wiki_inst, rvs, bounds=bounds)

# Show that negative log likelihood of fitted parameters is better
print(res.nllf((1, 2, 3, 4, 0, 1)))  # 454.43887891645863
print(res.nllf())  # 453.19908132025944
# The latter corresponds with the fitted parameters.
# It's lower, so it's a better fit according to MLE

型
一起来：

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, stats, special
from scipy.stats import sampling
from scipy.stats._distn_infrastructure import _ShapeInfo

class normal_ratio_wiki(stats.rv_continuous):
    def _pdf(self, z, mu_x, mu_y, sigma_x, sigma_y):
    
        a_z = np.sqrt(((1/(sigma_x**2))*(np.power(z,2))) + (1/(sigma_y**2)))
        b_z = ((mu_x/(sigma_x**2)) * z) + (mu_y/sigma_y**2)
        c = ((mu_x**2)/(sigma_x**2)) + ((mu_y**2)/(sigma_y**2))
        d_z = np.exp(((b_z**2)-((c*a_z**2))) / (2*(a_z**2)))
        pdf_z = ((b_z * d_z) / (a_z**3)) * (1/(np.sqrt(2*np.pi)*sigma_x*sigma_y)) * \
        (special.ndtr(b_z/a_z) - special.ndtr(-b_z/a_z)) + ((1/((a_z**2) * np.pi * sigma_x * sigma_y))*np.exp(-c/2))

        return pdf_z
        
    def _rvs(self, mu_x, mu_y, sigma_x, sigma_y, size=None, random_state=None):
        x = random_state.normal(mu_x, sigma_x, size)
        y = random_state.normal(mu_y, sigma_y, size)
        return x/y
    
    def _argcheck(self, mu_x, mu_y, sigma_x, sigma_y):
        return (sigma_x > 0) & (sigma_y > 0)

    def _shape_info(self):
        return [_ShapeInfo("mu_x", False, (-np.inf, np.inf), (False, False)),
                _ShapeInfo("mu_y", False, (-np.inf, np.inf), (False, False)),
                _ShapeInfo("sigma_x", False, (0, np.inf), (False, False)),
                _ShapeInfo("sigma_y", False, (0, np.inf), (False, False))]

rng = np.random.default_rng(11)
nr_wiki_inst = normal_ratio_wiki(name='normal_ratio_wiki', seed=rng)

mu_x = 1
mu_y = 2
sigma_x = 3
sigma_y = 4
dist = nr_wiki_inst(mu_x, mu_y, sigma_x, sigma_y)
rvs = dist.rvs(size = 100_000, random_state=rng)
x = np.linspace(-10, 10, 300)
pdf = dist.pdf(x)
plt.plot(x, pdf)
plt.hist(rvs, density=True, bins=np.linspace(-10, 10, 300))

dist2 = sampling.NumericalInversePolynomial(dist, center=0, random_state=rng)
dist2.ppf(dist.cdf(1.5))  # 1.5000000001244032

rvs = dist.rvs(size=200, random_state=rng)
bounds = dict(mu_x=(-10, 10), mu_y=(-10, 10), sigma_x=(1e-1, 10), sigma_y=(1e-1, 10))
res = stats.fit(nr_wiki_inst, rvs, bounds=bounds)
print(res.nllf((1, 2, 3, 4, 0, 1)))  # 399.2571952688255
print(res.nllf())  # 396.24324681778955

型