关闭

最小二乘法C#源码转

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

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:22889次
    • 积分:698
    • 等级:
    • 排名:千里之外
    • 原创:38篇
    • 转载:47篇
    • 译文:2篇
    • 评论:6条
    文章分类
    最新评论