Python性能优化实践—曼肯德尔趋势检验计算效率提升
曼—肯德尔算法可以进行序列趋势和突变点的检验，实践发现对趋势的检验有较高的准确性，对突变点的检验相对较差。本文参考相关资料通过Python实现了对序列的曼肯德尔检验并进行了运行优化，提升了生产中的应用价值。
本次实践主要使用了算法优化、函数缓存、循环遍历优化、多进程、numba优化。由于整个过程以numpy数组的数值计算为主进一步还可以使用向量化进行加速，或者使用numexpr对Z值计算部分进行优化，但是未实践。
本次优化前将使用第三方库的算法编写函数，并通过apply函数将函数应用到每一列，在我的电脑上每10000行上的计算时间为70秒以上。
方法一：优化循环后为40多秒可以运算完毕，在使用常量法优化然后利用三核多进程计算可以18秒运算完毕。
方法二：利用函数缓存法优化正态分布的分位点部分计算，然后使用numba加速，预热后运算可以达到4.5秒。
主要参考资料：
趋势检验方法（二）MK趋势检验_aaakirito的博客-CSDN博客_mk检验
程序的优化的主要步骤为发现性能瓶颈和对主要瓶颈进行优化，主要方法有算法优化、循环优化、计算优化和编译优化。
注意：由于算法的原因，下面例子中的numpy数组必须是一维行数组，即reshape(1,-1)。
文章目录
- 一、程序性能分析
- 二、程序性能优化
- 1、算法优化
- 第一部分：计算S值，趋势检验方法
- 第二部分：优化scipy.stats 模块中标准正态分布的分位点函数norm.ppf()
- 第三部分:计算sen斜率
- 使用函数缓存进行优化
- 常量法
- 2、循环优化
- 3、计算优化
- 4、编译优化
- 三、整体代码
一、程序性能分析
在jupyter中
粗分析：
使用%%timeit对cell进行分析，比较各种实现方法的效率
细分析：
使用line_profiler对函数逐行分析，找到耗时占比最大的行
%load_ext line_profiler 加载line_profiler功能
%lprun -f funcname funcname(argu) 对函数funcname进行分析,获取每行的运行时间
二、程序性能优化
1、算法优化
将使用第三方库的算法，改为手动优化过的算法
进行曼—肯德尔趋势检验可以使用现成的pymannkendall进行计算，这个库中包含了曼—肯德尔趋势检验和他的多种扩展，计算量也较为全面。使用较为方便但是运行较慢。
import pymannkendall as mk
def mk_ku(data):
result=mk.original_test(data)
return (result.trend,result.z,result.slope)
参考CSDN上aaakirito的博客_CSDN博客-ACM算法题,ACM简单题,python领域博主的方法进行改进，提高算法的运行效率，但是功能缩减仅能计算最常使用的trend、z和solope三种结果。
def mk_(y): alpha=0.1 n = len(y) # 计算S的值v = np.array(y[1:]).reshape((-1, 1))v = np.tile(v, len(y) - 1)for i in range(len(y) - 1):v[:, i] -= y[i]v[0:i, i] = 0v = np.sign(v)s = np.sum(v)# 判断x里面是否存在重复的数，输出唯一数队列unique_x,重复数数量队列tpunique_x, tp = np.unique(x, return_counts=True)g = len(unique_x)# 计算方差VAR(S)if n == g: # 如果不存在重复点var_s = (n * (n - 1) * (2 * n + 5)) / 18else:var_s = (n * (n - 1) * (2 * n + 5) - np.sum(tp * (tp - 1) * (2 * tp + 5))) / 18# 计算z_valueif n <= 10: # n<=10属于特例z = s / (n * (n - 1) / 2)else:if s > 0:z = (s - 1) / np.sqrt(var_s)elif s < 0:z = (s + 1) / np.sqrt(var_s)else:z = 0# 计算p_value，可以选择性先对p_value进行验证#p = 2 * (1 - norm.cdf(abs(z)))# 计算Z(1-alpha/2)q=normal_ppf(1 - alpha / 2)h = abs(z) >q# trend趋势判断if (z < 0) and h:trend ="decreasing"elif (z > 0) and h:trend ="increasing"else:trend = "no trend"# slope n=len(y) t=np.arange(0,n)slope, intercept = np.polyfit(y, t, deg=1)z=round(z,2)slope=round(slope,2)return trend,z,slope优化算法主要集中在两部分
第一部分：计算S值，趋势检验方法
（二）MK趋势检验_aaakirito的博客-CSDN博客_mk检验的作者aaakirito和pymannkendall库的作者都给出了基于numpy数组的优化计算方法。其中pymannkendall库的计算方法更为优秀
以下是pymannkendall库的方法，利用了全1数组的特性进行了符号化计算和矩阵想减的方法加速了运算。
def mk_score(x):s = 0n=len(x)demo = np.ones(n) for k in range(n-1):s = s + np.sum(demo[k+1:n][x[k+1:n] > x[k]]) - np.sum(demo[k+1:n][x[k+1:n] < x[k]]) return saaakirito的思路为，由双循环改为通过复制np数组错位想减实现，大大提高了效率。作者也提出将循环缩减为log2（len(y)）的方法，但是由于增加了开销，在短序列的时候表现较差，测试时在序列达到10000个样本时依然较差。
# 计算S的值v = np.array(y[1:]).reshape((-1, 1))v = np.tile(v, len(y) - 1)for i in range(len(y) - 1):v[:, i] -= y[i]v[0:i, i] = 0v = np.sign(v)s = np.sum(v)
第二部分：优化scipy.stats 模块中标准正态分布的分位点函数norm.ppf()
经过测试该函数效率较低，耗时占总计算时间的40%多，通过在网上找了一个俄罗斯人的写的函数进行替代，效率大大提高。
但是这个函数中exit(‘R8_NORMAL_CDF_INVERSE - Fatal error!’) 语句因为使用了exit()函数导致无法序列化进而不能开启多进程，用return 代替。(jupyter中有此问题，pycharm中没有。)
def r8_huge ( ):#返回理论最大值，当p小于等于0或者大于等于1时value = 1.79769313486231571E+308return valuedef r8poly_value ( n, a, x ):value = 0.0for i in range ( n - 1, -1, -1 ):value = value * x + a[i]return valuedef normal_ppf(p):a = np.array([ \3.3871328727963666080, 1.3314166789178437745e+2, \1.9715909503065514427e+3, 1.3731693765509461125e+4, \4.5921953931549871457e+4, 6.7265770927008700853e+4, \3.3430575583588128105e本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！