#includeiostream#includevector#includecmath#includestdexcepttypedefstd::vectordoubleVec;typedefstd::vectorVecMat;// Mat[i][j] 第 i 行第 j 列// 1. 特殊函数lgamma / 不完全 Beta // ln GammaLanczos 近似Numerical Recipes 系数doublegammln(doublexx){staticconstdoublecof[6]{76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};doublexxx,yxx;doubletmpx5.5;tmp-(x0.5)*std::log(tmp);doubleser1.000000000190015;for(intj0;j6;j){y1.0;sercof[j]/y;}return-tmpstd::log(2.5066282746310005*ser/x);}// 正则化不完全 Beta 的连分数部分doublebetacf(doublea,doubleb,doublex){constintMAXIT300;constdoubleEPS3e-15,FPMIN1e-300;doubleqabab,qapa1.0,qama-1.0;doublec1.0;doubled1.0-qab*x/qap;if(std::fabs(d)FPMIN)dFPMIN;d1.0/d;doublehd;for(intm1;mMAXIT;m){intm22*m;doubleaam*(b-m)*x/((qamm2)*(am2));d1.0aa*d;if(std::fabs(d)FPMIN)dFPMIN;c1.0aa/c;if(std::fabs(c)FPMIN)cFPMIN;d1.0/d;h*d*c;aa-(am)*(qabm)*x/((am2)*(qapm2));d1.0aa*d;if(std::fabs(d)FPMIN)dFPMIN;c1.0aa/c;if(std::fabs(c)FPMIN)cFPMIN;d1.0/d;doubledeld*c;h*del;if(std::fabs(del-1.0)EPS)break;}returnh;}// 正则化不完全 Beta 函数 I_x(a,b)doublebetai(doublea,doubleb,doublex){if(x0.0)return0.0;if(x1.0)return1.0;doublebtstd::exp(gammln(ab)-gammln(a)-gammln(b)a*std::log(x)b*std::log(1.0-x));if(x(a1.0)/(ab2.0))returnbt*betacf(a,b,x)/a;elsereturn1.0-bt*betacf(b,a,1.0-x)/b;}// 2. 最小二乘Householder QR // 解 A*beta ≈ bA 为 m×nm≥nVeclstsq(constMatA,constVecb){intm(int)A.size();if(m0)returnVec();intn(int)A[0].size();if(n0)returnVec();Mat RA;Vec yb;for(intk0;kn;k){doublenorm0.0;for(intik;im;i)normR[i][k]*R[i][k];normstd::sqrt(norm);if(norm1e-300)continue;doublealpha(R[k][k]0.0)?-norm:norm;Vecv(m-k);for(intik;im;i)v[i-k]R[i][k];v[0]-alpha;doublevnorm20.0;for(inti0;im-k;i)vnorm2v[i]*v[i];if(vnorm21e-300)continue;for(intjk;jn;j){doubledot0.0;for(intik;im;i)dotv[i-k]*R[i][j];doublef2.0*dot/vnorm2;for(intik;im;i)R[i][j]-f*v[i-k];}doubledot0.0;for(intik;im;i)dotv[i-k]*y[i];doublef2.0*dot/vnorm2;for(intik;im;i)y[i]-f*v[i-k];}Vecbeta(n,0.0);for(intin-1;i0;--i){doublesy[i];for(intji1;jn;j)s-R[i][j]*beta[j];beta[i](std::fabs(R[i][i])1e-300)?0.0:s/R[i][i];}returnbeta;}// 3. 统计核心 doublemean(constVecv){doubles0.0;for(size_t i0;iv.size();i)sv[i];returns/v.size();}doublepearson(constVeca,constVecb){doublemamean(a),mbmean(b);doublenum0.0,da0.0,db0.0;for(size_t i0;ia.size();i){doublexaa[i]-ma,xbb[i]-mb;numxa*xb;daxa*xa;dbxb*xb;}doubledenomstd::sqrt(da*db);if(denom1e-15)return0.0;returnnum/denom;}Vecresidualize(constVecv,constMatZ){intn(int)v.size();intk(Z.empty()||Z[0].empty())?0:(int)Z[0].size();// 增广常数列截距项与 MATLAB partialcorr 的回归口径一致// 残差 v - [1 Z] * beta而非过原点回归 v ~ ZMatZa(n,Vec(k1,0.0));for(inti0;in;i){Za[i][0]1.0;for(intj0;jk;j)Za[i][j1]Z[i][j];}Vec betalstsq(Za,v);Vecr(n);for(inti0;in;i){doublesbeta[0];for(intj0;jk;j)sZ[i][j]*beta[j1];r[i]v[i]-s;}returnr;}VecgetCol(constMatM,intj){Vecc(M.size());for(size_t i0;iM.size();i)c[i]M[i][j];returnc;}// 对应 MATLAB: partialcorr(X, Y, Z)Matpartialcorr(constMatX,constMatY,constMatZ){intn(int)X.size();if((int)Y.size()!n)throwstd::runtime_error(X,Y row mismatch);intpx(int)X[0].size(),py(int)Y[0].size();std::vectorVecrx(px),ry(py);for(intj0;jpx;j)rx[j]residualize(getCol(X,j),Z);for(intj0;jpy;j)ry[j]residualize(getCol(Y,j),Z);Matrho(px,Vec(py));for(inti0;ipx;i)for(intj0;jpy;j)rho[i][j]pearson(rx[i],ry[j]);returnrho;}// t 分布双侧 p 值p I_{df/(dft²)}(df/2, 1/2)doublet_pval(doubler,intdf){if(std::fabs(r)1.0)return0.0;doublet2r*r*df/(1.0-r*r);doublexdf/(dft2);returnbetai(df/2.0,0.5,x);}Matpartialcorr_pval(constMatrho,intn,intk){intpx(int)rho.size(),py(int)rho[0].size();intdfn-k-2;if(df0)throwstd::runtime_error(Not enough degrees of freedom);Matpv(px,Vec(py));for(inti0;ipx;i)for(intj0;jpy;j)pv[i][j]t_pval(rho[i][j],df);returnpv;}// 4. 输出 voidprintMat(constMatM){for(size_t i0;iM.size();i){for(size_t j0;jM[i].size();j){if(j)std::cout ;std::coutM[i][j];}std::cout\n;}}intmain(){// n 6 的确定性数据便于复现Mat Z{{1},{2},{3},{4},{5},{6}};Mat X{{2},{4},{5},{4},{5},{7}};Mat Y{{2},{5},{6},{7},{8},{11}};Mat rhopartialcorr(X,Y,Z);Mat pvalpartialcorr_pval(rho,(int)X.size(),(int)Z[0].size());std::coutpartial correlation:\n;printMat(rho);std::cout\np-value:\n;printMat(pval);return0;}
#includeiostream#includevector#includecmath#includestdexcepttypedefstd::vectordoubleVec;typedefstd::vectorVecMat;// Mat[i][j] 第 i 行第 j 列// 1. 特殊函数lgamma / 不完全 Beta // ln GammaLanczos 近似Numerical Recipes 系数doublegammln(doublexx){staticconstdoublecof[6]{76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};doublexxx,yxx;doubletmpx5.5;tmp-(x0.5)*std::log(tmp);doubleser1.000000000190015;for(intj0;j6;j){y1.0;sercof[j]/y;}return-tmpstd::log(2.5066282746310005*ser/x);}// 正则化不完全 Beta 的连分数部分doublebetacf(doublea,doubleb,doublex){constintMAXIT300;constdoubleEPS3e-15,FPMIN1e-300;doubleqabab,qapa1.0,qama-1.0;doublec1.0;doubled1.0-qab*x/qap;if(std::fabs(d)FPMIN)dFPMIN;d1.0/d;doublehd;for(intm1;mMAXIT;m){intm22*m;doubleaam*(b-m)*x/((qamm2)*(am2));d1.0aa*d;if(std::fabs(d)FPMIN)dFPMIN;c1.0aa/c;if(std::fabs(c)FPMIN)cFPMIN;d1.0/d;h*d*c;aa-(am)*(qabm)*x/((am2)*(qapm2));d1.0aa*d;if(std::fabs(d)FPMIN)dFPMIN;c1.0aa/c;if(std::fabs(c)FPMIN)cFPMIN;d1.0/d;doubledeld*c;h*del;if(std::fabs(del-1.0)EPS)break;}returnh;}// 正则化不完全 Beta 函数 I_x(a,b)doublebetai(doublea,doubleb,doublex){if(x0.0)return0.0;if(x1.0)return1.0;doublebtstd::exp(gammln(ab)-gammln(a)-gammln(b)a*std::log(x)b*std::log(1.0-x));if(x(a1.0)/(ab2.0))returnbt*betacf(a,b,x)/a;elsereturn1.0-bt*betacf(b,a,1.0-x)/b;}// 2. 最小二乘Householder QR // 解 A*beta ≈ bA 为 m×nm≥nVeclstsq(constMatA,constVecb){intm(int)A.size();if(m0)returnVec();intn(int)A[0].size();if(n0)returnVec();Mat RA;Vec yb;for(intk0;kn;k){doublenorm0.0;for(intik;im;i)normR[i][k]*R[i][k];normstd::sqrt(norm);if(norm1e-300)continue;doublealpha(R[k][k]0.0)?-norm:norm;Vecv(m-k);for(intik;im;i)v[i-k]R[i][k];v[0]-alpha;doublevnorm20.0;for(inti0;im-k;i)vnorm2v[i]*v[i];if(vnorm21e-300)continue;for(intjk;jn;j){doubledot0.0;for(intik;im;i)dotv[i-k]*R[i][j];doublef2.0*dot/vnorm2;for(intik;im;i)R[i][j]-f*v[i-k];}doubledot0.0;for(intik;im;i)dotv[i-k]*y[i];doublef2.0*dot/vnorm2;for(intik;im;i)y[i]-f*v[i-k];}Vecbeta(n,0.0);for(intin-1;i0;--i){doublesy[i];for(intji1;jn;j)s-R[i][j]*beta[j];beta[i](std::fabs(R[i][i])1e-300)?0.0:s/R[i][i];}returnbeta;}// 3. 统计核心 doublemean(constVecv){doubles0.0;for(size_t i0;iv.size();i)sv[i];returns/v.size();}doublepearson(constVeca,constVecb){doublemamean(a),mbmean(b);doublenum0.0,da0.0,db0.0;for(size_t i0;ia.size();i){doublexaa[i]-ma,xbb[i]-mb;numxa*xb;daxa*xa;dbxb*xb;}doubledenomstd::sqrt(da*db);if(denom1e-15)return0.0;returnnum/denom;}Vecresidualize(constVecv,constMatZ){intn(int)v.size();intk(Z.empty()||Z[0].empty())?0:(int)Z[0].size();// 增广常数列截距项与 MATLAB partialcorr 的回归口径一致// 残差 v - [1 Z] * beta而非过原点回归 v ~ ZMatZa(n,Vec(k1,0.0));for(inti0;in;i){Za[i][0]1.0;for(intj0;jk;j)Za[i][j1]Z[i][j];}Vec betalstsq(Za,v);Vecr(n);for(inti0;in;i){doublesbeta[0];for(intj0;jk;j)sZ[i][j]*beta[j1];r[i]v[i]-s;}returnr;}VecgetCol(constMatM,intj){Vecc(M.size());for(size_t i0;iM.size();i)c[i]M[i][j];returnc;}// 对应 MATLAB: partialcorr(X, Y, Z)Matpartialcorr(constMatX,constMatY,constMatZ){intn(int)X.size();if((int)Y.size()!n)throwstd::runtime_error(X,Y row mismatch);intpx(int)X[0].size(),py(int)Y[0].size();std::vectorVecrx(px),ry(py);for(intj0;jpx;j)rx[j]residualize(getCol(X,j),Z);for(intj0;jpy;j)ry[j]residualize(getCol(Y,j),Z);Matrho(px,Vec(py));for(inti0;ipx;i)for(intj0;jpy;j)rho[i][j]pearson(rx[i],ry[j]);returnrho;}// t 分布双侧 p 值p I_{df/(dft²)}(df/2, 1/2)doublet_pval(doubler,intdf){if(std::fabs(r)1.0)return0.0;doublet2r*r*df/(1.0-r*r);doublexdf/(dft2);returnbetai(df/2.0,0.5,x);}Matpartialcorr_pval(constMatrho,intn,intk){intpx(int)rho.size(),py(int)rho[0].size();intdfn-k-2;if(df0)throwstd::runtime_error(Not enough degrees of freedom);Matpv(px,Vec(py));for(inti0;ipx;i)for(intj0;jpy;j)pv[i][j]t_pval(rho[i][j],df);returnpv;}// 4. 输出 voidprintMat(constMatM){for(size_t i0;iM.size();i){for(size_t j0;jM[i].size();j){if(j)std::cout ;std::coutM[i][j];}std::cout\n;}}intmain(){// n 6 的确定性数据便于复现Mat Z{{1},{2},{3},{4},{5},{6}};Mat X{{2},{4},{5},{4},{5},{7}};Mat Y{{2},{5},{6},{7},{8},{11}};Mat rhopartialcorr(X,Y,Z);Mat pvalpartialcorr_pval(rho,(int)X.size(),(int)Z[0].size());std::coutpartial correlation:\n;printMat(rho);std::cout\np-value:\n;printMat(pval);return0;}
