[部分抄袭] 线性拟合，过原点的线性拟合，多项式曲线拟合，标准差，变异系数

曲线拟合 部分只是简单改了一下，没有验证算法，需要小心

#include
#include
#include
#include

// 求和
double StaSum( const double vals[], size_t n )
{
    double sum = 0;

    double C = 0.0;
    for( size_t i=0; i!=n; ++i )
    {
        double Y = vals[i] - C;
        double T = sum + Y;
        C = T - sum - Y;
        sum = T;
    }

    return sum;
}
// 求平均值
double StaAverage( const double vals[], size_t n )
{
 assert( n!=0 );
    if(n==0) return 0.0;
    return StaSum(vals,n)/n;
}
// 求平方和
double StaSumSq( const double vals[], size_t n )
{
    double sum = 0;

    double C = 0.0;
    for( size_t i=0; i!=n; ++i )
    {
        double Y = vals[i]*vals[i] - C;
        double T = sum + Y;
        C = T - sum - Y;
        sum = T;
    }

    return sum;
}
// 求a*b的和
double StaSumAB( const double a[], const double b[], size_t n )
{
    double sum = 0;

    double C = 0.0;
    for( size_t i=0; i!=n; ++i )
    {
        double Y = a[i]*b[i] - C;
        double T = sum + Y;
        C = T - sum - Y;
        sum = T;
    }

    return sum;
}
// 求总体标准差
double StaStdevp( const double vals[], size_t n )
{
 assert( n!=0 );
    if(n==0) return 0.0;
    double avg = StaAverage(vals,n);
    return sqrt( StaSumSq(vals,n)/n - avg*avg );
}
// 求样本标准差
double StaStdev( const double vals[], size_t n )
{
 assert( n>1 );
    if(n<=1) return 0.0;
    double sum = StaSum(vals,n);
    return sqrt( (StaSumSq(vals,n)-sum*sum/n) / (n-1) );
}
// 求变异系数
double StaCvrsd( double vals[], size_t n )
{
    double sum = StaSum(vals,n);
 assert( sum!=0 );
    if(sum==0) return 0.0;
    return sqrt( StaSumSq(vals,n)/(sum*sum)*n - 1 );
}

// 功能：从数组(vals,n)中选出d（不小于2）个元素求变异系数，使得变异系数最低
// input:
//     数组vals，数组元素数目n；选出元素数目d
// output:
//     数组下标起始索引idx，长度，最低变异系数cv
void StaExtLdla( double vals[], size_t n, size_t d, size_t& idx, size_t& num, double& cv )
{
 assert( n>1 );
 if( n <= 1 )
 {
  idx = 0;
  num = n;
  cv  = 0;
  return;
 }
 assert( d>=2 );
 num = std::max(d,2u);

    // 排序
    std::sort( vals+0, vals+n );

    // 只取num个元素，最小cv值为多少
    cv = std::numeric_limits::max();
    for( size_t i=0; i<=n-num; ++i )
    {
        double cv_tmp = StaCvrsd( vals+i, num );
        if( cv_tmp < cv )
        {
            cv = cv_tmp;
            idx = i;
        }
    }
}
// 功能：从数组(vals,n)中选出尽可能多（不少于d个）的元素，使得变异系数低于指定值(mcv)。若无法低于mcv，则不剔除
// input:
//     数组vals，数组元素数目n，至少元素数目d，期待变异系数不大于mcv
// output:
//     数组下标起始索引idx，数组元素数目，平均值，变异系数cv
void StaExtLdlb( double vals[], size_t n, size_t d, double mcv, size_t& idx, size_t& num, double& avg, double& cv )
{
 assert( n>1 );
 if( n <= 1 )
 {
  idx = 0;
  num = n;
  avg = (n==0?0:vals[0]);
  cv  = 0;
  return;
 }
 if( d > n )
 {
  idx = 0;
  num = n;
  avg = StaAverage(vals+idx,num);
  cv = StaCvrsd( vals+idx, num );
  return;
 }

    // 排序
    std::sort( vals+0, vals+n );

    idx = 0;
    num = n;
    // 分别去掉  0, 1, 2, ……, n-d 个元素
    for( size_t e=0; e<=n-d; ++e )
    {
        // 去除e个元素后，最小cv值为多少
        double cv_min = std::numeric_limits::max();
        size_t cv_idx = 0;
        for( size_t i=0; i<=e; ++i )
        {
            double cv_tmp = StaCvrsd( vals+i, n-e );
            if( cv_tmp < cv_min )
            {
                cv_min = cv_tmp;
                cv_idx = i;
            }
        }
        // 如果有符合要求的
        if( cv_min <= mcv )
        {
            idx = cv_idx;
            num = n-e;
            break;
        }
    }
    avg = StaAverage(vals+idx,num);
    cv = StaCvrsd( vals+idx, num );
}

// 线性拟合
void StaFittingLinear( double x[], double y[], size_t n, double& b, double& k )
{
    double sx = StaSum(x,n);
    double sy = StaSum(y,n);
    double sxx = StaSumSq(x,n);
    double sxy = StaSumAB(x,y,n);

    b = (sx*sxy - sy*sxx)/(sx*sx-n*sxx);
    k = (sx*sy - n*sxy)/(sx*sx-n*sxx);
    //相关系数 = (sxy-sx*sy/n) / sqrt((sxx-sx*sx/n)*(syy-sy*sy/n));
}

// 过原点的线性拟合
void StaFittingLinearOrigin( double x[], double y[], size_t n, double& k )
{
    double sxx = StaSumSq(x,n);
    double sxy = StaSumAB(x,y,n);
    k = sxy / sxx;
}

// 曲线拟合
// out: 输出从0次幂到order-1次幂的系数
void StaFittingCurve( const double x[], const double y[], size_t n, size_t order, double out[] )
{
    double* b = new double[order*n];
    double* A = new double[order*order];
    double* B = new double[order];
    double* Y = new double[order];
    double* X = new double[order];

    for( size_t i=0; i    {
        for( size_t j=0; j        {
            b[i*n+j] = pow(x[j],(int)i);
        }
    }

    for( size_t i=0; i    {
        for( size_t j=0; j        {
            double s = 0.0;
            for( size_t k=0; k            {
                s += b[i*n+k] * b[j*n+k];
            }
            A[i*order+j] = s;
        }
    }

    for( size_t i=0; i    {
        double s = 0.0;
        for( size_t k=0; k        {
            s += y[k]*b[i*n+k];
        }
        B[i] = s;
    }

    for( size_t i=1; i    {
        A[i*order+0] /= A[0*order+0];
    }

    for( size_t i=1; i    {
        for( size_t j=i; j        {
            double s = 0.0;
            for( size_t k=0; k                s += A[k*order+j]*A[i*order+k];
            A[i*order+j] -= s;
            if( j+1 != order )
            {
                s = 0.0;
                for( size_t k=0; k                    s += A[k*order+i]*A[(j+1)*order+k];
                A[(j+1)*order+i] = (A[(j+1)*order+i]-s)/A[i*order+i];
            }
        }
    }

    Y[0] = B[0];
    for( size_t i=1; i    {
        double s = 0.0;
        for( size_t j=0; j            s += A[i*order+j]*Y[j];
        Y[i] = B[i]-s;
    }

    X[order-1] = Y[order-1]/A[(order-1)*order+(order-1)];
    for( int i=(int)order-2; i>=0; --i )
    {
        double s = 0.0;
        for( size_t j=i+1; j            s += A[i*order+j]*X[j];
        X[i] = (Y[i]-s)/A[i*order+i];
    }

    for( size_t i=0; i        out[i] = X[i];

    delete[] b;
    delete[] A;
    delete[] B;
    delete[] Y;
    delete[] X;
}

// 计算拟合线的相关系数的平方
double StaFittingCurveCorrel2( const double x[], const double y[], size_t n, size_t order, const double pcs[] )
{
    double sa = 0.0; // y和曲线的差的平方和
    {
        double C = 0.0;
        for( size_t i=0; i        {
            double yi = 0.0;
            {
                double C = 0.0;
                double c = 1.0;
                for( size_t j=0; j                {
                    double Y = pcs[j]*c - C;
                    double T = yi + Y;
                    C = T - yi - Y;
                    yi = T;
                }
            }

            double Y = (y[i]-yi)*(y[i]-yi) - C;
            double T = sa + Y;
            C = T - sa - Y;
            sa = T;
        }
    }
    double sb = StaAverage(y,n); // y的平均值
    double sc = 0.0; // y和平均值的差的平方和
    {
        double C = 0.0;
        for( size_t i=0; i        {
            double Y = (y[i]-sb)*(y[i]-sb) - C;
            double T = sc + Y;
            C = T - sc - Y;
            sc = T;
        }
    }

    return 1-sa/sc;
}

// 计算拟合线的相关系数
double StaFittingCurveCorrel( const double x[], const double y[], size_t n, int order, const double pcs[] )
{
    return sqrt( StaFittingCurveCorrel2(x,y,n,order,pcs) );
}

//////////////////// 以下为测试 ////////////////////
//*
#include
#include
using namespace std;

int main(void)
{
    double buf[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
    size_t n = _countof(buf);

    cout << "和         = " << StaSum(buf,n)      << endl; // 55
    cout << "均值       = " << StaAverage(buf,n)  << endl; // 5.5
    cout << "平方和     = " << StaSumSq(buf,n)    << endl; // 385
    cout << "总体标准差 = " << StaStdevp(buf,n)   << endl; // 2.872281
    cout << "样本标准差 = " << StaStdev(buf,n)    << endl; // 3.02765
    cout << "变异系数   = " << StaCvrsd(buf,n)    << endl; // 0.522233
    cout << endl;
    {
        // 依次去掉部分离散值，看看cv变化情况
        double buf[] = { 89,90, 99,100,101, 110};
        size_t idx, num;
        double cv;
        for( size_t d=_countof(buf); d>=4; --d )
        {
            StaExtLdla( buf, _countof(buf), d, idx, num, cv  );

            for( size_t i=idx; i                cout << buf[i] << ',';
            cout << "------- cv=" << cv << endl;
        }
    }
    cout << endl;
    {
        // 如果要达到低于指定cv的要求，应当去掉哪些离散值
        double buf[] = { 89,90, 99,100,101, 110};
        size_t idx, num;
        double avg, cv;
        StaExtLdlb( buf, _countof(buf), 4, 0.05, idx, num, avg, cv );
        for( size_t i=idx; i            cout << buf[i] << ',';
        cout << "------- cv=" << cv << endl;
    }
    cout << endl;
    {
        double x[] = { -2,   -1,  0,   1,   2   };
        double y[] = { -0.1, 0.1, 0.4, 0.9, 1.6 };
        assert( _countof(x) == _countof(y) );

        double b, k;
        StaFittingLinear( x, y, _countof(x), b, k );
        cout << b << ", " << k << ", " << endl;

        double out[2];
        StaFittingCurve( x, y, _countof(x), _countof(out), out );
        for( size_t i=0; i<_countof(out); ++i )
            cout << out[i] << ", ";
        cout << endl;
        cout << "CorrelationCoefficient = " << StaFittingCurveCorrel( x, y, _countof(x), _countof(out), out ) << endl;
    }
 cout << endl;
    {
        double x[] = { 1, 2 };
        double y[] = { 1, 1 };
        assert( _countof(x) == _countof(y) );

        double b, k;
        StaFittingLinear( x, y, _countof(x), b, k );
        cout << b << ", " << k << ", " << endl;

        double out[2];
        StaFittingCurve( x, y, _countof(x), _countof(out), out );
        for( size_t i=0; i<_countof(out); ++i )
            cout << out[i] << ", ";
        cout << endl;
        cout << "CorrelationCoefficient = " << StaFittingCurveCorrel( x, y, _countof(x), _countof(out), out ) << endl;
    }

    return 0;
}
//*/