python实现Theil-Sen Median斜率估计和Mann-Kendall趋势分析

python实现Theil-Sen Median斜率估计和Mann-Kendall趋势分析

Theil-Sen回归

Theil-Sen回归是一种鲁棒线性回归方法，用于减小异常值对拟合结果的影响。与最小二乘法和一些其他回归方法不同，Theil-Sen回归使用了一种称为中位数斜率的统计量来进行参数估计，从而提高回归模型的鲁棒性。

Theil-Sen回归的步骤如下：

• 对于给定的自变量和因变量数据，计算所有点对（两两数据点）的斜率。

• 然后找出所有斜率的中位数，这个中位数就是Theil-Sen回归的估计斜率。

• 使用估计的斜率，可以计算出拟合的截距。

Theil-Sen回归的优点在于它对于数据中的异常值有很好的鲁棒性。因为它使用的是中位数而不是均值来估计斜率，中位数对于异常值不敏感，能够较好地抵抗异常值的影响。这使得Theil-Sen回归在处理有较多异常值的数据时表现良好。

另外，Theil-Sen回归是一种非参数方法，它不需要对数据的分布进行假设，更加灵活。

代码实现

我的输入数据长这样，直接上代码

# -*- codeing = utf-8 -*-

import numpy as np
from scipy.stats import norm
import pandas as pd

path = r'F:\UHI\驱动力分析\python.xlsx'#文件路径
df = pd.read_excel(path, sheet_name='uhi')#表单


def sen_slope(x):
    """
    计算 Sen's 斜率估计值
    """
    n = len(x)
    slopes = []
    for i in range(n):
        for j in range(i + 1, n):
            slopes.append((x[j] - x[i]) / (j - i))
    return np.median(slopes)


def mann_kendall(x):
    """
    计算 Mann-Kendall 检验的统计量、Z 值和 p 值
    """
    n = len(x)
    s = 0
    for i in range(n - 1):
        for j in range(i + 1, n):
            s += np.sign(x[j] - x[i])
    var_s = n * (n - 1) * (2 * n + 5) / 18
    if s > 0:
        z = (s - 1) / np.sqrt(var_s)
    elif s < 0:
        z = (s + 1) / np.sqrt(var_s)
    else:
        z = 0
    p = 2 * (1 - norm.cdf(abs(z)))
    return s, z, p


for y in range(0, 45):  # 45行excel记录
    xi = list(df.loc[y].values)
    # print(xi)
    xii = xi[1:]
    # print(xii)
    x = xii

    slop = sen_slope(x)
    s, z, p = mann_kendall(x)

    print(slop, s, z, p)
    # print("Sen's slope:", slop)
    # print("Mann-Kendall statistic:", s)
    # print("Mann-Kendall检验的Z值：", z)
    # print("Mann-Kendall p-value:", p)

若出现计算数值过大的问题，可能是数据范围较广导致的数值溢出或精度问题，可通过数据归一化解决

for y in range(0, 45):  # 45行excel记录
    xi = list(df.loc[y].values)
    # print(xi)
    xii = xi[1:]
    # print(xii)
    x = xii
    # 数据归一化

    max_value = max(x)
    min_value = min(x)
    x_normalized = [(val - min_value) / (max_value - min_value) for val in x]

    slop = sen_slope(x_normalized)
    z, p = mann_kendall(x_normalized)

    print(slop)
    # print("Sen's slope:", slop)
    # print("Mann-Kendall statistic:", s)
    # print("Mann-Kendall p-value:", p)

感谢大家花时间来阅读本文 作者水平有限 有失误之处请大家斧正！