本文研究通過C語言實(shí)現(xiàn)最小二乘法擬合直線。
1 引入
最小二乘法，簡單來說就是根據(jù)一組觀測得到的數(shù)值，尋找一個函數(shù)，使得函數(shù)與觀測點(diǎn)的誤差的平方和達(dá)到最小。在工程實(shí)踐中，這個函數(shù)通常是比較簡單的，例如一次函數(shù)或二次函數(shù)。
汽車上的毫米波雷達(dá)可以探測到其他目標(biāo)車輛，通過最小二乘法擬合目標(biāo)車輛歷史點(diǎn)，可以簡單地預(yù)測目標(biāo)汽車未來的走向。
后文會推導(dǎo)最小二乘法擬合直線，并通過C語言實(shí)現(xiàn)，最后進(jìn)行簡單的驗(yàn)證。對于二次及更高次多項(xiàng)式的擬合，采用類似的方法。
2 公式推導(dǎo)
作為工程應(yīng)用，推導(dǎo)公式不需要像數(shù)學(xué)上的那么抽象，只需要針對當(dāng)前需求推導(dǎo)即可。
首先，需要擬合的方程為一次函數(shù)：
y = a x + b y = ax + b y=ax+b
并且已知n個觀測點(diǎn)：
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
.
.
.
(
x
n
,
y
n
)
(x_{1}, y_{1}),(x_{2}, y_{2}),...(x_{n}, y_{n})
(x1?,y1?),(x2?,y2?),...(xn?,yn?)
則擬合的誤差的平方和為：
E
(
a
,
b
)
=
∑
i
=
1
n
(
a
x
i
+
b
?
y
i
)
2
E(a,b)=\sum_{i=1}^n\left({{}}a{{x}}{}_{{i}}+b-{{y}}{}_{{i}}\right)^2
E(a,b)=i=1∑n?(axi?+b?yi?)2
注意，x和y是已知量，該函數(shù)是關(guān)于a和b的二元函數(shù)。目標(biāo)是求誤差的最小值，因此需要分別對a和b求偏導(dǎo)數(shù)：
?
E
?
a
=
∑
i
=
1
n
2
?
(
a
x
i
+
b
?
y
i
)
?
x
i
=
2
∑
i
=
1
n
(
a
x
i
2
+
b
x
i
?
y
i
x
i
)
\frac{\partial E}{\partial a}=\sum_{i=1}^n{2\cdot\left(a{{x}}{}_{{i}}+b-{{y}}{}_{{i}}\right)}\cdot{{x}}{}_{{i}}=2\sum_{i=1}^n{\left(a{{x}}_{{i}}^{{2}}+b{{x}}{}_{{i}}-{{y}}{}_{{i}}{{x}}{}_{{i}}\right)}
?a?E?=i=1∑n?2?(axi?+b?yi?)?xi?=2i=1∑n?(axi2?+bxi??yi?xi?)
?
?
E
?
b
=
∑
i
=
1
n
2
?
(
a
x
i
+
b
?
y
i
)
=
2
∑
i
=
1
n
(
a
x
i
+
b
?
y
i
)
\:\frac{\partial E}{\partial b}=\sum_{i=1}^n{2\cdot\left(a{{x}}{}_{{i}}+b-{{y}}{}_{{i}}\right)}=2\sum_{i=1}^n\left(a{{x}}{}_{{i}}+b-{{y}}{}_{{i}}\right)
?b?E?=i=1∑n?2?(axi?+b?yi?)=2i=1∑n?(axi?+b?yi?)
對偏導(dǎo)數(shù)取值為0，可以得到線性方程組：
a
∑
i
=
1
n
x
i
2
?
+
b
∑
i
=
1
n
x
i
?
=
∑
i
=
1
n
y
i
x
i
a\sum_{i=1}^n{{{x}}_{{i}}^{{2}}}\:+b\sum_{i=1}^n{{x}}{}_{{i}}\:=\sum_{i=1}^n{{y}}{}_{{i}}{{x}}{}_{{i}}
ai=1∑n?xi2?+bi=1∑n?xi?=i=1∑n?yi?xi?
a
∑
i
=
1
n
x
i
?
+
?
b
n
?
=
∑
i
=
1
n
y
i
a\sum_{i=1}^n{{x}}{}_{{i}}\:+\:bn_{{}}\:=\sum_{i=1}^n{{y}}{}_{{i}}{{}}
ai=1∑n?xi?+bn?=i=1∑n?yi?
求解方程組可以用克拉默法則：
D
=
?
∣
∑
i
=
1
n
x
i
2
∑
i
=
1
n
x
i
∑
i
=
1
n
x
i
n
∣
?
=
?
n
∑
i
=
1
n
x
i
2
?
?
?
(
∑
i
=
1
n
x
i
)
2
D=\:\left|\begin{matrix}{\sum_{i=1}^n{{{x}}_{{i}}^{{2}}}{}} & {\sum_{i=1}^n{{{x}}_{{i}}^{{}}}{}} \\ {\sum_{i=1}^n{{{x}}_{{i}}^{{}}}{}} & n\end{matrix}\right|\:=\:n\sum_{i=1}^n{{x}}_{{i}}^{{2}}\:-\:\left(\sum_{i=1}^n{{x}}_{{i}}^{}\right)^2
D=
?∑i=1n?xi2?∑i=1n?xi??∑i=1n?xi?n?
?=ni=1∑n?xi2??(i=1∑n?xi?)2
D
a
=
?
∣
∑
i
=
1
n
x
i
y
i
∑
i
=
1
n
x
i
∑
i
=
1
n
y
i
n
∣
?
=
?
n
∑
i
=
1
n
x
i
y
i
?
?
?
∑
i
=
1
n
x
i
∑
i
=
1
n
y
i
{{D}}{}_{{a}}=\:\left|\begin{matrix}{\sum_{i=1}^n{{{{x{}}{}_{{i}}y}}_{{i}}}{}} & {\sum_{i=1}^n{{{x}}_{{i}}^{{}}}{}} \\ {\sum_{i=1}^n{{y{}}_{{i}}^{{}}}{}} & n\end{matrix}\right|\:=\:n\sum_{i=1}^n{{{x{}}{}_{{i}}y}}_{{i}}^{{}}\:-\:\sum_{i=1}^n{x{}}{}_{{i}}\sum_{i=1}^n{y{}}_{{i}}^{{}}
Da?=
?∑i=1n?xi?yi?∑i=1n?yi??∑i=1n?xi?n?
?=ni=1∑n?xi?yi??i=1∑n?xi?i=1∑n?yi?
D
b
=
?
∣
∑
i
=
1
n
x
i
2
∑
i
=
1
n
x
i
y
i
∑
i
=
1
n
x
i
∑
i
=
1
n
y
i
∣
?
=
?
∑
i
=
1
n
x
i
2
?
∑
i
=
1
n
y
i
?
?
?
∑
i
=
1
n
x
i
∑
i
=
1
n
x
i
y
i
{{D}}{}_{}=\:\left|\begin{matrix}{\sum_{i=1}^n{{{x}}_{{i}}^{{2}}}{}} & {\sum_{i=1}^n{{{{x{}}{}_{{i}}y}}_{{i}}^{{}}}{}} \\ {\sum_{i=1}^n{{{x}}_{{i}}^{{}}}{}} & \sum_{i=1}^n{{y}}{}_{{i}}\end{matrix}\right|\:=\:\sum_{i=1}^n{{x}}_{{i}}^{{2}}\:\sum_{i=1}^n{{y}}{}_{{i}}\:-\:\sum_{i=1}^n{{x}}_{{i}}^{{}}\sum_{i=1}^n{{{x{}}{}_{{i}}{y{}}{}_{{i}}}}
Db?=
?∑i=1n?xi2?∑i=1n?xi??∑i=1n?xi?yi?∑i=1n?yi??
?=i=1∑n?xi2?i=1∑n?yi??i=1∑n?xi?i=1∑n?xi?yi?
可以求得a和b：
a
?
=
?
D
a
D
?
=
?
n
∑
x
i
y
i
?
∑
x
i
∑
y
i
n
∑
x
i
2
?
(
∑
x
i
)
2
a\:=\:\frac{{{D}}{}_{{a}}}{D}\:=\:\frac{n\sum{{{x}}{}_{{i}}}{{y}}{}_{{i}}-\sum{{{x}}{}_{{i}}}\sum{{y{}}{}_{{i}}}}{n\sum{{x}}_{{i}}^{{2}}-\left(\sum{{x}}{}_{{i}}\right)^2}
a=DDa??=n∑xi2??(∑xi?)2n∑xi?yi??∑xi?∑yi??
b
?
=
?
D
b
D
?
=
?
∑
x
i
2
∑
y
i
?
∑
x
i
∑
x
i
y
i
n
∑
x
i
2
?
(
∑
x
i
)
2
b\:=\:\frac{{{D}}{}_{b{}}}{D}\:=\:\frac{\sum{{x}}_{{i}}^{{2}}{{}}{}_{{}}\sum{{y{}}{}_{{i}}}-\sum{{{x}}{}_{{i}}}\sum{{{{x}}{}_{{i}}y{}}{}_{{i}}}}{n\sum{{x}}_{{i}}^{{2}}-\left(\sum{{x}}{}_{{i}}\right)^2}
b=DDb??=n∑xi2??(∑xi?)2∑xi2??∑yi??∑xi?∑xi?yi??
3 C語言代碼實(shí)現(xiàn)
1）首先設(shè)計(jì)頭文件polyfit_types.h，用于定義一些基本類型和結(jié)構(gòu)體；
//polyfit_types.h
#ifndef POLYFIT_TYPES_H
#define POLYFIT_TYPES_H
typedef signed char int8;
typedef unsigned char uint8;
typedef short int16;
typedef unsigned short uint16;
typedef int int32;
typedef unsigned int uint32;
typedef float float32;
typedef struct POINT_TAG
{
float32 x_f32;
float32 y_f32;
}POINT_TYPE;
typedef struct LINE_TAG
{
float32 a_f32;
float32 b_f32;
}LINE_TYPE;
#endif
這里先定義一些基本類型，比如uint8，float32等。接著定義兩個結(jié)構(gòu)體：點(diǎn)和直線。點(diǎn)結(jié)構(gòu)體的成員是點(diǎn)的x和y坐標(biāo)，直線結(jié)構(gòu)體的成員是直線的系數(shù)a和b。
2）設(shè)計(jì)頭文件polyfit.h；
//polyfit.h
#ifndef POLYFIT_H
#define POLYFIT_H
//Include
#include "polyfit_types.h"
//Define
#define EPS 0.000001F
#define OK 1U
#define NOK 0U
//Function
uint8 polyfit(const POINT_TYPE* points_vs, const uint8 point_number_u8, LINE_TYPE* line_s, float32* square_error_f32);
#endif
頭文件首先包含了type頭文件。接著是宏定義EPS表示一個極小值，用于浮點(diǎn)數(shù)相等比較。然后就是函數(shù)polyfit的聲明，函數(shù)的前兩個參數(shù)是輸入的點(diǎn)的數(shù)組和點(diǎn)的數(shù)量，后兩個參數(shù)表示輸出的直線和殘差的平方和。
3）設(shè)計(jì)C文件polyfit.c，是實(shí)現(xiàn)算法的核心代碼；
//polyfit.c
#include <math.h>
#include <string.h>
#include "polyfit.h"
//Function
uint8 polyfit(const POINT_TYPE* points_vs, const uint8 point_number_u8, LINE_TYPE* line_s, float32* square_error_f32)
{
//var
float32 sum_x_f32 = 0.0F;
float32 sum_xy_f32 = 0.0F;
float32 sum_x2_f32 = 0.0F;
float32 sum_y_f32 = 0.0F;
float32 xi = 0.0F;
float32 yi = 0.0F;
float32 Det_f32 = 0.0F;
float32 Det_a_f32 = 0.0F;
float32 Det_b_f32 = 0.0F;
float32 Err = 0.0F;
uint8 index_u8;
//initialize
*square_error_f32 = 0.0F;
memset(line_s, 0.0F, sizeof(LINE_TYPE));
//calculate sum
for (index_u8 = 0U; index_u8 < point_number_u8; index_u8++)
{
xi = points_vs[index_u8].x_f32;
yi = points_vs[index_u8].y_f32;
sum_x_f32 += xi;
sum_xy_f32 += xi * yi;
sum_x2_f32 += xi * xi;
sum_y_f32 += yi;
}
//calculate determinant
Det_f32 = point_number_u8 * sum_x2_f32 - sum_x_f32 * sum_x_f32;
Det_a_f32 = point_number_u8 * sum_xy_f32 - sum_x_f32 * sum_y_f32;
Det_b_f32 = sum_x2_f32 * sum_y_f32 - sum_x_f32 * sum_xy_f32;
//determine if Det_f32 is 0
if (fabsf(Det_f32) < EPS)
{
return NOK;
}
//calculate coefficient
line_s->a_f32 = Det_a_f32 / Det_f32;
line_s->b_f32 = Det_b_f32 / Det_f32;
//calculate sum of square error
for (index_u8 = 0U; index_u8 < point_number_u8; index_u8++)
{
xi = points_vs[index_u8].x_f32;
yi = points_vs[index_u8].y_f32;
Err = yi - (line_s->a_f32 * xi + line_s->b_f32);
*square_error_f32 += Err * Err;
}
return OK;
}
代碼中，首先定義局部變量，并初始化輸出參數(shù)。接著，按照公式計(jì)算三個行列式，并判斷D是否為0，如果為零就停止計(jì)算并返回。最后通過克拉默法則計(jì)算系數(shù)，以及根據(jù)算出的直線計(jì)算殘差的平方和。
4 測試驗(yàn)證
在主函數(shù)里可以簡單寫一段測試代碼，驗(yàn)證一下是否正確。
1）參考某百科，假設(shè)輸入4個點(diǎn)為（1，6），（2，5），（3，7），（4，10），測試代碼如下；
#include <stdio.h>
#include "polyfit.h"
int main()
{
POINT_TYPE point_vs[4];
point_vs[0].x_f32 = 1.0F; point_vs[0].y_f32 = 6.0F;
point_vs[1].x_f32 = 2.0F; point_vs[1].y_f32 = 5.0F;
point_vs[2].x_f32 = 3.0F; point_vs[2].y_f32 = 7.0F;
point_vs[3].x_f32 = 4.0F; point_vs[3].y_f32 = 10.0F;
LINE_TYPE line_s;
float32 square_error_f32;
uint8 retVal;
retVal = polyfit(point_vs, 4U, &line_s, &square_error_f32);
printf("retVal = %d , a = %f , b = %f , err = %f \r\n", retVal, line_s.a_f32, line_s.b_f32, square_error_f32);
}
打印出來的結(jié)果為：
這表示擬合的直線方程為：
y = 1.4 x + 3.5 y = 1.4x + 3.5 y=1.4x+3.5
繪圖出來的結(jié)果是：
2）如果輸入的4個點(diǎn)是（1，6），（1，5），（1，7），（1，10），即四個點(diǎn)在一條垂直于x軸的直線上，這時(shí)行列式D=0，就無法得出擬合結(jié)果：
#include <stdio.h>
#include "polyfit.h"
int main()
{
POINT_TYPE point_vs[4];
point_vs[0].x_f32 = 1.0F; point_vs[0].y_f32 = 6.0F;
point_vs[1].x_f32 = 1.0F; point_vs[1].y_f32 = 5.0F;
point_vs[2].x_f32 = 1.0F; point_vs[2].y_f32 = 7.0F;
point_vs[3].x_f32 = 1.0F; point_vs[3].y_f32 = 10.0F;
LINE_TYPE line_s;
float32 square_error_f32;
uint8 retVal;
retVal = polyfit(point_vs, 4U, &line_s, &square_error_f32);
printf("a = %f , b = %f , err = %f \r\n", line_s.a_f32, line_s.b_f32, square_error_f32);
}
打印出來的結(jié)果為：
5 總結(jié)
本文本文研究通過C語言實(shí)現(xiàn)最小二乘法擬合直線。在工程應(yīng)用中，一次和二次多項(xiàng)式的擬合用的比較多。二次多項(xiàng)式擬合可以參考一次的推導(dǎo)過程和編程過程，需要求解三階行列式求解三個系數(shù)。文章來源：http://www.zghlxwxcb.cn/news/detail-651888.html
>>返回個人博客總目錄文章來源地址http://www.zghlxwxcb.cn/news/detail-651888.html
到了這里，關(guān)于C語言編程：最小二乘法擬合直線的文章就介紹完了。如果您還想了解更多內(nèi)容，請?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
