最小二乘法求直线y=kx+b的参数

#include "iostream.h" #include "math.h" typedef struct Point { int x; int y; } stPoint; bool FittingLine(stPoint* pp , int pointNum , double* k , double* b); void main() { stPoint p[10]; // p[0].x = 120 , p[0].y = 2; // p[1].x = 119 , p[1].y = 3; // p[2].x = 119 , p[2].y = 4; // p[3].x = 119 , p[3].y = 5; // p[4].x = 119 , p[4].y = 6; // p[5].x = 118 , p[5].y = 7; // p[6].x = 118 , p[6].y = 8; // p[7].x = 118 , p[7].y = 9; // p[8].x = 118 , p[8].y = 10; // p[9].x = 117 , p[9].y = 11; p[0].x = 2 , p[0].y = 120; p[1].x = 3 , p[1].y = 119; p[2].x = 4 , p[2].y = 119; p[3].x = 5 , p[3].y = 119; p[4].x = 6 , p[4].y = 119; p[5].x = 7 , p[5].y = 118; p[6].x = 8 , p[6].y = 118; p[7].x = 9 , p[7].y = 118; p[8].x = 10 , p[8].y = 118; p[9].x = 11 , p[9].y = 117; // p[0].x = -0 , p[0].y = 2; // p[1].x = -0 , p[1].y = 3; // p[2].x = -0 , p[2].y = 4; // p[3].x = -0 , p[3].y = 5; // p[4].x = -0 , p[4].y = 6; // p[5].x = -0 , p[5].y = 7; // p[6].x = -0 , p[6].y = 8; // p[7].x = -0 , p[7].y = 9; // p[8].x = -0 , p[8].y = 10; // p[9].x = -0 , p[9].y = 11; double k = 0; double b = 0; bool flag = FittingLine(p , 10 , &k , &b); if(flag) cout << "k = " << k << " , b = " << b << endl; else cout << "k不存在" << endl; } bool FittingLine(stPoint* pp , int pointNum , double* k , double* b) { /************************************** * 最小二乘法求直线方程参数 * pp: [in] 点数组首地址 * pointNum: [in] 点的个数 * k,b: [out] y = kx + b中的两个参数 * * 返回值: * false:垂直于x轴的直线 k、b无效 * true: k、b有效 * * 2011年6月8日 * * 数学原理: * * M = 累加 [yi - (k * xi + b)] ^ 2 * * 求使M取到最小的k 和 b * * 分别对k ，b求偏导，另其 = 0，得到如下的二元一次方程组 * * sumXY = k * sumX2 + b * sumX * sumY = k * sumX + b * n (n即为pointNum) * * 若解不唯一，则为垂直于x轴的直线 返回false * 否则，k b 即为对应直线方程的参数 返回true 及 k、b的值 **************************************/ double sumXY = 0.0; // x * y 的累加 double sumX2 = 0.0; // x^2 的累加 double sumX = 0.0; // x 的累加 double sumY = 0.0; // y 的累加 // 计算出方程组的系数 int i = 0; for(i = 0 ; i < pointNum ; i++) { sumXY += pp[i].x * pp[i].y; sumX2 += pp[i].x * pp[i].x; sumX += pp[i].x; sumY += pp[i].y; } // 解方程组 // L1 = k * Rk1 + b * Rb1 // L2 = k * Rk2 + b * Rb2 double L1 = sumXY * pointNum; double Rk1 = sumX2 * pointNum; double L2 = sumY * sumX; double Rk2 = sumX * sumX; // Rb1 = sumX * n // Rb2 = n * sumX // 两者相减后被消去 double flag = Rk1 - Rk2; if(fabs(flag) > 0.00001) { *k = (L1 - L2) / flag; *b = (sumY - *k * sumX) / pointNum; return true; } else return false; }