复合Sigmoid拟合导数曲线峰值统计代码修正请求

核心问题修正思路 #

• 针对标准差计算过大：采用加权标准差（以导数曲线值为权重，模拟概率密度）替代普通标准差，排除非峰值区间的干扰
• 针对无法区分两类峰值统计：先检测所有显著峰值，按峰值大小、x轴位置区分低温绝对最大值与高温局部最大值，再为每个峰值单独划定计算区间
• 针对标准差转换不准确：使用加权Z值计算，匹配偏态分布的特性
• 针对不符合偏态特性：所有统计量均采用导数曲线加权计算，拟合偏态分布的权重分布

修正后代码 #

import numpy as np
from scipy.signal import find_peaks
from scipy.integrate import trapz

def process_derivative_curve(x, dy_dx):
    # 1. 检测显著峰值，区分绝对最大值（低温区）与局部最大值（高温区）
    peaks, _ = find_peaks(dy_dx, height=np.percentile(dy_dx, 70))  # 筛选高置信度峰值
    peak_values = dy_dx[peaks]
    
    # 确定低温区绝对最大值（峰值最高的点）
    abs_peak_idx = np.argmax(peak_values)
    abs_peak_x = x[peaks[abs_peak_idx]]
    
    # 确定高温区局部最大值（剩余峰值中x最大的点）
    other_peaks = peaks[peaks != peaks[abs_peak_idx]]
    local_peak_x = x[other_peaks[np.argmax(x[other_peaks])]] if len(other_peaks) > 0 else None

    # 2. 为单个峰值划定计算区间（用左右谷底分割）
    def get_peak_interval(target_peak_x):
        peak_pos = np.where(x == target_peak_x)[0][0]
        # 找左侧谷底
        left_valley_idx = np.argmin(dy_dx[:peak_pos])
        # 找右侧谷底
        right_valley_idx = np.argmin(dy_dx[peak_pos:]) + peak_pos
        return x[left_valley_idx:right_valley_idx+1], dy_dx[left_valley_idx:right_valley_idx+1]

    results = {}
    # 处理低温区绝对最大值
    abs_x, abs_dy = get_peak_interval(abs_peak_x)
    # 归一化权重为概率密度
    abs_weights = abs_dy / trapz(abs_dy, abs_x)
    
    # 加权统计量计算
    abs_mean = np.average(abs_x, weights=abs_weights)
    # 加权中位数
    sorted_abs_idx = np.argsort(abs_x)
    sorted_abs_weights = abs_weights[sorted_abs_idx]
    cum_abs_weights = np.cumsum(sorted_abs_weights)
    abs_median = abs_x[sorted_abs_idx][np.argmax(cum_abs_weights >= 0.5)]
    # 加权标准差
    abs_var = np.average((abs_x - abs_mean)**2, weights=abs_weights)
    abs_std = np.sqrt(abs_var)
    # 均值与峰值距离的标准差等价（加权Z值）
    abs_dist_std = abs(abs_mean - abs_peak_x) / abs_std

    results["absolute_peak"] = {
        "peak_x": abs_peak_x,
        "mean": abs_mean,
        "median": abs_median,
        "std": abs_std,
        "distance_to_peak_std": abs_dist_std
    }

    # 处理高温区局部最大值（存在时）
    if local_peak_x is not None:
        local_x, local_dy = get_peak_interval(local_peak_x)
        local_weights = local_dy / trapz(local_dy, local_x)
        
        local_mean = np.average(local_x, weights=local_weights)
        # 加权中位数
        sorted_local_idx = np.argsort(local_x)
        sorted_local_weights = local_weights[sorted_local_idx]
        cum_local_weights = np.cumsum(sorted_local_weights)
        local_median = local_x[sorted_local_idx][np.argmax(cum_local_weights >= 0.5)]
        # 加权标准差
        local_var = np.average((local_x - local_mean)**2, weights=local_weights)
        local_std = np.sqrt(local_var)
        # 均值与峰值距离的标准差等价
        local_dist_std = abs(local_mean - local_peak_x) / local_std

        results["local_peak"] = {
            "peak_x": local_peak_x,
            "mean": local_mean,
            "median": local_median,
            "std": local_std,
            "distance_to_peak_std": local_dist_std
        }

    return results

# 示例调用（替换为你的6组数据集）
if __name__ == "__main__":
    # 模拟复合Sigmoid导数曲线数据
    x = np.linspace(0, 100, 1000)
    # 双Sigmoid导数叠加模拟目标曲线
    dy_dx = 0.3 * np.exp(-(x-20)**2/100) / (1 + np.exp(-(x-20)/5))**2 + 0.15 * np.exp(-(x-70)**2/80) / (1 + np.exp(-(x-70)/4))**2

    # 处理单组数据
    stats = process_derivative_curve(x, dy_dx)
    print("低温区绝对峰值统计：")
    for k, v in stats["absolute_peak"].items():
        print(f"{k}: {v:.2f}")
    if "local_peak" in stats:
        print("\n高温区局部峰值统计：")
        for k, v in stats["local_peak"].items():
            print(f"{k}: {v:.2f}")

    # 批量处理6组数据示例
    # for group_idx in range(6):
    #     # 加载第group_idx组的x和dy_dx数据
    #     x_group, dy_dx_group = load_your_dataset(group_idx)
    #     group_stats = process_derivative_curve(x_group, dy_dx_group)
    #     print(f"\n=== 第{group_idx+1}组数据统计 ===")
    #     # 输出统计结果

关键修正说明 #

• 峰值区间划分：通过寻找峰值左右的谷底作为区间边界，确保每个峰值仅对应自身的分布区间，避免跨区间干扰导致的标准差偏大
• 加权统计：将导数曲线值归一化为概率密度权重，所有统计量（均值、中位数、标准差）均采用加权计算，贴合偏态分布的特性
• 标准差转换：直接使用加权均值与峰值的距离除以加权标准差，得到等价的Z值，准确反映偏态分布下的相对距离
• 两类峰值区分：通过峰值大小筛选绝对最大值，再通过x轴位置筛选高温局部最大值，确保两类峰值的统计结果单独输出

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