一般多项式曲线的最小二乘回归(Linear Regression)

对于一般多项式：

K为多项式最高项次，a为不确定的常数项，共k+1个;

有离散数据集对应,其方差：

β为,方差函数S对β自变量第j个参数的梯度(偏导数):

当以上梯度为零时，S函数值最小，即：

中的每个每个偏导数构成一个等式：

...

则：

...

变为矩阵形式：

这样就变成线性方程求解形式，可用高斯消元等方法求得,注意在计算过程中要判断对角线上的值是否为零，如果等于零可以通过换行的方法解决;

        /// <summary>
        /// Function y = a0+a1*x+a2*x^2+ ... + an*x^n
        /// </summary>
        /// <param name="dataX">x values</param>
        /// <param name="Y">y values</param>
        /// <param name="N">Degree of a Polynomial</param>
        /// <returns>[a0,a1,...,an]</returns>
        public static double[] PolyRegress(double[] dataX, double[] Y,int N)
        {
            const double tiny = 0.00001;
            int M = dataX.Length;             // M = Number of data points
            int D = N + 1;
            double x = 0, y = 0;
            double[,] V = new double[D, D];
            double[] A = new double[D];
            for (int i = 0; i < M; i++)
            {
                x = dataX[i];
                y = Y[i];
                for(int j=0;j<D;j++)
                {
                    for(int k=0;k<D;k++)
                    {
                        V[j, k] += Math.Pow(x, j + k);
                    }
                    A[j] += y * Math.Pow(x, j);
                }

            }
            for (int i = 0; i < D; i++)
            {
                double m = V[i, i];
                if (Math.Abs(m) < tiny)
                {
                    for (int i2 = i + 1; i2 < D; i2++)
                    {
                        if (Math.Abs(V[i2, i]) > tiny)
                        {
                            double tmp = 0;
                            for (int c = 0; c < D; c++)
                            {
                                tmp = V[i, c];
                                V[i, c] = V[i2, c];
                                V[i2, c] = tmp;
                            }
                            tmp = A[i];
                            A[i] = A[i2];
                            A[i2] = tmp;
                            break;
                        }
                    }
                    m = V[i, i];
                }

                if (Math.Abs(m) > tiny)
                {
                    for (int j = i; j < D; j++)
                    {
                        V[i, j] /= m;
                    }
                    A[i] /= m;
                    for (int k = 0; k < D; k++)
                    {
                        if (k != i)
                        {
                            m = V[k, i];
                            for (int l = i; l < D; l++)
                            {
                                V[k, l] -= m * V[i, l];
                            }
                            A[k] -= m * A[i];
                        }
                    }
                }
            }
            return A;
        }

 

维基关于最小二乘的解释