// fitcgam.cpp — 完全自包含无需任何第三方库// 编译: g -stdc17 -O2 -o fitcgam fitcgam.cpp#includevector#includecmath#includenumeric#includealgorithm#includestdexcept#includeiostream#includerandom#includeiomanip// // 特殊函数实现 (替代 Boost.Math)// // digamma 函数: ψ(x) d/dx ln Γ(x)// 递推 ψ(x) ψ(x1) - 1/x 将 x 提升到 8再用渐近展开staticdoubledigamma(doublex){doubleresult0.0;while(x8.0){result-1.0/x;x1.0;}constdoubleinv1.0/x;constdoubleuinv*inv;resultstd::log(x)-0.5*inv-u*(1.0/12.0-u*(1.0/120.0-u*(1.0/252.0-u*(1.0/240.0-u*(1.0/132.0)))));returnresult;}// trigamma 函数: ψ(x) d²/dx² ln Γ(x)// 递推 ψ(x) ψ(x1) 1/x²再用渐近展开staticdoubletrigamma(doublex){doubleresult0.0;while(x8.0){result1.0/(x*x);x1.0;}constdoubleinv1.0/x;constdoubleuinv*inv;resultinv0.5*uinv*u*(1.0/6.0-u*(1.0/30.0-u*(1.0/42.0-u*(1.0/30.0-u*(5.0/66.0)))));returnresult;}// // GammaFitter: 伽马分布 MLE 拟合 (等效于 MATLAB gamfit)// structGammaFitResult{doubleshape;// alpha (形状参数)doublescale;// beta (尺度参数)doublelogLik;// 对数似然值boolconverged;intiterations;};classGammaFitter{public:staticconstexprintkMaxIter100;staticconstexprdoublekTol1e-8;staticGammaFitResultfit(conststd::vectordoubledata,doubleshape0-1.0){constintnstatic_castint(data.size());if(n2)throwstd::invalid_argument(需要至少 2 个数据点);for(doublex:data){if(x0.0)throwstd::invalid_argument(伽马分布要求所有数据 0);}constdoublexbarstd::accumulate(data.begin(),data.end(),0.0)/n;doublelogSum0.0;for(doublex:data)logSumstd::log(x);constdoublemeanLoglogSum/n;// 目标方程: ln(α) - ψ(α) ln(x̄) - mean(ln x)constdoublerhsstd::log(xbar)-meanLog;// 初始值: 矩估计 或 用户指定doublealpha;if(shape00.0){alphashape0;}else{doublesumSq0.0;for(doublex:data)sumSq(x-xbar)*(x-xbar);constdoublevarsumSq/(n-1);alpha(var1e-12)?(xbar*xbar/var):1.0;if(alpha0.0)alpha1.0;}// Newton-Raphson 迭代求 αboolconvergedfalse;intiter0;for(iter0;iterkMaxIter;iter){constdoublefvalstd::log(alpha)-digamma(alpha)-rhs;constdoublederiv1.0/alpha-trigamma(alpha);if(std::abs(deriv)1e-15){alpha*1.01;continue;}constdoublestepfval/deriv;doublealphaNewalpha-step;if(alphaNew1e-10)alphaNewalpha*0.5;// 阻尼if(std::abs(alphaNew-alpha)kTol*std::max(1.0,alpha)){alphaalphaNew;convergedtrue;iter;break;}alphaalphaNew;}constdoublebetaxbar/alpha;constdoublelogLikcomputeLogLik(data,alpha,beta);return{alpha,beta,logLik,converged,iter};}staticdoublecomputeLogLik(conststd::vectordoubledata,doublealpha,doublebeta){constintnstatic_castint(data.size());doublesumLogX0.0,sumX0.0;for(doublex:data){sumLogXstd::log(x);sumXx;}return(alpha-1.0)*sumLogX-sumX/beta-n*alpha*std::log(beta)-n*std::lgamma(alpha);}staticdoublecomputeAIC(doublelogLik,intnumParams2){return-2.0*logLik2.0*numParams;}staticdoublecomputeBIC(doublelogLik,intn,intnumParams2){return-2.0*logLiknumParams*std::log(static_castdouble(n));}};// // main// intmain(){// 生成模拟数据: Gamma(shape5.0, scale2.0)constdoubletrueShape5.0;constdoubletrueScale2.0;constintN1000;std::mt19937rng(42);std::gamma_distributiondoubledist(trueShape,trueScale);std::vectordoubledata(N);for(autox:data)xdist(rng);autoresultGammaFitter::fit(data);std::coutstd::fixedstd::setprecision(6);std::cout 伽马分布拟合结果 \n;std::cout形状参数 alpha result.shape (真实值: trueShape)\n;std::cout尺度参数 beta result.scale (真实值: trueScale)\n;std::cout对数似然 result.logLik\n;std::cout收敛 (result.converged?是:否) (迭代 result.iterations 次)\n;std::coutAIC GammaFitter::computeAIC(result.logLik)\n;std::coutBIC GammaFitter::computeBIC(result.logLik,N)\n;return0;}