#includevector#includecmath#includestdexcept#includelimits#includeiostream// 计算单个序列的峰度忽略 NaN// flag 1: 有偏估计MATLAB 默认// flag 0: 偏差校正估计doublekurtosis_vector(conststd::vectordoublex,intflag1){std::vectordoublevalid;valid.reserve(x.size());for(doublev:x){if(!std::isnan(v))valid.push_back(v);}constsize_t nvalid.size();if(n4flag0){returnstd::numeric_limitsdouble::quiet_NaN();}if(n0)returnstd::numeric_limitsdouble::quiet_NaN();doublemean0.0;for(doublev:valid)meanv;mean/n;doublem20.0,m40.0;for(doublev:valid){doubledv-mean;doubled2d*d;m2d2;m4d2*d2;}m2/n;m4/n;if(m20.0)returnstd::numeric_limitsdouble::quiet_NaN();doublek1m4/(m2*m2);if(flag1){returnk1;}doublek0(static_castdouble(n-1)/((n-2)*(n-3)))*((n1)*k1-3.0*(n-1))3.0;returnk0;}std::vectordoublekurtosis_matrix(conststd::vectordoubleX,size_t nrows,size_t ncols,intflag1,intdim-1){if(X.size()!nrows*ncols){throwstd::invalid_argument(X size does not match nrows * ncols);}if(dim-1){if(nrows1)dim1;elseif(ncols1)dim2;elsedim1;}if(dim!1dim!2){throwstd::invalid_argument(dim must be 1 or 2);}if(dim1){std::vectordoubleresult(ncols);std::vectordoublecolumn(nrows);for(size_t j0;jncols;j){for(size_t i0;inrows;i){column[i]X[i*ncolsj];}result[j]kurtosis_vector(column,flag);}returnresult;}else{std::vectordoubleresult(nrows);std::vectordoublerow(ncols);for(size_t i0;inrows;i){for(size_t j0;jncols;j){row[j]X[i*ncolsj];}result[i]kurtosis_vector(row,flag);}returnresult;}}voidprint_vec(conststd::vectordoublev){std::cout[;for(size_t i0;iv.size();i){if(i0)std::cout, ;std::coutv[i];}std::cout];}intmain(){std::vectordoublex{1,2,3,4,5};std::vectordoublexn{1,2,std::numeric_limitsdouble::quiet_NaN(),3,4,5};std::vectordoublexc{2,2,2,2};std::vectordoubleX{1,2,3,4,5,6};// 2x3行主序std::coutkurtosis_vector(x, 1) kurtosis_vector(x,1)\n;std::coutkurtosis_vector(x, 0) kurtosis_vector(x,0)\n;std::cout\n;std::coutkurtosis_vector(xn, 1) kurtosis_vector(xn,1)\n;std::coutkurtosis_vector(xn, 0) kurtosis_vector(xn,0)\n;std::cout\n;std::coutkurtosis_vector(xc, 1) kurtosis_vector(xc,1)\n;std::cout\n;std::coutkurtosis_matrix(X, 2, 3, 1, -1) ;print_vec(kurtosis_matrix(X,2,3,1,-1));std::cout\n;std::coutkurtosis_matrix(X, 2, 3, 0, -1) ;print_vec(kurtosis_matrix(X,2,3,0,-1));std::cout\n;std::cout\n;std::coutkurtosis_matrix(X, 2, 3, 1, 2) ;print_vec(kurtosis_matrix(X,2,3,1,2));std::cout\n;std::coutkurtosis_matrix(X, 2, 3, 0, 2) ;print_vec(kurtosis_matrix(X,2,3,0,2));std::cout\n;return0;}