用C语言实现最小二乘法算法



用C语言实现最小二乘法算法

本文博客链接:http://blog.csdn.net/jdh99,作者:jdh,转载请注明.

 

环境：

主机:WIN8

开发环境:MINGW

说明：

参考维基百科最小二乘法资料：

测试文中战列舰例子：用战列舰的长度预测宽度

简单线性模型 y = b₀ + b₁t 的例子[编辑]

随机选定10艘战舰，并分析它们的长度与宽度，寻找它们长度与宽度之间的关系。由下面的描点图可以直观地看出，一艘战舰的长度（t）与宽度（y）基本呈线性关系。散点图如下：

以下图表列出了各战舰的数据，随后步骤是采用最小二乘法确定两变量间的线性关系。

编号	长度 (m)	宽度 (m)	ti - t	yi - y
i	ti	yi	ti*	yi*	ti*yi*	ti*ti*	yi*yi*
1	208	21.6	40.2	3.19	128.238	1616.04	10.1761
2	152	15.5	-15.8	-2.91	45.978	249.64	8.4681
3	113	10.4	-54.8	-8.01	438.948	3003.04	64.1601
4	227	31.0	59.2	12.59	745.328	3504.64	158.5081
5	137	13.0	-30.8	-5.41	166.628	948.64	29.2681
6	238	32.4	70.2	13.99	982.098	4928.04	195.7201
7	178	19.0	10.2	0.59	6.018	104.04	0.3481
8	104	10.4	-63.8	-8.01	511.038	4070.44	64.1601
9	191	19.0	23.2	0.59	13.688	538.24	0.3481
10	130	11.8	-37.8	-6.61	249.858	1428.84	43.6921
总和（Σ）	1678	184.1	0.0	0.00	3287.820	20391.60	574.8490

仿照上面给出的例子

并得到相应的.

然后确定b₁

可以看出，战舰的长度每变化1m，相对应的宽度便要变化16cm。并由下式得到常数项b₀：

在这里随机理论不加阐述。可以看出点的拟合非常好，长度和宽度的相关性大约为96.03％。 利用Matlab得到拟合直

线：

算法结果：

mean_x = 167.800003,mean_y = 18.410000
a = -8.645075,b = 0.161234
300m长度的战舰预测宽度为为39.725140米

源代码：

#include <stdio.h>

/*********************************************************************
*							宏定义
**********************************************************************/

/*********************************************************************
*							数据fifo长度
**********************************************************************/

#define LEN		100

/*********************************************************************
*							数据结构
**********************************************************************/

/*********************************************************************
*							数据单元
**********************************************************************/

struct _Data_Unit
{
	float x;
	float y;
};

/*********************************************************************
*							数据fifo
**********************************************************************/

struct _Data
{
	struct _Data_Unit data[LEN];
	int len;
};

/*********************************************************************
*							全局变量
**********************************************************************/

/*********************************************************************
*							数据fifo
**********************************************************************/

struct _Data Data = 
{
	.len = 0
};

/*********************************************************************
*							函数
**********************************************************************/

/*********************************************************************
*							压入数据
*参数:x:测量数据x
*     y:测量数据y
**********************************************************************/

void push(float x,float y)
{
	int i = 0;
	
	if (Data.len < LEN)
	{
		Data.data[Data.len].x = x;
		Data.data[Data.len++].y = y;
		return;
	}
	
	//数据移动,去掉最后一个数据
	for (i = 0;i < LEN - 1;i++)
	{
		Data.data[i].x = Data.data[i + 1].x;
		Data.data[i].y = Data.data[i + 1].y;
	}
	Data.data[LEN].x = x;
	Data.data[LEN].y = y;
	Data.len = LEN;
}

/*********************************************************************
*							计算估值
*拟合曲线y = a * x + b
*参数:x:需要估值的数的x值
*返回:估值y
**********************************************************************/

float calc(float x)
{
	int i = 0;
	float mean_x = 0;
	float mean_y = 0;
	float num1 = 0;
	float num2 = 0;
	float a = 0;
	float b = 0;
	
	//求t,y的均值
	for (i = 0;i < Data.len;i++)
	{
		mean_x += Data.data[i].x;
		mean_y += Data.data[i].y;
	}
	mean_x /= Data.len;
	mean_y /= Data.len;
	
	printf("mean_x = %f,mean_y = %f\n",mean_x,mean_y);
	
	for (i = 0;i < Data.len;i++)
	{
		num1 += (Data.data[i].x - mean_x) * (Data.data[i].y - mean_y);
		num2 += (Data.data[i].x - mean_x) * (Data.data[i].x - mean_x); 
	}
	
	b = num1 / num2;
	a = mean_y - b * mean_x;
	
	printf("a = %f,b = %f\n",a,b);
	
	return (a + b * x);
}


int main()
{
	float length[10] = {208,152,113,227,137,238,178,104,191,130};
	float width[10] = {21.6,15.5,10.4,31.0,13.0,32.4,19.0,10.4,19.0,11.8};
	int i = 0;
	
	for (i = 0;i < 10;i++)
	{
		push(length[i],width[i]);
	}
	printf("300m长度的战舰预测宽度为为%f米\n",calc(300));
	
	getchar();
	return 0;
}