如何用Python的scipy和lmoments3拟合Log Pearson3分布？

问题描述 #

我正在尝试将Log-Pearson III（LP3）分布拟合到年最大流量（AMF）径流数据中，但遇到了工具链的问题：scipy和lmoments3库仅提供Pearson III分布的拟合函数，没有直接的LP3实现。我尝试对AMF取对数后拟合Pearson III再取反对数，以为这等效于直接拟合LP3，但实际参数估计存在差异，也找不到明确的指导文档。

现有代码 #

import os
import numpy as np
import pandas as pd
import scipy.stats as st
import lmoments3 as lm
from lmoments3 import distr, stats

# Define the folder path and get the list of CSV files
folder_path = '.../daily_ts/'
# folder_path = '.../all_sites/'

csv_files = sorted([f for f in os.listdir(folder_path) if f.endswith('.csv')])

# Iterate over each file and process the data
for i, file in enumerate(csv_files[5:6]):
    station_code = file.split('_')[0]
    # Read the data
    df = pd.read_csv(os.path.join(folder_path, file), skiprows=27, names=["Date", "Flow (ML)", "Bureau QCode"])
    station_code = file.split('_')[0]

    # Process the data
    df = df.dropna(subset=["Flow (ML)"]) # Ensure 'Flow (ML)' has no NaN values
    df['Date'] = pd.to_datetime(df['Date'])  # Ensure the 'Date' column is in datetime format
    df['Year'] = df['Date'].dt.year  # Extract the year from the 'Date' column and add it as a new column
    df_max_flow = df.loc[df.groupby('Year')['Flow (ML)'].idxmax()]  # Max daily per year
    df_max_flow_sorted = df_max_flow.sort_values(by='Flow (ML)', ascending=False)  # Sort by 'Flow (ML)' in descending order
    df_max_flow_sorted['Rank'] = range(1, len(df_max_flow_sorted) + 1)
    df_max_flow_sorted['ARI_Empir'] = (df_max_flow_sorted['Rank'].iloc[-1] + 1) / df_max_flow_sorted['Rank']
    df_max_flow_sorted['AEP_Empir'] = 1 / df_max_flow_sorted['ARI_Empir']
    data = df_max_flow_sorted['Flow (ML)'].values

# Fit Log-Pearson3 using Maximun Likelihood Estimation method (MLE)
dist_name = 'pearson3'
# Use getattr to dynamically get the fitting method command
scipy_dist_fit = getattr(st, dist_name)
# Calculate natural log of the data
log_data = np.log(data)
param = scipy_dist_fit.fit(log_data)

# Applying the Kolmogorov-Smirnov test
ks_stat, p_val = st.kstest(log_data, dist_name, args=param)

print(f"parameters: {param}")
print(f"ks_stat, p_val: {ks_stat, p_val}")

ARI_dict = {}
AEP_dict = {}
data_dict = {}
logdata_dict = {}

# test the results of ARI and AEP according to the fitted distribution
loc = param[1]
scale = param[2]
shape = param[0]
# get the attribute to run SciPy stat in the loop
scipy_dist_fit = getattr(st, dist_name)
# run the SciPy stat
fitted_dist = scipy_dist_fit(shape, loc=loc, scale=scale)

# Calculate the return period and AEP based on the fitted_dist and add relevant columns to the table
AEP = 1 - fitted_dist.cdf(log_data)
ARI = 1 / AEP

# Store AEP and ARI in a dictionary
AEP_dict['AEP_LP3_MLE'] = np.sort(AEP)
ARI_dict['ARI_LP3_MLE'] = np.sort(ARI)[::-1]
data_dict['data'] = data
logdata_dict['log_data'] = log_data

AEP_df = pd.DataFrame(AEP_dict)
ARI_df = pd.DataFrame(ARI_dict)
data_df = pd.DataFrame(data_dict)
logdata_df = pd.DataFrame(logdata_dict)

# Reset indices of all DataFrames to ensure they align correctly
AEP_df = AEP_df.reset_index(drop=True)
ARI_df = ARI_df.reset_index(drop=True)
df_max_flow_sorted = df_max_flow_sorted.reset_index(drop=True)
data_df = data_df.reset_index(drop=True)
logdata_df = logdata_df.reset_index(drop=True)

df_fitDist = pd.concat([df_max_flow_sorted, AEP_df, ARI_df], axis=1)
LP3_FFA = pd.concat([data_df, logdata_df, AEP_df, ARI_df], axis=1)

# data = fitted_dist.ppf(1 - AEP)  # Inverse of CDF
flood_100year = np.exp(fitted_dist.ppf(1 - 0.01))
# flood_100year =fitted_dist.ppf(1 - 0.01)
print(f"100-Year Flood = {flood_100year}")
flood_50year = np.exp(fitted_dist.ppf(1 - 0.02))
# flood_50year = fitted_dist.ppf(1 - 0.02)
print(f"50-Year Flood = {flood_50year}")

核心问题解析 #

你采用的"取对数拟合Pearson III再指数还原"的逻辑本质就是Log-Pearson III分布的标准拟合流程，所谓的参数估计差异主要来自以下几点：

1. 参数定义混淆：LP3没有独立的参数体系，它的参数完全对应对数后Pearson III的参数，但原始数据的统计量（均值、方差、偏度）与对数后数据的统计量是通过对数/指数变换关联的，并非直接一一对应。
2. 估计方法差异：scipy的pearson3.fit()默认用极大似然估计（MLE），而水文领域更常用L-moment方法（对极端值更稳健），两种方法的参数结果自然会有区别。
3. 验证逻辑偏差：你计算的经验ARI和拟合分布的ARI在排序后拼接，可能导致数据对应关系混乱，看起来差异被放大。

修正方案 #

方法1：优化Scipy的MLE拟合流程 #

你的代码逻辑正确，只需补充细节优化：

• 数据预处理：确保径流数据无0或负值（取对数前检查，可添加data = data[data > 0]过滤异常）
• 验证逻辑修正：不要单独排序AEP/ARI后拼接，直接对原始数据计算对应的拟合AEP/ARI，保持数据对应关系
• 尾部验证：KS检验对分布尾部不敏感，建议增加QQ图验证对数后数据与Pearson III分布的贴合度

优化后的关键代码片段：

# 过滤异常值（确保无0/负值）
data = data[data > 0]
log_data = np.log(data)

# 拟合对数后的Pearson III
param = st.pearson3.fit(log_data)
fitted_dist = st.pearson3(*param)

# 直接对原始数据计算拟合的AEP/ARI，保持索引对应
df_max_flow_sorted['AEP_LP3_MLE'] = 1 - fitted_dist.cdf(log_data)
df_max_flow_sorted['ARI_LP3_MLE'] = 1 / df_max_flow_sorted['AEP_LP3_MLE']

# 绘制QQ图验证拟合效果
import matplotlib.pyplot as plt
st.probplot(log_data, dist=fitted_dist, plot=plt)
plt.title("Q-Q Plot of Log-Transformed AMF vs Pearson III")
plt.show()

方法2：基于lmoments3的L-moment拟合（水文领域推荐） #

L-moment方法对极端值的稳健性优于MLE，更适合径流频率分析，步骤与MLE一致：

# 过滤异常值
data = data[data > 0]
log_data = np.log(data)

# 计算对数数据的L-moments
lmoms = lm.lmom_ratios(log_data)
# 用L-moment拟合Pearson III分布
pearson3_params = distr.pearson3.lmom_fit(lmoms)
pearson3_dist = distr.pearson3(**pearson3_params)

# 计算设计洪水
flood_100year_lmom = np.exp(pearson3_dist.ppf(1 - 0.01))
flood_50year_lmom = np.exp(pearson3_dist.ppf(1 - 0.02))
print(f"100-Year Flood (L-moment): {flood_100year_lmom}")
print(f"50-Year Flood (L-moment): {flood_50year_lmom}")

# 计算拟合的AEP/ARI并入DataFrame
df_max_flow_sorted['AEP_LP3_LMOM'] = 1 - pearson3_dist.cdf(log_data)
df_max_flow_sorted['ARI_LP3_LMOM'] = 1 / df_max_flow_sorted['AEP_LP3_LMOM']

验证建议 #

• QQ图分析：观察对数后数据的经验分位数与Pearson III理论分位数的贴合程度，重点关注尾部（右上角）是否拟合良好
• 返回期曲线对比：绘制经验ARI（y轴）与流量（x轴）的散点图，叠加LP3拟合的曲线，验证极端流量（高ARI）的拟合效果
• 参数合理性检查：对数后Pearson III的偏度一般在±2之间，如果偏度绝对值过大，需检查数据是否包含异常值（如特大洪水未单独处理）

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