求多元线性回归的算法

能够解决形如y=a0+a1*x1+a2*x2+a3*x3的算法

如若用最小二乘法将怎样推导？

最好有公式或原代码（C++）.


拜托各位高手！急用！！！

4 个解决方案

#1

/*
多元线性回归分析
*/
#include <iostream>
#include <fstream>
#include <math>
#include "matrix.h"

////////////////////////////////

// Note: The following conditional compilation statements are included

//       so that you can (likely to) compile this sample, without any

//       modification, using a compiler which does not support any or

//       all of the ANSI C++ features like NAMEPACE, TEMPLATE, and

//       EXCEPTION, used in this class.

//

//       If you have a compiler, such as C++ Builder, Borland C++ 5.0,

//       MS Visual C++ 5.0 or higher, etc., which supports most of the ANSI

//       C++ features used in this class, you do not need to include these

//       statements in your program.

//
#ifndef _NO_NAMESPACE
using namespace         std;
using namespace         math;
#define STD std
#else
#define STD
#endif
#ifndef _NO_TEMPLATE
typedef matrix<double>  Matrix;
#else
typedef matrix          Matrix;
#endif
#ifndef _NO_EXCEPTION
#define TRYBEGIN() \
    try \
    {
#define CATCHERROR() \
} \
 \
/* */ \
catch(const STD::exception &e) \
{ \
    cerr << "Error: " << e.what() << endl; \
}

#else
#define TRYBEGIN()
#define CATCHERROR()
#endif
static Matrix           m_x, m_y, m_a;

/*
从数据文件中读取数据，格式如下
第一行是观测点的个数n
第二行是自变量的个数m
然后是n行数据x，自变量x(0)～x(m-1)
然后是n行应变量y
*/
long read_data(char *datfile)
{
    long    n, m, i, j;
    fstream fs(datfile);
    fs >> n;
    fs >> m;

    //cout<<"n="<<n<<",m="<< m<<endl;
    m_x.SetSize(n, m);  //row=n, col=m
    m_y.SetSize(n, 1);  //row=n, col=1
    fs >> m_x;
    fs >> m_y;

    //cout << "x=" << endl << m_x << endl;

    //cout << "y=" << endl << m_y << endl;
    return n;
}

/* */
long solve(void)
{
    Matrix  c, cct, cy;
    long    x, y;
    x=m_x.ColNo();
    y=m_x.RowNo();
    m_a.SetSize(x + 1, 1);
    c.SetSize(x + 1, y);
    for(long i=0; i < x; i++)
    {
        for(long j=0; j < y; j++)
        {
            c(i, j)=m_x(j, i);
        }
    }

    for(long j=0; j < y; j++)
    {
        c(x, j)=1;
    }

    //cout << "c=" << endl << c << endl;
    cct=c *~c;
    cy=c * m_y;

    //cct*a=cy
    m_a=cct.Solve(cy);

    //cout << "m_a=" << endl << m_a << endl;
    return 0;
}

/* */
void eval(void)
{
    double  q, s, r, u, y, t, ye, ax;
    Matrix  v, qm;
    long    n, m;
    m=m_x.ColNo();
    n=m_x.RowNo();
    v.SetSize(m, 1);
    qm.SetSize(m, 1);
    q=0;
    for(long i=0; i < n; i++)
    {
        y=0;
        for(long j=0; j < m; j++)
        {
            y+=m_a(j, 0) * m_x(i, j);
        }

        y+=m_a(m, 0);
        q+=(m_y(i, 0) - y) * (m_y(i, 0) - y);
    }

    cout << "偏差平方和, q=" << q << endl;  //偏差平方和
    s=sqrt(q / n);
    cout << "平均标准差, s=" << s << endl;  //平均标准偏差
    ye=0;
    for(long i=0; i < n; i++)
    {
        ye+=m_y(i, 0) / n;
    }

    t=0;
    for(long i=0; i < n; i++)
    {
        t+=(m_y(i, 0) - ye) * (m_y(i, 0) - ye);
    }

    r=sqrt(1 - q / t);
    cout << "复相关系数, r=" << r << endl;  //复相关系数
    u=0;
    for(long i=0; i < n; i++)
    {
        y=0;
        for(long j=0; j < m; j++)
        {
            y+=m_a(j, 0) * m_x(i, j);
        }

        y+=m_a(m, 0);
        u+=(ye - y) * (ye - y);
    }

    cout << "回归平方和, u=" << u << endl;  //回归平方和
    for(long i=0; i < m; i++)
    {
        qm(i, 0)=0;
        for(long j=0; j < n; j++)
        {
            ax=0;
            for(long k=0; k < m; k++)
            {
                if(k != j)
                {
                    ax+=m_a(k, 0) * m_x(j, k);
                }
            }

            qm(i, 0)+=(m_y(i, 0) - (m_a(m, 0) + ax)) * (m_y(i, 0) - (m_a(m, 0) + ax));
        }

        v(i, 0)=sqrt(1 - q / qm(i, 0));
    }

    cout << "偏相关系数, v=" << endl << v << endl;  //偏相关系数
}

/* */
int main(int argc, char **argv)
{
    if(argc != 2)
    {
        cout << "mls <datafile name>" << endl;
        return 0;
    }

    TRYBEGIN()  //try
    read_data(argv[1]);
    solve();
    cout << "回归系数, a=" << endl << m_a << endl;  //回归系数
    eval();
    CATCHERROR()    //catch
    return 0;
}

程序中的Matrix是矩阵类。
是不是这种东西？

#2

其实上网搜索一下就出来了。或者用MATLAB也行。

#3

先谢谢 NowCan(城市浪人) 
我调试后在结贴。

我不想用Matlab

#4

我有完整的代码，需要可以发给你