基于Scipy实现独立正态分布比分布遇收敛与迭代失败问题求助

先修正现有代码的核心错误 #

你的_cdf方法存在逻辑错误：当前积分区间是(-inf, inf)，返回值始终为1，完全不符合CDF的定义（CDF应为从负无穷到当前z的积分）。先修正这个基础问题：

def _cdf(self, z, mu_x, mu_y, sigma_x, sigma_y):
    # CDF是从负无穷到z的积分，而非到正无穷
    cdf_z, _ = integrate.quad(self._pdf, -np.inf, z, args=(mu_x, mu_y, sigma_x, sigma_y))
    return cdf_z

即使修正CDF，默认的rvs和fit效率依然极低，且容易出现收敛问题。以下是针对性的自定义实现：

1. 自定义_rvs方法：直接抽样法 #

由于z = x/y（x~N(mu_x, sigma_x²)，y~N(mu_y, sigma_y²)），最高效的抽样方式是先生成x和y的样本再做除法，完全避开数值积分的低效问题：

def _rvs(self, mu_x, mu_y, sigma_x, sigma_y, size=1, random_state=None):
    rng = random_state or np.random.default_rng()
    x_samples = rng.normal(mu_x, sigma_x, size=size)
    y_samples = rng.normal(mu_y, sigma_y, size=size)
    # 避免除以0，给y添加极小扰动
    y_samples = np.where(y_samples == 0, 1e-10, y_samples)
    return x_samples / y_samples

2. 自定义_ppf方法：数值求逆CDF #

如果需要分位数函数（PPF），可以基于修正后的CDF用二分法数值求解，避免默认方法的收敛问题：

def _ppf(self, q, mu_x, mu_y, sigma_x, sigma_y):
    # 先通过抽样确定z的大致范围，缩小二分法区间
    rvs_sample = self._rvs(mu_x, mu_y, sigma_x, sigma_y, size=1000)
    z_min, z_max = np.percentile(rvs_sample, [0.1, 99.9])
    
    # 二分法求解单个分位数
    def find_z(target_q):
        low, high = z_min, z_max
        for _ in range(100):
            mid = (low + high) / 2
            current_cdf = self._cdf(mid, mu_x, mu_y, sigma_x, sigma_y)
            if current_cdf < target_q:
                low = mid
            else:
                high = mid
        return (low + high) / 2
    
    # 处理批量q值
    if np.isscalar(q):
        return find_z(q)
    else:
        return np.array([find_z(q_i) for q_i in q])

3. 自定义_fit方法：极大似然估计+合理初始值 #

默认拟合方法不适合这种复合分布，我们基于极大似然估计，结合样本矩给出合理初始值，解决收敛问题：

def _fit(self, data, **kwargs):
    # 定义负对数似然函数，用于极小化
    def neg_log_likelihood(params):
        mu_x, mu_y, sigma_x, sigma_y = params
        # 避免标准差为负
        if sigma_x <=0 or sigma_y <=0:
            return np.inf
        pdf_vals = self._pdf(data, mu_x, mu_y, sigma_x, sigma_y)
        # 避免log(0)的情况
        pdf_vals = np.where(pdf_vals < 1e-20, 1e-20, pdf_vals)
        return -np.sum(np.log(pdf_vals))
    
    # 基于样本矩给出初始值猜测
    mu_z = np.mean(data)
    sigma_z = np.std(data)
    init_mu_y = kwargs.get('init_mu_y', 300)
    init_mu_x = mu_z * init_mu_y
    init_sigma_y = kwargs.get('init_sigma_y', 60)
    init_sigma_x = np.sqrt( sigma_z**2 * init_mu_y**2 + (init_mu_x**2 * init_sigma_y**2) / init_mu_y**2 )
    # 确保初始标准差为正
    init_sigma_x = max(init_sigma_x, 10)
    
    # 调用优化器求解
    from scipy.optimize import minimize
    result = minimize(neg_log_likelihood, x0=[init_mu_x, init_mu_y, init_sigma_x, init_sigma_y],
                      bounds=[(None, None), (None, None), (1e-5, None), (1e-5, None)])
    
    return result.x

完整修正后的代码 #

import numpy as np
from scipy.stats import norm, rv_continuous
import scipy.integrate as integrate
from scipy.optimize import minimize

class normal_ratio_wiki(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)) * \
        (norm.cdf(b_z/a_z) - norm.cdf(-b_z/a_z)) + ((1/((a_z**2) * np.pi * sigma_x * sigma_y))*np.exp(-c/2))

        return pdf_z

    def _cdf(self, z, mu_x, mu_y, sigma_x, sigma_y):
        cdf_z, _ = integrate.quad(self._pdf, -np.inf, z, args=(mu_x, mu_y, sigma_x, sigma_y))
        return cdf_z
    
    def _rvs(self, mu_x, mu_y, sigma_x, sigma_y, size=1, random_state=None):
        rng = random_state or np.random.default_rng()
        x_samples = rng.normal(mu_x, sigma_x, size=size)
        y_samples = rng.normal(mu_y, sigma_y, size=size)
        y_samples = np.where(y_samples == 0, 1e-10, y_samples)
        return x_samples / y_samples
    
    def _ppf(self, q, mu_x, mu_y, sigma_x, sigma_y):
        rvs_sample = self._rvs(mu_x, mu_y, sigma_x, sigma_y, size=1000)
        z_min, z_max = np.percentile(rvs_sample, [0.1, 99.9])
        
        def find_z(target_q):
            low, high = z_min, z_max
            for _ in range(100):
                mid = (low + high) / 2
                current_cdf = self._cdf(mid, mu_x, mu_y, sigma_x, sigma_y)
                if current_cdf < target_q:
                    low = mid
                else:
                    high = mid
            return (low + high) / 2
        
        if np.isscalar(q):
            return find_z(q)
        else:
            return np.array([find_z(q_i) for q_i in q])
    
    def _fit(self, data, **kwargs):
        def neg_log_likelihood(params):
            mu_x, mu_y, sigma_x, sigma_y = params
            if sigma_x <=0 or sigma_y <=0:
                return np.inf
            pdf_vals = self._pdf(data, mu_x, mu_y, sigma_x, sigma_y)
            pdf_vals = np.where(pdf_vals < 1e-20, 1e-20, pdf_vals)
            return -np.sum(np.log(pdf_vals))
        
        mu_z = np.mean(data)
        sigma_z = np.std(data)
        init_mu_y = kwargs.get('init_mu_y', np.mean(y) if 'y' in kwargs else 300)
        init_mu_x = mu_z * init_mu_y
        init_sigma_y = kwargs.get('init_sigma_y', np.std(y) if 'y' in kwargs else 60)
        init_sigma_x = np.sqrt( sigma_z**2 * init_mu_y**2 + (init_mu_x**2 * init_sigma_y**2) / init_mu_y**2 )
        init_sigma_x = max(init_sigma_x, 10)
        
        result = minimize(neg_log_likelihood, x0=[init_mu_x, init_mu_y, init_sigma_x, init_sigma_y],
                          bounds=[(None, None), (None, None), (1e-5, None), (1e-5, None)])
        
        return result.x

# 测试代码
rng1 = np.random.default_rng(99)
rng2 = np.random.default_rng(88)

# Sample Data 1
x = rng1.normal(141739.951, 1.223808e+06, 1000)
y = rng2.normal(333.91, 64.494571, 1000)
z = x / y

mu_x = x.mean()
mu_y = y.mean()
sigma_x = x.std()
sigma_y = y.std()

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

# 测试rvs
nr_wiki_vars = nr_wiki_inst.rvs(mu_x, mu_y, sigma_x, sigma_y, size=100)
print("rvs样本生成成功，前5个样本：", nr_wiki_vars[:5])

# 测试fit
nr_wiki_params = nr_wiki_inst.fit(z, init_mu_y=mu_y, init_sigma_y=sigma_y)
print("拟合参数：", nr_wiki_params)

关键说明 #

• _rvs直接基于原始正态分布抽样，效率远高于数值积分逆变换，彻底解决收敛问题。
• _ppf用二分法结合抽样确定的区间，避免了默认方法的发散问题。
• _fit通过极大似然估计+合理初始值，适配了复合分布的拟合需求，解决了默认方法的不适用性。

内容的提问来源于stack exchange，提问作者Smith