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

               


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


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

 

环境:

主机:WIN8

开发环境:MINGW


说明:

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


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

简单线性模型 y = b0 + b1t 的例子[编辑]

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

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

编号长度 (m)宽度 (m)ti - tyi - y   
itiyiti*yi*ti*yi*ti*ti*yi*yi*
120821.640.23.19128.2381616.0410.1761
215215.5-15.8-2.9145.978249.648.4681
311310.4-54.8-8.01438.9483003.0464.1601
422731.059.212.59745.3283504.64158.5081
513713.0-30.8-5.41166.628948.6429.2681
623832.470.213.99982.0984928.04195.7201
717819.010.20.596.018104.040.3481
810410.4-63.8-8.01511.0384070.4464.1601
919119.023.20.5913.688538.240.3481
1013011.8-37.8-6.61249.8581428.8443.6921
总和(Σ)1678184.10.00.003287.82020391.60574.8490


仿照上面给出的例子

\bar t = \frac {\sum_{i=1}^n t_i}{n} = \frac {1678}{10} = 167{.}8 并得到相应的\bar y = 18{.}41.

然后确定b1

b_1 = \frac{\sum_{i=1}^n (t_i- \bar {t})(y_i - \bar y)}{\sum_{i=1}^n (t_i- \bar t)^2}
= \frac{3287{.}820} {20391{.}60} = 0{.}1612 \;,

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

b_0 = \bar y - b_1 \bar t = 18{.}41 - 0{.}1612 \cdot 167{.}8 = -8{.}6394 \;,

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

线:

最小二乘法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 = 0float mean_x = 0float mean_y = 0float num1 = 0float num2 = 0float a = 0float 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;}

           
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值