zoukankan      html  css  js  c++  java
  • 用矩阵运算实现最小二乘法曲线拟合算法

    1. 多项式拟合函数:

    y= a0 + a1x + a2x^2 + ... + akx^k    (其中k为拟合次数)

    当k=1 为线性拟合 ,k=2 为二次多项式 ...  三次多项式。

    2. 最小二乘原理矩阵算法原理:

    X*A=Y
    A=((X'*X)-1)*X'*Y

    |1  X1  X1^2  ... X1^k|            |y0|
    |1  X2  X2^2  ... X2^k| |a0|     |y1|
    |...                              | |a1|     |.  |
    |...                              | |.   |  = |.  |
    |...                              | |ak|     |.  |
    |1  Xn  Xn^2  ... Xn^k|            |yn|

    其中 X为初始范德蒙矩阵,A 为系数向量, Y 为因变量值向量

    3.计算单元算法:几个矩阵运算算法函数

    1 ) 矩阵的转置运算; 2)矩阵的逆运算;3 )矩阵乘法运算

    4. C++代码实现

    1)矩阵乘法运算函数:

    BOOL CMatrix::Mul(const Matrix &a, const Matrix &b, Matrix &c)
    {
        if(a.n!=b.m)
        {
            return FALSE;
        }
        
        c.m = a.m;
        c.n = b.n;
        int i, j, k;
        for (i=0; i<a.m; i++)
        {
            for (j=0; j<b.n; j++)
            {
                for (c.p[i * c.n + j] = 0, k =0; k<a.n; k++)
                {
                    c.p[i * c.n + j] += a.p[i * a.n + k] * b.p[k * b.n + j];
                }
            }
        }
        return TRUE;
    }
    2)矩阵转置运算函数:
    void CMatrix::Trs(Matrix &a, Matrix &trs)
    {
        trs.m = a.n;
        trs.n = a.m;
        
        for (int i = 0; i < a.m; i++)
        {
            for (int j = 0; j < a.n; j++)
            {
                trs.p[j * a.m + i] = a.p[i * a.n + j];
            }
        }
    }

    3)矩阵逆运算函数:

    // 工具函数

    long double CMatrix::Det(Matrix &a)
    {
        long double det = 0;
        if (a.m != a.n)
        {   
            //...
            return 0;
        }
        if (a.n == 1)
        {
            det = a.p[0];
            return det;
        }
        else
        {
            for (int i = 0; i < a.n; i++)
            {
                if (i % 2 == 0)
                {
                    Matrix mat;
                    Adjunct(a, i, 0,mat);
                    det += a.p[i * a.n] * Det(mat);
                    delete mat.p;
                }
                else  
                {
                    Matrix mat;
                    Adjunct(a, i, 0,mat);
                    det -= a.p[i * a.n] * Det(mat);
                    delete mat.p;
                }
            }
        }
     
        return det;
    }
     

    //工具函数
    void CMatrix::Adjunct(Matrix a, int indexm, int indexn,Matrix &adj)
    {
        adj.SetSize(a.n - 1, a.n - 1);
        for (int i = 0; i < indexm; i++)
        {
            for (int j = 0; j < indexn; j++)
            {
                adj.p[i * (a.n - 1) + j] = a.p[i * a.n + j];
            }
            for (int k = indexn + 1; k < a.n; k++)
            {
                adj.p[i *(a.n - 1) + k -1] = a.p[i * a.n + k];
            }
        }
     
        for (int m = indexm + 1; m < a.n; m++)
        {
            for (int j = 0; j < a.n - 1; j++)
            {
                adj.p[(m - 1) * (a.n - 1) + j] = a.p[m * a.n + j];
            }
            for (int k = indexn + 1; k < a.n; k++)
            {
                adj.p[(m - 1) * (a.n - 1) + k - 1] = a.p[m * a.n + k];
            }
        }
     
    }
     
    // 逆运算函数
    BOOL CMatrix::Inv(Matrix &a,Matrix &inv)
    {
        Matrix Temp(a.m,a.n);
     
        if(a.m!=a.n)
        {
            return FALSE;
        }
     
     
        long double det = Det(a);
        if (det == 0)
        {
            return FALSE;
        }
        
        for (int i = 0; i < Temp.m; i++)
        {
            for (int j = 0; j < Temp.n; j++)
            {
                if ((i +j) % 2 == 0)
                {    
                    Matrix mat;        
                    Adjunct(a, i, j,mat);
                    Temp.p[i * Temp.m + j] = Det(mat) / det;
                    delete mat.p;
                }
                else
                {
                    Matrix mat;
                    Adjunct(a, i, j,mat);
                    Temp.p[i * Temp.m + j] = -Det(mat) / det;
                    delete mat.p;
                }
            }
        }
                
        Trs(Temp,inv);    
        return TRUE;
    }

    4)多项式拟合函数可以根据以上运算单元,结合矩阵运算公式:A=((X'*X)-1)*X'*Y

    自由实现。

    5)数据结构定义:

    struct Matrix{  
        int m, n;   
        long double *p; 
        Matrix()
        {
            m = 0;
            n = 0;
            p = NULL;
        }
        Matrix(int t_m ,int t_n)
        {
            m = t_m;
            n = t_n;
            p = new long double[m*n];
        }
        void SetSize(int t_m,int t_n)
        {
            m = t_m;
            n = t_n;
            p = new long double[m*n];
        }
        ~Matrix()
        {
            if(p != NULL)
            {
            //    delete p;
            }
        }
    };

  • 相关阅读:
    BEC listen and translation exercise 44
    中译英12
    BEC listen and translation exercise 43
    中译英11
    BEC listen and translation exercise 42
    中译英10
    BEC listen and translation exercise 41
    中译英9
    BEC listen and translation exercise 40
    中译英8
  • 原文地址:https://www.cnblogs.com/Esperanto/p/5647918.html
Copyright © 2011-2022 走看看