最小二乘法拟合二元多次曲线

来源：互联网发布：欧弟知乎编辑：程序博客网时间：2024/04/28 04:56
原文出处不详，数学原理大概不懂，代码有效。作用为已知一堆样本数据，拟合一个近似的2元n次函数表达式。
 public class Line    {        ///<summary>        ///用最小二乘法拟合二元多次曲线        ///</summary>        ///<param name="arrX">已知点的x坐标集合</param>        ///<param name="arrY">已知点的y坐标集合</param>        ///<param name="length">已知点的个数</param>        ///<param name="dimension">方程的最高次数</param>        public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线        {            int n = dimension + 1;                  //dimension次方程需要求 dimension+1个 系数            double[,] Guass = new double[n, n + 1];      //高斯矩阵 例如：y=a0+a1*x+a2*x*x            for (int i = 0; i < n; i++)            {                int j;                for (j = 0; j < n; j++)                {                    Guass[i, j] = SumArr(arrX, j + i, length);                }                Guass[i, j] = SumArr(arrX, i, arrY, 1, length);            }            return ComputGauss(Guass, n);        }        public static double SumArr(double[] arr, int n, int length) //求数组的元素的n次方的和        {            double s = 0;            for (int i = 0; i < length; i++)            {                if (arr[i] != 0 || n != 0)                    s = s + Math.Pow(arr[i], n);                else                    s = s + 1;            }            return s;        }        public static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)        {            double s = 0;            for (int i = 0; i < length; i++)            {                if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))                    s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);                else                    s = s + 1;            }            return s;        }        public static double[] ComputGauss(double[,] Guass, int n)        {            int i, j;            int k, m;            double temp;            double max;            double s;            double[] x = new double[n];            for (i = 0; i < n; i++) x[i] = 0.0;//初始化            for (j = 0; j < n; j++)            {                max = 0;                k = j;                for (i = j; i < n; i++)                {                    if (Math.Abs(Guass[i, j]) > max)                    {                        max = Guass[i, j];                        k = i;                    }                }                if (k != j)                {                    for (m = j; m < n + 1; m++)                    {                        temp = Guass[j, m];                        Guass[j, m] = Guass[k, m];                        Guass[k, m] = temp;                    }                }                if (0 == max)                {                    // "此线性方程为奇异线性方程"                     return x;                }                for (i = j + 1; i < n; i++)                {                    s = Guass[i, j];                    for (m = j; m < n + 1; m++)                    {                        Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);                    }                }            }//结束for (j=0;j<n;j++)            for (i = n - 1; i >= 0; i--)            {                s = 0;                for (j = i + 1; j < n; j++)                {                    s = s + Guass[i, j] * x[j];                }                x[i] = (Guass[i, n] - s) / Guass[i, i];            }            return x;        }    }