前几天，清理出一些十年以前
DOS
下的程序及代码，看来目前也没什么用了，想打个包刻在光碟上，却发现有些代码现在可能还能起作用，其中就有计算一元回归和多元回归的代码，一看代码文件时间，居然是
1993
年的，于是稍作整理，存放在这，分析虽不十分完整，但一般应用是没问题的，最起码，可提供给那些刚学
C
的学生们参考。
先看看一元线性回归函数代码：
//
求线性回归方程：Y=a+bx
//
dada[rows*2]数组：X,Y；rows：数据行数；a,b：返回回归系数
//
SquarePoor[4]：返回方差分析指标:回归平方和，剩余平方和，回归平方差，剩余平方差
//
返回值：0求解成功，-1错误
int
LinearRegression(
double
*
data,
int
rows,
double
*
a,
double
*
b,
double
*
SquarePoor)
{
int
m;
double
*
p,Lxx
=
0.0
,Lxy
=
0.0
,xa
=
0.0
,ya
=
0.0
;
if
(data
==
0
||
a
==
0
||
b
==
0
||
rows
<
1
)
return
-
1
;
for
(p
=
data,m
=
0
;m
<
rows;m
++
)
{
xa
+=
*
p
++
;
ya
+=
*
p
++
;
}
xa
/=
rows;
//
X平均值
ya
/=
rows;
//
Y平均值
for
(p
=
data,m
=
0
;m
<
rows;m
++
,p
+=
2
)
{
Lxx
+=
((
*
p
-
xa)
*
(
*
p
-
xa));
//
Lxx=Sum((X-Xa)平方)
Lxy
+=
((
*
p
-
xa)
*
(
*
(p
+
1
)
-
ya));
//
Lxy=Sum((X-Xa)(Y-Ya))
}
*
b
=
Lxy
/
Lxx;
//
b=Lxy/Lxx
*
a
=
ya
-
*
b
*
xa;
//
a=Ya-b*Xa
if
(SquarePoor
==
0
)
return
0
;
//
方差分析
SquarePoor[
0
]
=
SquarePoor[
1
]
=
0.0
;
for
(p
=
data,m
=
0
;m
<
rows;m
++
,p
++
)
{
Lxy
=
*
a
+
*
b
*
*
p
++
;
SquarePoor[
0
]
+=
((Lxy
-
ya)
*
(Lxy
-
ya));
//
U(回归平方和)
SquarePoor[
1
]
+=
((
*
p
-
Lxy)
*
(
*
p
-
Lxy));
//
Q(剩余平方和)
}
SquarePoor[
2
]
=
SquarePoor[
0
];
//
回归方差
SquarePoor[
3
]
=
SquarePoor[
1
]
/
(rows
-
2
);
//
剩余方差
return
0
;
}
为了理解代码，把几个与代码有关的公式写在下面(回归理论和公式推导就免了，网上搜索到处是，下面的公式图片也是网上搜的，有些公式图形网上没找到或者不合适，可参见后面多元回归中的公式)：
1、回归方程式：
2、回归系数：
其中：
3、回归平方和：
4、剩余平方和：
实例计算：
double
data1[
12
][
2
]
=
{
//
XY
{
187.1
,
25.4
},
{
179.5
,
22.8
},
{
157.0
,
20.6
},
{
197.0
,
21.8
},
{
239.4
,
32.4
},
{
217.8
,
24.4
},
{
227.1
,
29.3
},
{
233.4
,
27.9
},
{
242.0
,
27.8
},
{
251.9
,
34.2
},
{
230.0
,
29.2
},
{
271.8
,
30.0
}
};
void
Display(
double
*
dat,
double
*
Answer,
double
*
SquarePoor,
int
rows,
int
cols)
{
double
v,
*
p;
int
i,j;
printf(
"
回归方程式: Y=%.5lf
"
,Answer[
0
]);
for
(i
=
1
;i
<
cols;i
++
)
printf(
"
+%.5lf*X%d
"
,Answer[i],i);
printf(
"
"
);
printf(
"
回归显著性检验:
"
);
printf(
"
回归平方和：%12.4lf回归方差：%12.4lf
"
,SquarePoor[
0
],SquarePoor[
2
]);
printf(
"
剩余平方和：%12.4lf剩余方差：%12.4lf
"
,SquarePoor[
1
],SquarePoor[
3
]);
printf(
"
离差平方和：%12.4lf标准误差：%12.4lf
"
,SquarePoor[
0
]
+
SquarePoor[
1
],sqrt(SquarePoor[
3
]));
printf(
"
F检验：%12.4lf相关系数：%12.4lf
"
,SquarePoor[
2
]
/
SquarePoor[
3
],
sqrt(SquarePoor[
0
]
/
(SquarePoor[
0
]
+
SquarePoor[
1
])));
printf(
"
剩余分析:
"
);
printf(
"
观察值估计值剩余值剩余平方
"
);
for
(i
=
0
,p
=
dat;i
<
rows;i
++
,p
++
)
{
v
=
Answer[
0
];
for
(j
=
1
;j
<
cols;j
++
,p
++
)
v
+=
*
p
*
Answer[j];
printf(
"
%12.2lf%12.2lf%12.2lf%12.2lf
"
,
*
p,v,
*
p
-
v,(
*
p
-
v)
*
(
*
p
-
v));
}
system(
"
pause
"
);
}
int
main()
{
double
Answer[2
],SquarePoor[
4
];
if
(LinearRegression((
double
*
)data1,
12
,
&
Answer[
0
],
&
Answer[
1
],SquarePoor)
==
0
)
Display((
double
*
)data1,Answer,SquarePoor,
12
,
2
);
return
0
;
}
运行结果：
上面的函数和例子程序不仅计算了回归方程式，还计算了显著性检验指标，例如F检验指标，我们可以在统计F分布表上查到F0.01(1,10)=10.04(注：括号里的1,10分别为回归平方和和剩余平方和所拥有的自由度)，小于计算的F检验值25.94，可以认为该回归例子高度显著。
如果使用图形界面，可以根据原始数据和计算结果绘制各种图表，如散点图、趋势图、控制图等。很多非线性方程可以借助数学计算，转化为直线方程进行回归分析。
同一元线性回归相比，多元线性回归分析代码可就复杂多了，必须求解线性方程，因此本代码中包含一个可独立使用的线性方程求解函数：
void
FreeData(
double
**
dat,
double
*
d,
int
count)
{
int
i,j;
free(d);
for
(i
=
0
;i
<
count;i
++
)
free(dat[i]);
free(dat);
}
//
解线性方程。data[count*(count+1)]矩阵数组；count：方程元数；
//
Answer[count]：求解数组。返回：0求解成功，-1无解或者无穷解
int
LinearEquations(
double
*
data,
int
count,
double
*
Answer)
{
int
j,m,n;
double
tmp,
**
dat,
*
d
=
data;
dat
=
(
double
**
)malloc(count
*
sizeof
(
double
*
));
for
(m
=
0
;m
<
count;m
++
,d
+=
(count
+
1
))
{
dat[m]
=
(
double
*
)malloc((count
+
1
)
*
sizeof
(
double
));
memcpy(dat[m],d,(count
+
1
)
*
sizeof
(
double
));
}
d
=
(
double
*
)malloc((count
+
1
)
*
sizeof
(
double
));
for
(m
=
0
;m
<
count
-
1
;m
++
)
{
//
如果主对角线元素为0，行交换
for
(n
=
m
+
1
;n
<
count
&&
dat[m][m]
==
0.0
;n
++
)
{
if
(dat[n][m]
!=
0.0
)
{
memcpy(d,dat[m],(count
+
1
)
*
sizeof
(
double
));
memcpy(dat[m],dat[n],(count
+
1
)
*
sizeof
(
double
));
memcpy(dat[n],d,(count
+
1
)
*
sizeof
(
double
));
}
}
//
行交换后，主对角线元素仍然为0，无解，返回-1
if
(dat[m][m]
==
0.0
)
{
FreeData(dat,d,count);
return
-
1
;
}
//
消元
for
(n
=
m
+
1
;n
<
count;n
++
)
{
tmp
=
dat[n][m]
/
dat[m][m];
for
(j
=
m;j
<=
count;j
++
)
dat[n][j]
-=
tmp
*
dat[m][j];
}
}
for
(j
=
0
;j
<
count;j
++
)
d[j]
=
0.0
;
//
求得count-1的元
Answer[count
-
1
]
=
dat[count
-
1
][count]
/
dat[count
-
1
][count
-
1
];
//
逐行代入求各元
for
(m
=
count
-
2
;m
>=
0
;m
--
)
{
for
(j
=
count
-
1
;j
>
m;j
--
)
d[m]
+=
Answer[j]
*
dat[m][j];
Answer[m]
=
(dat[m][count]
-
d[m])
/
dat[m][m];
}
FreeData(dat,d,count);
return
0
;
}
//
求多元回归方程：Y=B0+B1X1+B2X2+...BnXn
//
data[rows*cols]二维数组；X1i,X2i,...Xni,Yi(i=0torows-1)
//
rows：数据行数；cols数据列数；Answer[cols]：返回回归系数数组(B0,B1...Bn)
//
SquarePoor[4]：返回方差分析指标:回归平方和，剩余平方和，回归平方差，剩余平方差
//
返回值：0求解成功，-1错误
int
MultipleRegression(
double
*
data,
int
rows,
int
cols,
double
*
Answer,
double
*
SquarePoor)
{
int
m,n,i,count
=
cols
-
1
;
double
*
dat,
*
p,a,b;
if
(data
==
0
||
Answer
==
0
||
rows
<
2
||
cols
<
2
)
return
-
1
;
dat
=
(
double
*
)malloc(cols
*
(cols
+
1
)
*
sizeof
(
double
));
dat[
0
]
=
(
double
)rows;
for
(n
=
0
;n
<
count;n
++
)
//
n=0tocols-2
{
a
=
b
=
0.0
;
for
(p
=
data
+
n,m
=
0
;m
<
rows;m
++
,p
+=
cols)
{
a
+=
*
p;
b
+=
(
*
p
*
*
p);
}
dat[n
+
1
]
=
a;
//
dat[0,n+1]=Sum(Xn)
dat[(n
+
1
)
*
(cols
+
1
)]
=
a;
//
dat[n+1,0]=Sum(Xn)
dat[(n
+
1
)
*
(cols
+
1
)
+
n
+
1
]
=
b;
//
dat[n+1,n+1]=Sum(Xn*Xn)
for
(i
=
n
+
1
;i
<
count;i
++
)
//
i=n+1tocols-2
{
for
(a
=
0.0
,p
=
data,m
=
0
;m
<
rows;m
++
,p
+=
cols)
a
+=
(p[n]
*
p[i]);
dat[(n
+
1
)
*
(cols
+
1
)
+
i
+
1
]
=
a;
//
dat[n+1,i+1]=Sum(Xn*Xi)
dat[(i
+
1
)
*
(cols
+
1
)
+
n
+
1
]
=
a;
//
dat[i+1,n+1]=Sum(Xn*Xi)
}
}
for
(b
=
0.0
,m
=
0
,p
=
data
+
n;m
<
rows;m
++
,p
+=
cols)
b
+=
*
p;
dat[cols]
=
b;
//
dat[0,cols]=Sum(Y)
for
(n
=
0
;n
<
count;n
++
)
{
for
(a
=
0.0
,p
=
data,m
=
0
;m
<
rows;m
++
,p
+=
cols)
a
+=
(p[n]
*
p[count]);
dat[(n
+
1
)
*
(cols
+
1
)
+
cols]
=
a;
//
dat[n+1,cols]=Sum(Xn*Y)
}
n
=
LinearEquations(dat,cols,Answer);
//
计算方程式
//
方差分析
if
(n
==
0
&&
SquarePoor)
{
b
=
b
/
rows;
//
b=Y的平均值
SquarePoor[
0
]
=
SquarePoor[
1
]
=
0.0
;
p
=
data;
for
(m
=
0
;m
<
rows;m
++
,p
++
)
{
for
(i
=
1
,a
=
Answer[
0
];i
<
cols;i
++
,p
++
)
a
+=
(
*
p
*
Answer[i]);
//
a=Ym的估计值
SquarePoor[
0
]
+=
((a
-
b)
*
(a
-
b));
//
U(回归平方和)
SquarePoor[
1
]
+=
((
*
p
-
a)
*
(
*
p
-
a));
//
Q(剩余平方和)(*p=Ym)
}
SquarePoor[
2
]
=
SquarePoor[
0
]
/
count;
//
回归方差
if (rows - cols > 0.0)
SquarePoor[3] = SquarePoor[1] / (rows - cols); // 剩余方差
else
SquarePoor[3] = 0.0;
}
free(dat);
return
n;
}
为了理解代码，同样贴几个主要公式在下面，其中回归平方和和剩余平方和公式和一元回归相同：
1、回归方程式：,
2、回归系数方程组：
3、F检验：
4、相关系数：，其中，Syy是离差平方和(回归平方和与剩余平方和之和)。该公式其实就是U/(U+Q)的平方根(没找到这个公式的图)。
5、回归方差：U / m，m为回归方程式中自变量的个数(没找到图)。
6、剩余方差：Q / (n - m - 1)，n为观察数据的样本数，m同上(没找到图)。
7、标准误差：也叫标准误，就是剩余方差的平方根(没找到图)。
下面是一个多元回归的例子：
double
data[
15
][
5
]
=
{
//
X1X2X3X4Y
{
316
,
1536
,
874
,
981
,
3894
},
{
385
,
1771
,
777
,
1386
,
4628
},
{
299
,
1565
,
678
,
1672
,
4569
},
{
326
,
![京东抢购助手_[Windows] 京东极速抢购助手V2.0,支持京东健康+扫货抢购](https://img-blog.csdnimg.cn/img_convert/6833e7cbc4b6d54f2a9a429dbfa5b10b.png)
