最近做的项目有用到非线性拟合这个东西,然后查阅多方面资料,发现大部分人的做法都是通过矩阵的线性求解去完成这个任务。并且,要想通过矩阵求解这个问题,还需要用到一些矩阵的库或者是OpenCV的库。显然这种做法是有点麻烦的。于是,通过多方面的查找,终于在GitHub找到了一个可实现非线性拟合的C++代码,并且不需要下载任何其他的库。
经过测试,该代码拟合效果与excel保持一致,精度要优于excel的非线性拟合。
(无法科学上网了,无法提供源网址了,如有侵权请联系删除!)
以下便是实现非线性拟合的C++代码:
class Matrix {
int N;
double a[100][100];
double b[100];
double x[2][100];
double eps;
public:
void init(int _N, double _a[][100], double _b[], double _x[], double _eps) {
N = _N;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++)
a[i][j] = _a[i][j];
b[i] = _b[i];
x[0][i] = _x[i];
x[1][i] = x[0][i];
}
eps = _eps;
// printf("N : %d, eps : %f\n\n", N, eps);
// for(int i = 0; i < N; i ++)
// for(int j = 0; j < N; j ++)
// printf((j == N - 1)? "%f\n" : "%f ", a[i][j]);
// printf("b:\n");
// for(int i = 0; i < N; i ++) printf((i == N - 1) ? "%f\n" : "%f ", b[i]);
// printf("x:\n");
// for(int i = 0; i < N; i ++) printf((i == N - 1) ? "%f\n" : "%f ", x[0][i]);
}
void showN() {
printf("N = %d\n", N);
}
void showA() {
printf("A:\n");
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
printf((j == N - 1) ? "%.3f\n" : "%.3f ", a[i][j]);
}
void showb() {
printf("b:\n");
for (int i = 0; i < N; i++) printf("%.3f\n", b[i]);
}
void showx(int type) {
// if(type > 1 || type < 0) return printf("Parameter [int type] is invalid in function showx().\n") * 0;
printf("x:\n");
for (int i = 0; i < N; i++) printf("x%d = %.10f\n", i, x[type][i]);
// for(int i = 0; i < N; i++) printf("x%d = %.3f\n", i, x[1][i]);
}
void showEps() {
printf("Epsilon = %f\n", eps);
}
void show() {
printf("--\nN = %d:\n", N);
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) printf("%.3f ", a[i][j]);
printf("%.3f\n", b[i]);
}
}
double getEps() {
return eps;
}
bool finish() {
double R = 0;
for (int i = 0; i < N; i++) R = max(R, abs(x[0][i] - x[1][i]));
return R <= eps;
}
// Cal x, last recursion result store in x[0], new recursion result store in x[1].
void Recursion(int id, int type, double omega) {
// if(type < 0 || type > 2) return printf("Parameter [int type] is invalid in function Recursion().\n") * 0;
if (a[id][id] == 0) {
printf("Div by 0. a[%d][%d] == 0.\n", id, id);
return;
}
double sum = 0;
for (int i = 0; i < N; i++) if (i != id) sum += a[id][i] * x[0][i];
if (type == 0) x[1][id] = (b[id] - sum) / a[id][id];
else if (type == 1) x[0][id] = (b[id] - sum) / a[id][id];
else if (type == 2) x[0][id] = (1 - omega) * x[0][id] + omega * (b[id] - sum) / a[id][id];
//printf("x[%d]:%.3f ", id, x[1][id]);
}
void Jacobi() {
do {
for (int i = 0; i < N; i++) Recursion(i, 0, 0);
for (int i = 0; i < N; i++) swap(x[0][i], x[1][i]);
} while (!finish());
showx(0);
}
void Gauss_Seidel() {
do {
for (int i = 0; i < N; i++) x[1][i] = x[0][i];
for (int i = 0; i < N; i++) Recursion(i, 1, 0);
} while (!finish());
showx(0);
}
double* SOR(double omega) {
do {
for (int i = 0; i < N; i++) x[1][i] = x[0][i];
for (int i = 0; i < N; i++) Recursion(i, 2, omega);
} while (!finish());
return x[0];
// showx(0);
}
};
class LeastSquare {
int n, order;
double x[100], y[100], p[100];
Matrix equations;
public:
void init(double _eps, double X[], double Y[], int ParaNum = 4, int Order = 3)
{
n = ParaNum;//输入数据的个数
for (int i = 0; i < n; i++) x[i] = X[i];
for (int i = 0; i < n; i++) y[i] = Y[i];
for (int i = 0; i < n; i++) p[i] = 1;
order = Order;//阶次
int _N = order;
double _a[100][100], _x[100], _b[100];
for (int i = 0; i < order; i++) _x[i] = 0.0;
//
for (int i = 0; i < order; i++) {
for (int j = 0; j < order; j++) {
_a[i][j] = 0;
for (int _i = 0; _i < n; _i++) _a[i][j] += p[_i] * phi(i, x[_i]) * phi(j, x[_i]);
}
}
for (int i = 0; i < order; i++) {
_b[i] = 0;
for (int j = 0; j < n; j++) _b[i] += p[j] * phi(i, x[j]) * y[j];
}
equations.init(_N, _a, _b, _x, _eps);
}
double* solve() {
double* a = equations.SOR(1.5);
double* Coeff = a;
//printf("\n\nf(x) = ");
//for (int i = 0; i < order; i++) {
// if (i == 0) printf("%.10f", a[i]);
// else {
// if (a[i] > 0) printf(" + %.10f x^%d", a[i], i);
// else printf(" - %.10f x^%d", -a[i], i);
// }//此处输出回归值
//}
//printf("\n\nx\t\tf(x)\t\ty\t\t|f(x) - y|\n");
//for (int i = 0; i < n; i++) {
// double fx = 0;
// for (int j = 0; j < order; j++) fx += a[j] * phi(j, x[i]);
// printf("%f\t%f\t%f\t%f\n", x[i], fx, y[i], fx - y[i] > 0 ? fx - y[i] : y[i] - fx);
//}//输出显示
return Coeff;
}
double phi(int pow, double x) {
double res = 1;
for (int i = 0; i < pow; i++) res *= x;
return res;
}
};
类内的具体作用暂时还没去深入了解(好像是高斯的方法?),感兴趣的朋友可以先行理解一下,然后评论区教一下博主,感谢!!!
以下是我使用上述代码实现非线性拟合的函数代码:
double* NonLinearCal(double X_Point[], double Y_Point[], int Para, int Order)
{
int Order = Order; // 需要返回多少个系数,如果拟合二次项 则有3个系数,线性拟合则2个,三次项4个...
int ParaNum = Para; // 拟合的数据对 有多少对X和Y
double* Temp = NULL;//用来接收solve返回的指针 指向的是拟合的系数
LeastSquare *LS = new LeastSquare;
LS->init(0.000000001, X_Point, Y_Point, ParaNum, Order);
Temp = LS->solve();
double *Coeff = new double[Order];
for(int i=-0;i<Order;i++)
Coeff[i] = temp[i];//temp是局部变量 这个函数结束了,原来的值也就没了,所以需要把这些数据转移到堆上
delete LS;
return Coeff;
}
拿Order = 3为例,此时非线性拟合的方程为一元二次方程,方程可理解为 y = a1 + a2*x +a3 *x^2,那么返回的Coeff[0] 就是a1的值,Coeff[1] 就是a2的值,以此类推。
提供一个计算拟合优度(R^2)的函数:
double R2Cal(double X_Point[], double Y_Point[], double Coeff[], int ParaNum, int Order)
{
double Total = 0;
for (int i = 0; i < ParaNum; i++)
{
Total += Y_Point[i];
}
double Average = Total / ParaNum;
double SSR = 0;
double SST = 0;
//求得平均值
switch (Order)
{
case 2:
{
for (int i = 0; i < ParaNum; i++)
{
SST += pow((Y_Point[i] - Average), 2);
}
for (int i = 0; i < ParaNum; i++)
{
SSR += pow((WlsCoeff[0] + WlsCoeff[1] * X_Point[i] - Average), 2);
}
break;
}
case 3:
{
for (int i = 0; i < ParaNum; i++)
{
SST += pow((Y_Point[i] - Average), 2);
}
for (int i = 0; i < ParaNum; i++)
{
SSR += pow((WlsCoeff[0] + WlsCoeff[1] * X_Point[i] + WlsCoeff[2] * X_Point[i] * X_Point[i] - Average), 2);
}
break;
}
case 4:
{
for (int i = 0; i < ParaNum; i++)
{
SST += pow((Y_Point[i] - Average), 2);
}
for (int i = 0; i < ParaNum; i++)
{
SSR += pow((WlsCoeff[0] + WlsCoeff[1] * X_Point[i] + WlsCoeff[2] * X_Point[i] * X_Point[i] + WlsCoeff[3] * X_Point[i] * X_Point[i] * X_Point[i] - Average), 2);
}
break;
}
default:
break;
}
return SSR / SST;
}