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


    引用

    http://blog.sina.com.cn/s/blog_6e51df7f0100thie.html

    对代码稍作修改和注释,防止链接失效。

      1         ///<summary>
      2         ///用最小二乘法拟合二元多次曲线
      3         ///例如y=ax+b
      4         ///其中MultiLine将返回a,b两个参数。
      5         ///a对应MultiLine[1]
      6         ///b对应MultiLine[0]
      7         ///</summary>
      8         ///<param name="arrX">已知点的x坐标集合</param>
      9         ///<param name="arrY">已知点的y坐标集合</param>
     10         ///<param name="length">已知点的个数</param>
     11         ///<param name="dimension">方程的最高次数</param>
     12         public  double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线
     13         {
     14             int n = dimension + 1;                  //dimension次方程需要求 dimension+1个 系数
     15             double[,] Guass = new double[n, n + 1];      //高斯矩阵 例如:y=a0+a1*x+a2*x*x
     16             for (int i = 0; i < n; i++)
     17             {
     18                 int j;
     19                 for (j = 0; j < n; j++)
     20                 {
     21                     Guass[i, j] = SumArr(arrX, j + i, length);
     22                 }
     23                 Guass[i, j] = SumArr(arrX, i, arrY, 1, length);
     24             }
     25 
     26             return ComputGauss(Guass, n);
     27 
     28         }
     29         private  double SumArr(double[] arr, int n, int length) //求数组的元素的n次方的和
     30         {
     31             double s = 0;
     32             for (int i = 0; i < length; i++)
     33             {
     34                 if (arr[i] != 0 || n != 0)
     35                     s = s + Math.Pow(arr[i], n);
     36                 else
     37                     s = s + 1;
     38             }
     39             return s;
     40         }
     41         private  double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)
     42         {
     43             double s = 0;
     44             for (int i = 0; i < length; i++)
     45             {
     46                 if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
     47                     s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);
     48                 else
     49                     s = s + 1;
     50             }
     51             return s;
     52 
     53         }
     54         private  double[] ComputGauss(double[,] Guass, int n)
     55         {
     56             int i, j;
     57             int k, m;
     58             double temp;
     59             double max;
     60             double s;
     61             double[] x = new double[n];
     62 
     63             for (i = 0; i < n; i++) x[i] = 0.0;//初始化
     64 
     65 
     66             for (j = 0; j < n; j++)
     67             {
     68                 max = 0;
     69 
     70                 k = j;
     71                 for (i = j; i < n; i++)
     72                 {
     73                     if (Math.Abs(Guass[i, j]) > max)
     74                     {
     75                         max = Guass[i, j];
     76                         k = i;
     77                     }
     78                 }
     79 
     80 
     81 
     82                 if (k != j)
     83                 {
     84                     for (m = j; m < n + 1; m++)
     85                     {
     86                         temp = Guass[j, m];
     87                         Guass[j, m] = Guass[k, m];
     88                         Guass[k, m] = temp;
     89 
     90                     }
     91                 }
     92 
     93                 if (0 == max)
     94                 {
     95                     // "此线性方程为奇异线性方程" 
     96 
     97                     return x;
     98                 }
     99 
    100 
    101                 for (i = j + 1; i < n; i++)
    102                 {
    103 
    104                     s = Guass[i, j];
    105                     for (m = j; m < n + 1; m++)
    106                     {
    107                         Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);
    108 
    109                     }
    110                 }
    111 
    112 
    113             }//结束for (j=0;j<n;j++)
    114 
    115 
    116             for (i = n - 1; i >= 0; i--)
    117             {
    118                 s = 0;
    119                 for (j = i + 1; j < n; j++)
    120                 {
    121                     s = s + Guass[i, j] * x[j];
    122                 }
    123 
    124                 x[i] = (Guass[i, n] - s) / Guass[i, i];
    125 
    126             }
    127 
    128             return x;
    129         }//返回值是函数的系数

    例如:y=a0+a1*x 返回值则为a0 a1

    例如:y=a0+a1*x+a2*x*x 返回值则为a0 a1 a2

  • 相关阅读:
    python定义函数的三种形式
    python函数的返回值
    python函数的调用
    python函数的定义
    python文件操作
    Python2和3字符编码的区别
    python的字符编码
    python异常处理
    python深浅拷贝
    python色彩缤纷的python(改变字体颜色以及样式)
  • 原文地址:https://www.cnblogs.com/crhdyl/p/5320046.html
Copyright © 2020-2023  润新知