2乘.txt

是一个用vc编写的非线形最小二乘估计

                m_val[i-1][j-1]=m_val[i-1][j-1]*y;    
                m[i-1][j-1]=m[i-1][j-1]*y;    
            }    
            for (k=1;k<=n;++k)    
                if (k!=i)     
                {    
                    y=m_val[k-1][i-1];    
                    for(j=n;j>=1;--j)    
                    {    
                        m[k-1][j-1]-=m[i-1][j-1]*y;    
                        m_val[k-1][j-1]-=m_val[i-1][j-1]*y;    
                    }    
                }    
         
        }    
        return m;      
    }    
}    
    
CMatrix CMatrix::Transpos()    
{    
    CMatrix m;    
    m.Resize(m_ne,m_nv);    
    for(int i=0;i<m_ne;i++)    
        for(int j=0;j<m_nv;j++)    
            m[i][j]=m_val[j][i];    
    return m;    
}    
    
CMatrix CMatrix::Standard()    
{    
    CMatrix m;    
    m.Resize(m_nv,m_ne);    
    CVector s(m_ne);    
    CVector aver(m_ne);    
    for(int i=0;i<m_ne;i++)    
    {    
        aver[i]=0.0;    
        s[i]=0.0;    
    }    
    for( i=0;i<m_nv;i++)    
        for(int j=0;j<m_ne;j++)    
            aver[j]+=m_val[i][j]/m_nv;    
    for(i=0;i<m_nv;i++)    
        for(int j=0;j<m_ne;j++)    
            s[j]+=(m_val[i][j]-aver[j])*(m_val[i][j]-aver[j]);    
    for(i=0;i<m_ne;i++)    
    {    
        s[i]=sqrt((s[i])/m_nv);    
    }    
    for( i=0;i<m_nv;i++)    
        for(int j=0;j<m_ne;j++)    
            m[i][j]=(m_val[i][j]-aver[j])/s[j];    
    return m;    
}    
    
CMatrix CMatrix::StandardNonLinear()    
{    
    CMatrix M(*this) ;    
    for ( int i = 0; i < M.m_nv; i++ )     
    {    
        for ( int j = 0; j < M.m_nv; j++ )     
        {    
            M[i][j] = M[i][j]/sqrt(M[i][i])/sqrt(M[j][j]) ;    
        }    
    }    
    return M ;    
}    
                
void CMatrix::Unit()    
{    
    if(m_nv!=m_ne)    
        throw("rows!=cols errro");    
    else    
    {    
        int n=m_nv;    
        for(register i=0;i<n;i++)    
            for(register j=0;j<n;j++)    
                m_val[i][j]= (i==j) ? 1.0 : 0.0;    
    }    
    
}    
    
CMatrix CMatrix::CaclCov()    
{    
    register int i,j,k;    
    CMatrix S(m_ne,m_ne);    
    double *x=new double[m_ne];    
    
    memset(x,0,m_ne*sizeof(float));    
    for(i=0;i<m_ne;i++)    
    {    
        for(j=0;j<m_nv;j++)    
            x[i]+=m_val[j][i];    
        x[i]=x[i]/m_nv;    
    }    
    
    for(i=0;i<m_ne;i++)    
    {    
        for(j=0;j<m_ne;j++)    
        {    
            S[i][j]=0;    
            for( k=0;k<m_nv;k++)    
                S[i][j]+=(m_val[k][i]-x[i])*(m_val[k][j]-x[j]);    
            S[i][j]=S[i][j]/(m_nv-1);    
        }    
    }    
    delete[] x;    
    return S;    
}    
    
int CMatrix::GetRows()    
{    
    return m_nv;    
}    
    
int CMatrix::GetCols()    
{    
    return m_ne;    
}    
    
void CMatrix::Resize(int nv,int ne)    
{    
    m_nv=nv;    
    m_ne=ne;    
    delete []m_val;    
    m_val=new CVector[nv];    
    for(int i=0;i<nv;i++)    
        m_val[i].Resize(ne);    
    
}    
    
/*   
CMatrix::CMatrix(CMatrix &amt;M)   
{   
    m_nv=M.m_nv;   
    m_ne=M.m_ne;   
    m_browvector=M.m_browvector;   
    m_val=new CVector[m_nv];   
    for(int i=0;i<m_nv;i++)   
        m_val[i]=M.m_val[i];   
   
}   
*/    
    
    
CMatrix CMatrix::ProduceMatrix( CMatrix B)    
{    
        
    if(m_ne==B.m_nv)    
    {    
        int n=m_nv;    
        int p=B.m_ne;    
        int m=m_ne;    
        CMatrix C(n,p);    
        for(int i=0;i<n;i++)    
            for(int j=0;j<p;j++)    
            {    
                C[i][j]=0;    
                for(int k=0;k<m;k++)    
                    C[i][j]+=m_val[i][k]*B[k][j];    
    
            }    
        return C;    
    }    
    else    
        throw("no match size");    
}    
ifstream &amt; operator >>(ifstream &amt;in,CMatrix &amt;M)    
{    
    int n,e;    
    in>>n;    
    in>>e;    
    M.Resize(n,e);    
    for(int i=0;i<M.m_nv;i++)    
        for(int j=0;j<M.m_ne;j++)    
            in>>M.m_val[i][j];    
    return in;    
}    
ofstream &amt; operator <<(ofstream &amt;out,CMatrix &amt;M)    
{    
    out<<M.m_nv<<" ";    
    out<<M.m_ne<<" ";    
    for(int i=0;i<M.m_nv;i++)    
        for(int j=0;j<M.m_ne;j++)    
            out<<M.m_val[i][j]<<" ";    
    return out;    
}    
ostream &amt; operator <<(ostream &amt;out,CMatrix &amt;M)    
{    
    for(int i=0;i<M.m_nv;i++)    
        for(int j=0;j<M.m_ne;j++)    
            out<<M.m_val[i][j]<<" ";    
    return out;    
}    
void CMatrix::DeleteRow(int r)    
{    
    m_nv--;    
    for(int i=r;i<m_nv;i++)    
        m_val[i]=m_val[i+1];    
}    
    
void CMatrix::DeleteCol(int c)    
{    
    m_ne--;    
    for(int i=0;i<m_nv;i++)    
    {    
        for(int j=c;j<m_ne;j++)    
            m_val[i][j]=m_val[i][j+1];    
    }    
}    
    
CMatrix CMatrix::operator +(const CMatrix &amt;B)    
{    
    CMatrix mat(m_nv,m_ne);    
    CMatrix M=B;    
    if(M.m_nv==m_nv&amt;&amt;M.m_ne==m_ne)    
    {    
        for(int i=0;i<m_nv;i++)    
            for(int j=0;j<m_ne;j++)    
                mat[i][j]=M[i][j]+m_val[i][j];    
        return mat;         
    }    
    else    
        throw("errow");    
}    
CMatrix CMatrix::operator -(const CMatrix &amt;B)    
{    
    CMatrix mat(m_nv,m_ne);    
    CMatrix M=B;    
    if(M.m_nv==m_nv&amt;&amt;M.m_ne==m_ne)    
    {    
        for(int i=0;i<m_nv;i++)    
            for(int j=0;j<m_ne;j++)    
                mat[i][j]=M[i][j]-m_val[i][j];    
        return mat;    
    }    
    else    
        throw("errow");    
}    
    
CVector CMatrix::operator *(const CVector &amt;V)    
{    
    CVector vec(m_nv,false);    
    for(int i=0;i<m_nv;i++)    
        vec[i]=0;    
    CVector B=V;    
    if(B.m_p==m_ne)    
    {    
        for(int i=0;i<m_nv;i++)    
            for(int j=0;j<m_ne;j++)    
                vec[i]+=m_val[i][j]*B[j];    
        return vec;    
    }    
    else     
        throw("errow");    
}    
    
    
void CMatrix::FreeSpace()    
{    
    delete[] m_val;    
}    
    
CVector CMatrix::ColMax()//增广矩阵M    
{    
    int iRow = m_nv ;    
    CVector X( iRow, false ) ;//解向量    
    CMatrix C(iRow,iRow);//消去时所乘的系数    
    CVector* pr=NULL;//存放r行首地址    
    CVector* pk=NULL;//存放k行首地址    
    
    for(int k=0;k<iRow;k++)    
    {    
        if(m_val[k][k]==0)    
        {    
            cout<<"主对角线元素为0"<<endl;    
        }    
        //选取矩阵Ａ最大值    
        int r=k;//记录列主元的行数    
        double temp=0;//存放最大值    
        for(int i1=k;i1<iRow;i1++)    
        {    
            if(fabs(m_val[i1][k])>temp)    
            {    
                temp=fabs(m_val[i1][k]);    
                r=i1;    
            }    
        }    
    
        //如果当前行主对角线元素不是最大，则交换增广矩阵Ａ的r行k行    
        if(r!=k)    
        {    
            pr=&amt;m_val[r];    
            pk=&amt;m_val[k];    
            exchange(m_val[r],m_val[k],iRow+1);    
        }    
        for(int i=k+1;i<iRow;i++)    
        {    
            C[i][k]=m_val[i][k]/m_val[k][k];    
            m_val[i][iRow]=m_val[i][iRow]-C[i][k]*m_val[k][iRow];    
            for(int j=k+1;j<iRow;j++)    
            {    
                m_val[i][j]=m_val[i][j]-C[i][k]*m_val[k][j];//化为三角型矩阵    
            }    
        }    
    }    
//  cout<<A[0][0]<<endl<<A[0][1]<<endl<<A[0][2]<<endl<<A[1][0]<<endl<<A[1][1]<<endl<<A[1][2]<<endl;    
        
    //回代求解    
    X[iRow-1]=m_val[iRow-1][iRow]/m_val[iRow-1][iRow-1];//最后一个x值    
    //cout<<"X"<<iRow<<"="<<X[iRow-1]<<endl;    
    for(int l=iRow-2;l>=0;l--)//回代求解    
    {    
        double c=0;    
        for(int p=l+1;p<iRow;p++)    
        {    
            c=c+m_val[l][p]*X[p];    
        }    
        X[l]=(m_val[l][iRow]-c)/m_val[l][l];    
        //cout<<"X"<<l+1<<"="<<X[l]<<endl;    
    }    
    return X;    
}    
    
void CMatrix::exchange(CVector&amt; pr,CVector&amt; pk,int col)         //交换增广矩阵Ａ的r行k行的函数，pr：r行首地址，pk：k行首地址    
{    
    double temp;    
    for(int i=0;i<col;i++)    
    {    
        temp=pr[i];    
        pr[i]=pk[i];    
        pk[i]=temp;    
    }    
}    
<SCRIPT src="/inc/gg_read2.js"></SCRIPT>  
// Matrix.h: interface for the CMatrix class.    
//    
//////////////////////////////////////////////////////////////////////    
    
#if !defined(AFX_MATRIX_H__8FC17A77_A602_40FB_A1D1_05958784C8E6__INCLUDED_)    
#define AFX_MATRIX_H__8FC17A77_A602_40FB_A1D1_05958784C8E6__INCLUDED_    
    
#if _MSC_VER > 1000    
#pragma once    
#endif // _MSC_VER > 1000    
    
#include "Matrix.h"    
#include "Vector.h"    
#include <fstream.h>    
class CVector;    
class CMatrix      
{    
public:    
    CMatrix StandardNonLinear();    
    void exchange(CVector&amt; pr,CVector&amt; pk,int col) ;        //交换增广矩阵Ａ的r行k行的函数，pr：r行首地址，pk：k行首地址    
    CVector ColMax();    
    void FreeSpace();    
    CMatrix(int rows=1,int cols=1,bool browvector=true);    
    CMatrix(const CMatrix &amt;M);    
    void DeleteCol(int c);    
    void DeleteRow(int r);    
    CMatrix ProduceMatrix( CMatrix b);    
    void Resize(int nv,int ne);    
    int GetCols();    
    int GetRows();    
    CMatrix CaclCov();//协差阵    
    void Unit();//单位矩阵    
    CMatrix Inverse();// 求逆    
    CMatrix Transpos();//求转置    
    CMatrix Standard();// 标准化    
    CMatrix operator +(const CMatrix &amt;M);    
    CMatrix operator -(const CMatrix &amt;M);    
    CVector operator *(const CVector &amt;V);    
    CMatrix operator *( const double &amt;n);    
    CMatrix &amt;operator =(const CMatrix &amt;M);    
    CMatrix &amt;operator =(const CVector &amt;M);    
    CMatrix operator *(const CMatrix &amt;M);    
    CVector &amt;operator [](int index);    
    virtual ~CMatrix();    
    friend ifstream &amt; operator >>(ifstream &amt;in,CMatrix &amt;M);    
    friend ofstream &amt; operator <<(ofstream &amt;out,CMatrix &amt;M);    
    friend ostream &amt; operator <<(ostream &amt;out,CMatrix &amt;M);    
    friend class CPls;    
    friend class CDosc;    
    friend class CVector;