如何在数据最优拟合分布图像上绘制RSS及RSS²/loc²函数

你需要先留存所有候选分布的RSS、loc参数，再通过双Y轴实现拟合结果与RSS的关联展示，以下是可直接运行的修改后代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy
import scipy.stats
from numpy.random import randn

# 加载模拟数据
L = np.linspace(5,1.0,210)
V = L + 0.02*randn(len(L))
data = np.exp(-3.0*L**2) + 0.1*randn(len(L))

### 统计数据直方图
y, x = np.histogram(data, bins=50, density=True)
# 计算每个bin的中点
x = (x + np.roll(x, -1))[:-1] / 2.0

plt.figure(figsize=(12,8))
plt.hist(data, bins=50, density=True)
plt.title("MOE")
plt.show()

# 候选分布列表
dist_names = ['norm', 'beta','gamma', 'pareto', 't', 'lognorm', 'invgamma',
              'invgauss',  'loggamma', 'alpha', 'chi', 'chi2']

sse = np.inf
sse_thr = 0.10
# 新增：存储所有分布的统计结果，用于后续绘制RSS²/loc²
dist_stats = []

# 遍历所有候选分布
for name in dist_names:
    dist = getattr(scipy.stats, name)
    param = dist.fit(data)
    # 提取参数
    loc = param[-2]
    scale = param[-1]
    arg = param[:-2]
    # 计算PDF
    pdf = dist.pdf(x, *arg, loc=loc, scale=scale)
    # 计算RSS（也就是代码里的SSE）
    model_sse = np.sum((y - pdf)**2)
    # 存储当前分布的统计值
    dist_stats.append({
        "name": name,
        "sse": model_sse,
        "loc": loc
    })
    # 更新最优分布
    if model_sse < sse :
        best_pdf = pdf
        sse = model_sse
        best_loc = loc
        best_scale = scale
        best_arg = arg
        best_name = name
    # 达到阈值提前退出
    if model_sse < sse_thr :
        break

# --------------------------
# 第一幅图：拟合结果结合RSS展示
# --------------------------
plt.figure(figsize=(12,8))
ax1 = plt.gca()
# 左Y轴：原始数据和拟合PDF
ax1.plot(x, y, label="原始数据分布", color="tab:blue")
ax1.plot(x, best_pdf, label=f"最优拟合：{best_name}", linewidth=3, color="tab:orange")
ax1.set_ylabel("概率密度", fontsize=12)
ax1.legend(loc="upper left")
# 右Y轴：每个bin的残差平方
ax2 = ax1.twinx()
residual_square = (y - best_pdf)**2
ax2.bar(x, residual_square, width=np.diff(x)[0]*0.8, alpha=0.3, color="tab:red", label="单bin残差平方")
ax2.set_ylabel("残差平方值", fontsize=12)
# 标注全局总RSS
plt.text(0.65, 0.9, f"总RSS = {sse:.4f}", transform=ax1.transAxes, 
         bbox=dict(facecolor="white", alpha=0.8), fontsize=12)
# 合并图例
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc="upper right")
plt.title("拟合分布与RSS关联图", fontsize=14)
plt.show()

# --------------------------
# 第二幅图：RSS²/loc² 图像
# --------------------------
plt.figure(figsize=(12,6))
# 提取计算指标，加极小值避免loc为0时报错
dist_names_list = [d["name"] for d in dist_stats]
metric = [ (d["sse"]**2) / (d["loc"]**2 + 1e-8) for d in dist_stats ]
plt.bar(dist_names_list, metric, color="tab:green", alpha=0.7)
plt.xticks(rotation=45, fontsize=10)
plt.ylabel("RSS² / loc²", fontsize=12)
plt.xlabel("候选分布名称", fontsize=12)
plt.title("不同候选分布的RSS²/loc²指标对比", fontsize=14)
# 给每个柱子加数值标签
for i, v in enumerate(metric):
    plt.text(i, v+max(metric)*0.01, f"{v:.2e}", ha="center", fontsize=9)
plt.tight_layout()
plt.show()

# 打印最优分布信息
print("Selected Model : ", best_name)
print("Loc. param. : ", best_loc)
print("Scale param. : ", best_scale)
print("Other arguments : ", best_arg)
print("SSE : ", sse)

修改说明 #

• 新增了dist_stats列表存储所有候选分布的名称、RSS、loc参数，解决了原代码仅留存最优值无法绘制多分布对比指标的问题
• 拟合结果图采用双Y轴设计：左轴展示原始数据分布和最优拟合PDF，右轴用柱状图展示每个bin对应的残差平方，全局总RSS直接标注在图内
• 新增的RSS²/loc²图采用柱状图展示所有候选分布的对应指标，加了数值标签方便对比，同时添加了极小值1e-8避免loc为0时的除以0报错

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