企业官网建设流程全解析

// fitcgam.cpp — 完全自包含，无需任何第三方库// 编译: g++ -std=c++17 -O2 -o fitcgam fitcgam.cpp#include<vector>#include<cmath>#include<numeric>#include<algorithm>#include<stdexcept>#include<iostream>#include<random>#include<iomanip>// ============================================================// 特殊函数实现 (替代 Boost.Math)// ============================================================// digamma 函数: ψ(x) = d/dx ln Γ(x)// 递推 ψ(x) = ψ(x+1) - 1/x 将 x 提升到 >= 8，再用渐近展开staticdoubledigamma(doublex){doubleresult=0.0;while(x<8.0){result-=1.0/x;x+=1.0;}constdoubleinv=1.0/x;constdoubleu=inv*inv;result+=std::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) = ψ'(x+1) + 1/x²，再用渐近展开staticdoubletrigamma(doublex){doubleresult=0.0;while(x<8.0){result+=1.0/(x*x);x+=1.0;}constdoubleinv=1.0/x;constdoubleu=inv*inv;result+=inv+0.5*u+inv*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:staticconstexprintkMaxIter=100;staticconstexprdoublekTol=1e-8;staticGammaFitResultfit(conststd::vector<double>&data,doubleshape0=-1.0){constintn=static_cast<int>(data.size());if(n<2)throwstd::invalid_argument("需要至少 2 个数据点");for(doublex:data){if(x<=0.0)throwstd::invalid_argument("伽马分布要求所有数据 > 0");}constdoublexbar=std::accumulate(data.begin(),data.end(),0.0)/n;doublelogSum=0.0;for(doublex:data)logSum+=std::log(x);constdoublemeanLog=logSum/n;// 目标方程: ln(α) - ψ(α) = ln(x̄) - mean(ln x)constdoublerhs=std::log(xbar)-meanLog;// 初始值: 矩估计 或 用户指定doublealpha;if(shape0>0.0){alpha=shape0;}else{doublesumSq=0.0;for(doublex:data)sumSq+=(x-xbar)*(x-xbar);constdoublevar=sumSq/(n-1);alpha=(var>1e-12)?(xbar*xbar/var):1.0;if(alpha<=0.0)alpha=1.0;}// Newton-Raphson 迭代求 αboolconverged=false;intiter=0;for(iter=0;iter<kMaxIter;++iter){constdoublefval=std::log(alpha)-digamma(alpha)-rhs;constdoublederiv=1.0/alpha-trigamma(alpha);if(std::abs(deriv)<1e-15){alpha*=1.01;continue;}constdoublestep=fval/deriv;doublealphaNew=alpha-step;if(alphaNew<=1e-10)alphaNew=alpha*0.5;// 阻尼if(std::abs(alphaNew-alpha)<kTol*std::max(1.0,alpha)){alpha=alphaNew;converged=true;++iter;break;}alpha=alphaNew;}constdoublebeta=xbar/alpha;constdoublelogLik=computeLogLik(data,alpha,beta);return{alpha,beta,logLik,converged,iter};}staticdoublecomputeLogLik(conststd::vector<double>&data,doublealpha,doublebeta){constintn=static_cast<int>(data.size());doublesumLogX=0.0,sumX=0.0;for(doublex:data){sumLogX+=std::log(x);sumX+=x;}return(alpha-1.0)*sumLogX-sumX/beta-n*alpha*std::log(beta)-n*std::lgamma(alpha);}staticdoublecomputeAIC(doublelogLik,intnumParams=2){return-2.0*logLik+2.0*numParams;}staticdoublecomputeBIC(doublelogLik,intn,intnumParams=2){return-2.0*logLik+numParams*std::log(static_cast<double>(n));}};// ============================================================// main// ============================================================intmain(){// 生成模拟数据: Gamma(shape=5.0, scale=2.0)constdoubletrueShape=5.0;constdoubletrueScale=2.0;constintN=1000;std::mt19937rng(42);std::gamma_distribution<double>dist(trueShape,trueScale);std::vector<double>data(N);for(auto&x:data)x=dist(rng);autoresult=GammaFitter::fit(data);std::cout<<std::fixed<<std::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::cout<<"AIC = "<<GammaFitter::computeAIC(result.logLik)<<"\n";std::cout<<"BIC = "<<GammaFitter::computeBIC(result.logLik,N)<<"\n";return0;}