Python频率分析：拟合良好但多分布出现无穷大AIC问题

在积雪年最大值序列的频率分析中，数据存在合法零值。计算AIC/BIC拟合优度时，零值对应的概率密度为0，其对数为负无穷，导致AIC/BIC结果为无穷大；即使将零值替换为极小非零值，该问题仍无法解决。

最小可复现代码：

import numpy as np
import lmoments3.distr as ld
import matplotlib.pyplot as plt

data=np.array([279, 244, 226, 216, 300, 208, 267, 239, 201, 277, 135, 79, 201, 15, 198, 231, 251, 277, 168,
      43, 53, 160, 366, 163, 239, 287, 196, 117, 0, 208])

data[data==0]=1E-07

return_periods = np.array([1.1111, 1.2500, 1.4286,2,3,5,10,20,50,100,200,1000])

# Fit GEV distribution to the data
params = ld.gev.lmom_fit(data)
# Calculate quantiles for the given return periods
quantiles = ld.gev.ppf(1 - 1 / return_periods, **params)

# Calculate Cunnane plotting positions for the data
n = len(data)
sorted_data = np.sort(data)
ranks = np.arange(1, n + 1)
empirical_cdf = (ranks - 0.4) / (n + 0.2)  # Cunnane plotting positions
data_return_periods = 1 / (1 - empirical_cdf)

# Plot the quantile plot
plt.figure(figsize=(8, 6))
plt.plot(return_periods, quantiles, marker='o', linestyle='-', color='b', label='GEV Quantiles')
plt.scatter(data_return_periods, sorted_data, color='r', label='Input Data (Cunnane)')
plt.xscale('log')
plt.yscale('linear')
plt.xlabel('Return Period (Years)', fontsize=12)
plt.ylabel('Quantile Value', fontsize=12)
plt.title('Quantile Plot for GEV Distribution with Input Data (Cunnane)', fontsize=14)
plt.grid(True, which="both", ls="--")
plt.legend()
plt.show()

n = len(data)
empirical_cdf = (np.arange(1, n + 1) - 0.4) / (n + 0.2)  # Cunnane plotting positions
return_periods = 1 / (1 - empirical_cdf)  # Convert CDF to return periods

logpdf_values = ld.gev.logpdf(sorted_data, **params)
pdf_values = ld.gev.pdf(sorted_data, **params)
log_likelihood = np.sum(logpdf_values)
print(f"log likelihood: {log_likelihood}")
k = len(params)  # Number of parameters in the model

aic = 2 * k - 2 * log_likelihood
bic = k * np.log(n) - 2 * log_likelihood

print(f"AIC: {aic}")
print(f"BIC: {bic}")

1. 零膨胀分布（Zero-Inflated Distributions） #

将数据建模为二元混合模型：一部分概率生成零值，另一部分概率生成服从GEV分布的非零值。似然函数由零值的概率对数和非零值的GEV对数密度组成，彻底避免负无穷问题。

实现步骤： #

• 统计数据中零值和非零值的数量
• 估计零值发生概率
• 对非零值拟合GEV分布并计算其对数似然
• 总对数似然为零值部分与非零值部分的对数似然之和
• 模型参数数为GEV的3个参数加上零概率参数，共4个

代码示例： #

# 还原原始数据，保留真实零值
data_original = np.array([279, 244, 226, 216, 300, 208, 267, 239, 201, 277, 135, 79, 201, 15, 198, 231, 251, 277, 168,
      43, 53, 160, 366, 163, 239, 287, 196, 117, 0, 208])

n_zero = np.sum(data_original == 0)
n_nonzero = len(data_original) - n_zero
nonzero_data = data_original[data_original != 0]

# 拟合非零数据的GEV分布
params_gev = ld.gev.lmom_fit(nonzero_data)
# 计算非零数据的对数似然
log_likelihood_nonzero = np.sum(ld.gev.logpdf(nonzero_data, **params_gev))

# 估计零概率
p_zero = n_zero / len(data_original)
# 总对数似然
total_log_likelihood = n_zero * np.log(p_zero) + n_nonzero * np.log(1 - p_zero) + log_likelihood_nonzero

# 计算AIC/BIC
k = 4
aic = 2 * k - 2 * total_log_likelihood
bic = k * np.log(len(data_original)) - 2 * total_log_likelihood

print(f"零膨胀GEV总对数似然: {total_log_likelihood}")
print(f"AIC: {aic}")
print(f"BIC: {bic}")

2. 截断分布修正（Truncated Distribution Adjustment） #

如果零值代表积雪深度的物理下限（实际值不可能小于0），可使用截断在0处的GEV分布。通过修正似然函数，将零值的概率纳入分布的累积密度中，非零值的对数密度则扣除截断带来的影响。

实现步骤： #

• 用含零值的数据拟合GEV分布
• 计算零值处的累积密度，作为截断修正项
• 分别计算零值和非零值的对数似然并求和
• 模型参数数为GEV本身的3个参数

代码示例： #

params = ld.gev.lmom_fit(data_original)  # 直接用含零值的数据拟合GEV
cdf_zero = ld.gev.cdf(0, **params)
trunc_correction = np.log(1 - cdf_zero)

# 分离零值和非零值
zero_mask = data_original == 0
nonzero_mask = ~zero_mask

log_likelihood_zero = np.sum(zero_mask) * np.log(cdf_zero)
log_likelihood_nonzero = np.sum(ld.gev.logpdf(data_original[nonzero_mask], **params)) - np.sum(nonzero_mask) * trunc_correction

total_log_likelihood = log_likelihood_zero + log_likelihood_nonzero

k = 3  # GEV本身的3个参数
aic = 2 * k - 2 * total_log_likelihood
bic = k * np.log(len(data_original)) - 2 * total_log_likelihood

print(f"截断GEV总对数似然: {total_log_likelihood}")
print(f"AIC: {aic}")
print(f"BIC: {bic}")

3. 似然函数的数值稳定修正 #

如果不想引入额外模型参数，可对极小值的对数似然做近似处理：将零值替换为接近0的阈值，用该阈值处的生存函数近似概率密度，避免直接计算pdf导致的数值下溢。

代码示例： #

epsilon = 1e-7
data_adjusted = data_original.copy()
data_adjusted[data_adjusted == 0] = epsilon

params = ld.gev.lmom_fit(data_adjusted)
cdf_zero = ld.gev.cdf(0, **params)
cdf_epsilon = ld.gev.cdf(epsilon, **params)

# 计算修正后的对数似然
logpdf_values = []
for x in data_adjusted:
    if x == epsilon:
        # 用生存函数的近似代替直接logpdf
        approx_logpdf = np.log((cdf_epsilon - cdf_zero)/epsilon)
        logpdf_values.append(approx_logpdf)
    else:
        logpdf_values.append(ld.gev.logpdf(x, **params))

log_likelihood = np.sum(logpdf_values)
k = 3
aic = 2 * k - 2 * log_likelihood
bic = k * np.log(len(data_adjusted)) - 2 * log_likelihood

print(f"数值修正后对数似然: {log_likelihood}")
print(f"AIC: {aic}")
print(f"BIC: {bic}")

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