最小二乘法拟合曲线：二次函数

void myLMS_poly2(const std::vector<double> src_x, const std::vector<double> src_y, int size, std::vector<double>& dst_y)
{
double a, b, c;


//Mat A = Mat_<double>(3, 3);
//Mat B = Mat_<double>(3, 1);
//Mat C = Mat_<double>(3, 1);


Mat A(3, 3, CV_64F);
Mat B(3, 1, CV_64F);
Mat X(3, 1, CV_64F);


double sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += pow(src_x[i], 4);
}
A.at<double>(0, 0) = sum;


sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += pow(src_x[i], 3);
}
A.at<double>(0, 1) = A.at<double>(1, 0) = sum;


sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += pow(src_x[i], 2);
}
A.at<double>(0, 2) = A.at<double>(1, 1) = A.at<double>(2, 0) = sum;


sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += src_x[i];
}
A.at<double>(1, 2) = A.at<double>(2, 1) = sum;


sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += 1;
}
A.at<double>(2, 2) = sum;


sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += pow(src_x[i], 2) * src_y[i];
}
B.at<double>(0, 0) = sum;


sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += src_x[i] * src_y[i];
}
B.at<double>(1, 0) = sum;


sum = 0;
for (size_t i = 0; i < size; ++i)
{
sum += src_y[i];
}
B.at<double>(2, 0) = sum;


//X = A.inv() * B;


solve(A, B, X);// A * X = B


a = X.at<double>(0, 0);
b = X.at<double>(1, 0);
c = X.at<double>(2, 0);


printf("\n%f, %f, %f\n", a, b, c);


for (size_t i = 0; i < size; ++i)
{
double ret = a * pow(src_x[i], 2) + b * src_x[i] + c;
dst_y.push_back(ret);
}
}