语言版的线性回归分析函数

1、C语言版的线性回归分析函数 分类： C/C+ 2007-08-03 23:39 13840人阅读 评论(31) 收藏 举报 语言c数学计算delphisystem        前几天，清理出一些十年以前DOS下的程序及代码，看来目前也没什么用了，想打个包刻在光碟上，却发现有些代码现在可能还能起作用，其中就有计算一元回归和多元回归的代码，一看代码文件时间，居然是1993年的，于是稍作整理，存放在这，分析虽不十分完整，但一般应用是没问题的，最起码，可提供给那些刚学C的学生们参考。先看看一元线性回归函数代码： /

2、60;求线性回归方程：Y = a + bx/ dadarows*2数组：X, Y；rows：数据行数；a, b：返回回归系数/ SquarePoor4：返回方差分析指标: 回归平方和，剩余平方和，回归平方差，剩余平方差/ 返回值：0求解成功，-1错误int LinearRegression(double *data, int rows, double *a, double *b, double *

3、SquarePoor)    int m;    double *p, Lxx = 0.0, Lxy = 0.0, xa = 0.0, ya = 0.0;    if (data = 0 | a = 0 | b = 0 

4、| rows < 1)        return -1;    for (p = data, m = 0; m < rows; m +)            xa += *p

5、 +;        ya += *p +;        xa /= rows;                         

6、;            / X平均值    ya /= rows;                            

7、60;        / Y平均值    for (p = data, m = 0; m < rows; m +, p += 2)            Lxx += (*p 

8、;- xa) * (*p - xa);             / Lxx = Sum(X - Xa)平方)        Lxy += (*p - xa) * (*(p + 1) - ya

9、);       / Lxy = Sum(X - Xa)(Y - Ya)        *b = Lxy / Lxx;                  

10、60;              / b = Lxy / Lxx    *a = ya - *b * xa;                 

11、;             / a = Ya - b*Xa    if (SquarePoor = 0)        return 0;    / 方差分析    Sq

12、uarePoor0 = SquarePoor1 = 0.0;    for (p = data, m = 0; m < rows; m +, p +)            Lxy = *a + *b * 

13、;*p +;        SquarePoor0 += (Lxy - ya) * (Lxy - ya); / U(回归平方和)        SquarePoor1 += (*p - Lxy) * (*p - Lxy); / Q

14、(剩余平方和)        SquarePoor2 = SquarePoor0;                  / 回归方差    SquarePoor3 = SquarePoor1 / (rows - 

15、2);     / 剩余方差    return 0; 为了理解代码，把几个与代码有关的公式写在下面（回归理论和公式推导就免了，网上搜索到处是，下面的公式图片也是网上搜的，有些公式图形网上没找到或者不合适，可参见后面多元回归中的公式）：1、回归方程式：2、回归系数：       其中：            

16、       3、回归平方和：4、剩余平方和：实例计算：double data1122 = /    X      Y    187.1, 25.4,    179.5, 22.8,    157.0, 20.6,    

17、;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, 

18、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(&q

19、uot;回归方程式:    Y = %.5lf", Answer0);    for (i = 1; i < cols; i +)        printf(" + %.5lf*X%d", Answeri, i);    

20、printf(" ");    printf("回归显著性检验: ");    printf("回归平方和：%12.4lf  回归方差：%12.4lf ", SquarePoor0, SquarePoor2);    printf("剩余平方和：%12.4lf  剩余方差：%12.4lf ", SquarePoor1, 

21、;SquarePoor3);    printf("离差平方和：%12.4lf  标准误差：%12.4lf ", SquarePoor0 + SquarePoor1, sqrt(SquarePoor3);    printf("F   检  验：%12.4lf  相关系数：%12.4lf ", SquarePoor2 /Squa

22、rePoor3,           sqrt(SquarePoor0 / (SquarePoor0 + SquarePoor1);    printf("剩余分析: ");    printf("      观察值     &#

23、160;估计值      剩余值    剩余平方 ");    for (i = 0, p = dat; i < rows; i +, p +)            v = 

24、;Answer0;        for (j = 1; j < cols; j +, p +)            v += *p * Answerj;        p

25、rintf("%12.2lf%12.2lf%12.2lf%12.2lf ", *p, v, *p - v, (*p - v) * (*p - v);        system("pause");int main()    double Answer2, SquarePoor4; 

26、;   if (LinearRegression(double*)data1, 12, &Answer0, &Answer1, SquarePoor) = 0)        Display(double*)data1, Answer, SquarePoor, 12, 2);    return 0; 

27、   运行结果： 上 面的函数和例子程序不仅计算了回归方程式，还计算了显著性检验指标，例如F检验指标，我们可以在统计F分布表上查到 F0.01(1,10)=10.04（注：括号里的1,10分别为回归平方和和剩余平方和所拥有的自由度），小于计算的F检验值25.94，可以认为该回 归例子高度显著。如果使用图形界面，可以根据原始数据和计算结果绘制各种图表，如散点图、趋势图、控制图等。很多非线性方程可以借助数学计算，转化为直线方程进行回归分析。同一元线性回归相比，多元线性回归分析代码可就复杂多了，必须求解线性方程，因此本代码中包含一个可独立使用的线性方程求解函数：

28、0;void FreeData(double *dat, double *d, int count)    int i, j;    free(d);    for (i = 0; i < count; i +)        fre

29、e(dati);    free(dat);/ 解线性方程。datacount*(count+1)矩阵数组；count：方程元数；/ Answercount：求解数组 。返回：0求解成功，-1无解或者无穷解int LinearEquations(double *data, int count, double *Answer)    int j, m, n;   

30、60;double tmp, *dat, *d = data;    dat = (double*)malloc(count * sizeof(double*);    for (m = 0; m < count; m +, d += (count + 1)   

31、         datm = (double*)malloc(count + 1) * sizeof(double);        memcpy(datm, d, (count + 1) * sizeof(double);       

32、 d = (double*)malloc(count + 1) * sizeof(double);    for (m = 0; m < count - 1; m +)            / 如果主对角线元素为0，行交换  

33、60;     for (n = m + 1; n < count && datmm = 0.0; n +)                    if ( datnm

34、 != 0.0)                            memcpy(d, datm, (count + 1) * sizeof(double);      

35、60;         memcpy(datm, datn, (count + 1) * sizeof(double);                memcpy(datn, d, (count + 1) * sizeof(do

36、uble);                            / 行交换后，主对角线元素仍然为0，无解，返回-1        if (datmm = 0.0)   

37、;                 FreeData(dat, d, count);            return -1;            

38、0;   / 消元        for (n = m + 1; n < count; n +)                    tmp = dat

39、nm / datmm;            for (j = m; j <= count; j +)                datnj -= tmp * datmj

40、;                for (j = 0; j < count; j +)        dj = 0.0;    / 求得count - 1的元  

41、  Answercount - 1 = datcount - 1count / datcount - 1count - 1;    / 逐行代入求各元    for (m = count - 2; m >= 0; m -)   

42、60;        for (j = count - 1; j > m; j -)            dm += Answerj * datmj;        Ans

43、werm = (datmcount - dm) / datmm;        FreeData(dat, d, count);    return 0;/ 求多元回归方程：Y = B0 + B1X1 + B2X2 + .BnXn/ datarows*cols二维数组；X1i,X2i,.X

44、ni,Yi (i=0 to rows-1)/ rows：数据行数；cols数据列数；Answercols：返回回归系数数组(B0,B1.Bn)/ SquarePoor4：返回方差分析指标: 回归平方和，剩余平方和，回归平方差，剩余平方差/ 返回值：0求解成功，-1错误int MultipleRegression(double *data, int rows, int cols, double *Answer, double *

45、SquarePoor)    int m, n, i, count = cols - 1;    double *dat, *p, a, b;    if (data = 0 | Answer = 0 | rows < 2 |

46、0;cols < 2)        return -1;    dat = (double*)malloc(cols * (cols + 1) * sizeof(double);    dat0 = (double)rows;    for (n

47、60;= 0; n < count; n +)                     / n = 0 to cols - 2          

48、0; a = b = 0.0;        for (p = data + n, m = 0; m < rows; m +, p += cols)             

49、       a += *p;            b += (*p * *p);                datn + 1 = a;

50、0;                             / dat0, n+1 = Sum(Xn)        dat(n + 1) * (c

51、ols + 1) = a;               / datn+1, 0 = Sum(Xn)        dat(n + 1) * (cols + 1) + n + 1 =

52、 b;       / datn+1,n+1 = Sum(Xn * Xn)        for (i = n + 1; i < count; i +)           

53、  / i = n+1 to cols - 2                    for (a = 0.0, p = data, m = 0; m < rows; m&#

54、160;+, p += cols)                a += (pn * pi);            dat(n + 1) * (cols + 1) + 

55、;i + 1 = a;   / datn+1, i+1 = Sum(Xn * Xi)            dat(i + 1) * (cols + 1) + n + 1 = a;   / 

56、dati+1, n+1 = Sum(Xn * Xi)                for (b = 0.0, m = 0, p = data + n; m < rows; m +, p += 

57、cols)        b += *p;    datcols = b;                              

58、60;    / dat0, cols = Sum(Y)    for (n = 0; n < count; n +)            for (a = 0.0, p = data, m 

59、= 0; m < rows; m +, p += cols)            a += (pn * pcount);        dat(n + 1) * (cols + 1) +

60、60;cols = a;        / datn+1, cols = Sum(Xn * Y)        n = LinearEquations(dat, cols, Answer);          /

61、0;计算方程式    / 方差分析    if (n = 0 && SquarePoor)            b = b / rows;            

62、                    / b = Y的平均值        SquarePoor0 = SquarePoor1 = 0.0;        p 

63、;= data;        for (m = 0; m < rows; m +, p +)                    for (i = 1, a&#

64、160;= Answer0; i < cols; i +, p +)                a += (*p * Answeri);              

65、60;/ a = Ym的估计值            SquarePoor0 += (a - b) * (a - b);    / U(回归平方和)            SquarePoor1

66、60;+= (*p - a) * (*p - a);  / Q(剩余平方和)(*p = Ym)                SquarePoor2 = SquarePoor0 / count;       /

67、 回归方差  if (rows - cols > 0.0)    SquarePoor3 = SquarePoor1 / (rows - cols); / 剩余方差  else    SquarePoor3 = 0.0;        free(dat);    return n;为了理解代码，同样贴几个主要公式在下面，其中回归平方和和剩余平方和公式和

68、一元回归相同：1、回归方程式：, 2、回归系数方程组：3、F检验：4、相关系数：，其中，Syy是离差平方和（回归平方和与剩余平方和之和）。该公式其实就是U/(U+Q)的平方根（没找到这个公式的图）。5、回归方差：U / m，m为回归方程式中自变量的个数（没找到图）。6、剩余方差：Q / (n - m - 1)，n为观察数据的样本数，m同上（没找到图）。7、标准误差：也叫标准误，就是剩余方差的平方根（没找到图）。下面是一个多元回归的例子：double data155 = /   X1   X2 &

69、#160;  X3   X4    Y   316, 1536, 874, 981, 3894 ,   385, 1771, 777, 1386, 4628 ,   299, 1565, 678, 1672, 4569 ,   326,&#

70、160;1970, 785, 1864, 5340 ,   441, 1890, 785, 2143, 5449 ,   460, 2050, 709, 2176, 5599 ,   470, 1873, 673, 1769, 5010 ,   504, 1955, 7

71、93, 2207, 5694 ,   348, 2016, 968, 2251, 5792 ,   400, 2199, 944, 2390, 6126 ,   496, 1328, 749, 2287, 5025 ,   497, 1920, 952, 2388,&#

72、160;5924 ,   533, 1400, 1452, 2093, 5657 ,   506, 1612, 1587, 2083, 6019 ,   458, 1613, 1485, 2390, 6141 ,;void Display(double *dat, double *Answer, d

73、ouble *SquarePoor, int rows, int cols)    double v, *p;    int i, j;    printf("回归方程式:    Y = %.5lf", Answer0);    for (i =&#

74、160;1; i < cols; i +)        printf(" + %.5lf*X%d", Answeri, i);    printf(" ");    printf("回归显著性检验: ");    printf("回归平方和：%12.4lf  回归方差：%12.4lf ", SquarePoor0, SquarePoor2);    printf("剩余平方和：%12.4lf  剩余方差：%12.4lf ",