KW检验
克鲁斯卡尔-沃利斯检验是一种用于比较两个或两个以上独立样本是否来自相同分布的非参数统计方法。该检验不要求样本数据满足正态分布假设,因而适用于不满足正态性或方差齐性条件的数据。K-W 检验通过对样本数据进行排序并分析其秩次来评估各组间的差异。
在这篇论文当中,通过计算每个样本的特征(例如均值),把这个特征当作一条新的样本,利用KW检验去检验其分布特征。
一、导入相关的库
import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.stats.multicomp as multi
2. 读取数据
# 适配你上传的文件路径,无需修改
file_path = "../附件一材料1.xlsx"
df = pd.read_excel(file_path)
# 自动识别磁通采样列(B0~B1023,第5列开始)
wave_cols = df.columns[4:]
3. 定义特征提取函数
这个部分用于提取特征,之后用for循环可以多次调用该函数提取特征
# 时域特征提取
def extract_time_features(signal):
mean_val = np.mean(signal) # 均值
max_val = np.max(np.abs(signal)) # 最大值(峰值)
std_val = np.std(signal) # 标准差
skew_val = stats.skew(signal) # 偏度
kurt_val = stats.kurtosis(signal) # 峰度
rms_val = np.sqrt(np.mean(signal ** 2)) # 有效值
cf_val = max_val / rms_val if rms_val != 0 else 0 # 峰值因子
return mean_val, max_val, std_val, skew_val, kurt_val, cf_val
# 频域特征提取(带宽、谐波比)
def extract_freq_features(signal, fs, n, thd_ratio):
# FFT变换
fft_complex = np.fft.fft(signal)
amp_spec = np.abs(fft_complex) / (n / 2)
amp_spec = amp_spec[:n//2]
freq_axis = np.fft.fftfreq(n, 1/fs)[:n//2]
# 基波与谐波拆分
A1 = amp_spec[1] # 基波幅值
harm_amp = amp_spec[2:] # 高次谐波幅值
# 谐波比计算
if A1 == 0:
harm_ratio = 0.0
else:
harm_total = np.sqrt(np.sum(harm_amp ** 2))
harm_ratio = harm_total / A1
# 带宽计算
threshold = A1 * thd_ratio
valid_freq = freq_axis[amp_spec > threshold]
bandwidth = valid_freq.max() - valid_freq.min() if len(valid_freq) > 0 else 0.0
return bandwidth, harm_ratio
4. 批量提取所有样本特征
从idx到df.index就是遍历所有的样本,提取所有样本的各个特征
# 存储特征结果
feature_list = []
for idx in df.index:
# 提取单条波形数据
wave_sig = df.loc[idx, wave_cols].values.astype(float)
# 提取时域特征
mean_val, max_val, std_val, skew_val, kurt_val, cf_val = extract_time_features(wave_sig)
# 提取频域特征
bandwidth, harm_ratio = extract_freq_features(wave_sig, Fs, N, band_threshold_ratio)
# 存入结果
feature_list.append([mean_val, kurt_val, max_val, cf_val, skew_val, std_val, bandwidth, harm_ratio])
# 构建特征表
feature_names = ["均值", "峰度", "最大值", "峰值因子", "偏度", "标准差", "带宽", "谐波比"]
feature_df = pd.DataFrame(feature_list, columns=feature_names)
# 拼接原始标签列
df_full = pd.concat([df[["励磁波形"]], feature_df], axis=1)
5. 执行Kruskal-Wallis全局检验
该部分的主要代码是两个for循环,第一个for循环是遍历所有的特征,第二个for循环是提取不同励磁波(正弦波,三角波,梯形波)形的数据进行检验,把它们存在一个group_data里面,后面用*group_data进行解包处理,进行KW检验
# 分组变量
group_col = "励磁波形"
groups = df_full[group_col].unique()
# 存储检验结果
kw_result = []
for feat in feature_names:
# 提取各组特征数据
group_data = []
for g in groups:
data_g = df_full[df_full[group_col] == g][feat].dropna().values #.values转换成series
group_data.append(data_g)
# KW检验
H_stat, p_kw = stats.kruskal(*group_data)
# 存储结果
kw_result.append({
"特征": feat,
"KW检验p值": round(p_kw, 6) if p_kw >= 1e-6 else p_kw,
"H统计量": round(H_stat, 4)
})
# 转为DataFrame
kw_df = pd.DataFrame(kw_result)
print("===== Kruskal-Wallis全局检验结果 =====")
print(kw_df)

浙公网安备 33010602011771号