// removeoutliers.cpp// 对应 MATLAB rmoutliers(A) 的 C++ 实现 + 演示程序// 编译: g++ -std=c++17 -O2 -o removeoutliers_demo#include<iostream>#include<iomanip>#include<vector>#include<algorithm>#include<cmath>#include<numeric>#include<stdexcept>#include<string>/** * 对应 MATLAB rmoutliers(A) 的 C++ 实现 * 默认方法：MAD（中位数绝对偏差），阈值因子 3 * 支持方法：median(MAD)、mean(均值±标准差)、quartiles(四分位数) */classRemoveOutliers{public:// 检测方法enumclassMethod{Median,// 默认，MADMean,// 均值 ± 3σQuartiles// 四分位数法};structOptions{Method method=Method::Median;doublethresholdFactor=3.0;// median/mean 默认 3；quartiles 默认 1.5booluseScaledMAD=true;// 是否使用换算 MAD（乘常数 c）intdim=1;// 1 按行移除，2 按列移除（矩阵）};structResult{std::vector<std::vector<double>>cleanedData;// 清洗后数据std::vector<bool>outlierMask;// 离群值逻辑掩码（展平）std::vector<int>removedIndices;// 被移除的行/列索引};public:/** * 对向量或矩阵进行离群值检测与移除 * @param data 输入数据（二维：行×列） * @param opts 选项 * @return 清洗后的数据及元信息 */staticResultprocess(conststd::vector<std::vector<double>>&data){returnprocess(data,Options{});}staticResultprocess(conststd::vector<std::vector<double>>&data,constOptions&opts){if(data.empty()||data[0].empty())throwstd::invalid_argument("Input data is empty.");constsize_t rows=data.size();constsize_t cols=data[0].size();// 校验：所有行长度一致for(constauto&row:data){if(row.size()!=cols)throwstd::invalid_argument("Rows have inconsistent lengths.");}Result result;if(opts.dim==1){// 按列检测，移除整行（MATLAB 默认行为）result=processByColumn(data,opts);}elseif(opts.dim==2){// 按行检测，移除整列result=processByRow(data,opts);}else{throwstd::invalid_argument("dim must be 1 or 2.");}returnresult;}// ========== 核心统计函数 ==========/** * 计算中位数（自动排序副本，不修改原数据） */staticdoublemedian(std::vector<double>v){if(v.empty())return0.0;size_t n=v.size();std::sort(v.begin(),v.end());if(n%2==1){returnv[n/2];}else{return0.5*(v[n/2-1]+v[n/2]);}}/** * 计算换算 MAD（Median Absolute Deviation） * scaledMAD = c * median(|x - median(x)|) * c = 1 / (sqrt(2) * erfcinv(3/2)) ≈ 1.482602218505602 */staticdoublescaledMAD(conststd::vector<double>&data){if(data.size()<2)return0.0;doublemed=median(data);std::vector<double>absDev;absDev.reserve(data.size());for(doublex:data){absDev.push_back(std::fabs(x-med));}doublerawMAD=median(absDev);// 换算常数 cconstdoublec=1.482602218505602;returnc*rawMAD;}/** * 计算算术均值 */staticdoublemean(conststd::vector<double>&data){if(data.empty())return0.0;doublesum=std::accumulate(data.begin(),data.end(),0.0);returnsum/static_cast<double>(data.size());}/** * 计算样本标准差（无偏估计，除以 n-1，与 MATLAB std 默认一致） */staticdoublestddev(conststd::vector<double>&data){size_t n=data.size();if(n<2)return0.0;doublem=mean(data);doublesqSum=0.0;for(doublex:data){doubled=x-m;sqSum+=d*d;}returnstd::sqrt(sqSum/static_cast<double>(n-1));}/** * 计算指定百分位数（线性插值法，与 MATLAB prctile 一致） */staticdoublepercentile(std::vector<double>v,doublep){if(v.empty())return0.0;if(p<=0.0)return*std::min_element(v.begin(),v.end());if(p>=100.0)return*std::max_element(v.begin(),v.end());std::sort(v.begin(),v.end());doublerank=(p/100.0)*(static_cast<double>(v.size())-1.0);size_t lo=static_cast<size_t>(std::floor(rank));size_t hi=static_cast<size_t>(std::ceil(rank));doublefrac=rank-static_cast<double>(lo);if(lo==hi)returnv[lo];returnv[lo]*(1.0-frac)+v[hi]*frac;}private:// ========== 按列处理（移除行） ==========staticResultprocessByColumn(conststd::vector<std::vector<double>>&data,constOptions&opts){constsize_t rows=data.size();constsize_t cols=data[0].size();std::vector<bool>rowOutlier(rows,false);// 对每一列独立检测for(size_t j=0;j<cols;++j){std::vector<double>colData;colData.reserve(rows);for(size_t i=0;i<rows;++i)colData.push_back(data[i][j]);std::vector<bool>colMask=detectOutliers(colData,opts);for(size_t i=0;i<rows;++i){if(colMask[i])rowOutlier[i]=true;// 任一列检出则移除整行}}// 构建结果Result res;res.outlierMask=rowOutlier;for(size_t i=0;i<rows;++i){if(rowOutlier[i]){res.removedIndices.push_back(static_cast<int>(i));}else{res.cleanedData.push_back(data[i]);}}returnres;}// ========== 按行处理（移除列） ==========staticResultprocessByRow(conststd::vector<std::vector<double>>&data,constOptions&opts){constsize_t rows=data.size();constsize_t cols=data[0].size();std::vector<bool>colOutlier(cols,false);for(size_t i=0;i<rows;++i){std::vector<bool>rowMask=detectOutliers(data[i],opts);for(size_t j=0;j<cols;++j){if(rowMask[j])colOutlier[j]=true;}}Result res;for(size_t i=0;i<rows;++i){std::vector<double>newRow;for(size_t j=0;j<cols;++j){if(!colOutlier[j])newRow.push_back(data[i][j]);}res.cleanedData.push_back(std::move(newRow));}for(size_t j=0;j<cols;++j){if(colOutlier[j])res.removedIndices.push_back(static_cast<int>(j));}res.outlierMask=colOutlier;returnres;}// ========== 离群值检测核心 ==========staticstd::vector<bool>detectOutliers(conststd::vector<double>&data,constOptions&opts){size_t n=data.size();std::vector<bool>mask(n,false);if(n<2)returnmask;switch(opts.method){caseMethod::Median:{doublemed=median(data);doublemadVal=opts.useScaledMAD?scaledMAD(data):[&](){std::vector<double>absDev;for(doublex:data)absDev.push_back(std::fabs(x-med));returnmedian(absDev);}();if(madVal<=0.0)returnmask;// 所有值相同，无离群值doublethr=opts.thresholdFactor*madVal;for(size_t i=0;i<n;++i){if(std::fabs(data[i]-med)>thr)mask[i]=true;}break;}caseMethod::Mean:{doublem=mean(data);doubles=stddev(data);if(s<=0.0)returnmask;doublethr=opts.thresholdFactor*s;for(size_t i=0;i<n;++i){if(std::fabs(data[i]-m)>thr)mask[i]=true;}break;}caseMethod::Quartiles:{doubleq1=percentile(data,25.0);doubleq3=percentile(data,75.0);doubleiqr=q3-q1;doublelower=q1-opts.thresholdFactor*iqr;doubleupper=q3+opts.thresholdFactor*iqr;for(size_t i=0;i<n;++i){if(data[i]<lower||data[i]>upper)mask[i]=true;}break;}}returnmask;}};// ============================ 演示主程序 ============================intmain(){// 示例数据（含离群值）std::vector<std::vector<double>>data={{1.0,1.0},// 0{34.0,0.0},// 1{35.0,35.0},// 2{35.0,35.0},// 3{33.0,33.0},// 4{34.0,34.0},// 5{160.0,160.0},// 6 离群值{38.0,38.0},// 7{35.0,35.0},// 8{95.0,35.0},// 9 第1列离群{36.0,36.0},// 10{150.0,150.0},// 11 离群值{39.0,39.0},// 12{41.0,41.0},// 13{42.0,40.0},// 14{0.0,0.0},// 15 离群值{200.0,200.0}// 16 离群值};// 默认 MAD 方法autoresult=RemoveOutliers::process(data);std::cout<<"=== 默认 MAD 方法（thresholdFactor=3）===\n";std::cout<<"被移除的行索引：";for(intidx:result.removedIndices)std::cout<<idx<<" ";std::cout<<"\n\n清洗后数据：\n";for(constauto&row:result.cleanedData){for(doublev:row)std::cout<<std::setw(6)<<v<<" ";std::cout<<"\n";}// 均值 ± 3σ 方法RemoveOutliers::Options meanOpt;meanOpt.method=RemoveOutliers::Method::Mean;autoresultMean=RemoveOutliers::process(data,meanOpt);std::cout<<"\n=== Mean 方法 ===\n移除行数："<<resultMean.removedIndices.size()<<"\n";// 四分位数方法RemoveOutliers::Options quartOpt;quartOpt.method=RemoveOutliers::Method::Quartiles;quartOpt.thresholdFactor=1.5;autoresultQuart=RemoveOutliers::process(data,quartOpt);std::cout<<"\n=== Quartiles 方法（1.5×IQR）===\n移除行数："<<resultQuart.removedIndices.size()<<"\n";return0;}