2乘.txt

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

  
#include <iostream.h>    
#include <math.h>    
#include "Matrix.h"    
    
double fun(double t, double m, double a, double&amt; dm, double&amt; da)    
{    
    double  fenzi = 282 ;    
    
    double  fenmu_1 = 1 ;    
    double  fenmu_2 = m*exp(0-a*t) ;    
    double  fenmu = 1+fenmu_2 ;    
    dm = 0-(1/(fenmu*fenmu))*282*exp(0-a*t) ;    
    da = 1/(fenmu*fenmu)* 282*m*t*exp(0-a*t) ;    
//  cout<<fenzi/fenmu<<endl;    
    
    return  fenzi / fenmu ;    
}    
    
#define M  2//方程的维数（系数矩阵的维数）    
#define F  9    
void main()    
{    
    double w = 0.5 ; //当ｄ>100时的乘数ｗ值    
    //自变量t, 函数值y,　参数b，ｂ0，残差ｅ，残差平方和Ｑ，偏导数ｄＢ，系数矩阵Ａ,校正值Ｘ    
    double t[F] = { 11, 21, 31, 41, 51, 61, 66, 69, 72 } ;    
    double y[F] = { 20.5, 47.1, 83.5, 171, 226.7, 259.7, 271.3, 275.3, 280.3 } ;    
    double b[M] = { 80, 0.2 };    
    double b0[M] = { 80, 0.2 };    
    double b1[M] ;    
    //double y[9] = { 14.75057909,42.00247589,100.6422118,179.8121533,239.1396291,266.913152,273.3240223,275.8071105,277.5909959 } ;    
    //double b[2] = { 64.4764473660417, 0.1154 };    
    //double b[2] = { 64.47645, 0.1154 };    
        
//  double t[F] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31 } ;    
//  double y[F] = { 0.67,0.85,1.28,1.75,2.27,2.75,3.69,4.71,6.36,7.73,9.91,12.75,16.55,20.1,27.35,32.55,37.55,44.75,53.38,71.61,83.89,97.46,112.73,135.12,153.6,160.32,167.05,174.9,177.87,180.19,180.79} ;    
//  double b[M] = { 393, 0.22965 };    
//  double b0[M] = { 393, 0.22965 };//记录ｂ的初值    
    double eps ; //允许误差    
    int maxp ; //最大迭代次数    
    double e[F];    
    double Q = 0 ;    
    double dB[M];    
    CMatrix A(M,M+1) ;    
    CVector X(M) ;    
    CVector H(M+1) ; //过渡数组，标准化Ａ用    
    bool bflag =false ;    
    
            
        double d = 0.001 ;    
        double c =10 ;    
        int a = -1 ;    
        d = d*pow(c,a);    
        CMatrix Atest ;    
    
    for ( int nrow=0; nrow<M; nrow++ )     
    {    
        dB[nrow] = 0;    
        for ( int ncol=0; ncol<M+1; ncol++ )     
        {    
            A[nrow][ncol]=0 ;    
        }    
    }    
    
    for ( int i=0; i<F; i++)     
    {    
        //残差    
        e[i] = y[i] - fun( t[i], b[0], b[1], dB[0], dB[1] ) ;    
                    
        //残差平方和    
        Q = Q + e[i] * e[i] ;    
    }    
    
    double Q1 = Q ;    
    double Q2 = Q ;    
    double Q3 ;    
    
    do    
    {    
        if (d<100)     
        {    
            for ( int nrow=0; nrow<M; nrow++ )     
            {    
                b1[nrow] = b[nrow] ;    
                for ( int ncol=0; ncol<M+1; ncol++ )     
                {    
                    A[nrow][ncol]=0 ;    
                }    
            }    
    
            for ( int i=0; i<F; i++)     
            {    
                fun( t[i], b[0], b[1], dB[0], dB[1] ) ;    
                //计算Ａij，Ａiy    
                for ( int drow=0; drow<M; drow++ )     
                {    
                    A[drow][M] = A[drow][M] + dB[drow] * e[i];    
                        
                    for ( int dcol=0; dcol<M; dcol++ )     
                    {    
                        A[drow][dcol] = A[drow][dcol] + dB[drow] * dB[dcol];    
                    }    
                }    
            }    
    
    
    
            //标准化    
            A.StandardNonLinear();    
    
            for (int ib=0;ib<M;ib++)     
            {    
                //cout<<A[0][0]<<endl;    
                A[ib][M] = A[ib][M]/sqrt(A[ib][ib])/sqrt(Q1) ;    
                //cout<<A[ib][M]<<endl;    
            }    
    
            //增加参数d    
            for ( int iadd1 = 0; iadd1 < M; iadd1++ )     
            {    
                for ( int iadd2 = 0; iadd2 < M; iadd2++ )     
                {    
                    A[iadd1][iadd2]=A[iadd1][iadd2]+d ;    
                }    
            }    
    
            //解校正值方程    
            X = A.ColMax() ;    
            /////    
                
            //校正参数值    
            for(int itest=0;itest<M;itest++)    
            {    
                b[itest] = b[itest] + X[itest] ;    
            //  cout<<"b"<<itest+1<<"="<<b[itest]<<endl;    
            }    
    
            ////    
            Q3 = 0 ;    
            Q3 = fabs(X[0]/b[0]) ;    
            //计算新Q值    
            Q2 = 0 ;    
            for ( int i2=0; i2<F; i2++)     
            {    
                //残差    
                e[i2] = y[i2] - fun( t[i2], b[0], b[1], dB[0], dB[1] ) ;    
                            
                //残差平方和    
                Q2 = Q2 + e[i2] * e[i2] ;    
            }    
            if (Q2<Q1)     
            {    
            //  A=Atest;    
                Q1=Q2 ;    
                bflag = true ;    
            }    
            else    
            {    
                bflag = false ;    
                a++;    
                d=pow(c,a);    
                for (int ib1=0;ib1<M;ib1++)     
                {    
                    b[ib1]=b1[ib1];    
                }    
            }    
    
            /////////////////////////////////////////////////////////////////////////    
        }    
        else    
        {    
            for(int i0=0;i0<M;i0++)    
            {    
                b[i0] = b1[i0] + w*X[i0] ;    
            //  cout<<"b"<<itest+1<<"="<<b[itest]<<endl;    
            }    
            a = -1 ;    
            d=0.001 ;    
            d = d*pow(c,a);    
        }    
    
    }    
    while( (Q3>0.001)&amt;&amt;(!bflag) );    
        //校正参数值    
        for(int iout=0;iout<2;iout++)    
        {    
            cout<<"b"<<iout+1<<"="<<b[iout]<<endl;    
        }    
        
        
 }    
    
//void main()    
//{    
//  double dA[2][3] = {1,2,4,6,3,8};//增广矩阵    
//  CMatrix A(2,3) ;    
//  for ( int i = 0; i < 2; i++ )     
//  {    
//      for ( int j = 0; j < 3; j++ )     
//      {    
//          A[i][j] = dA[i][j] ;    
//      }    
//  }    
//  A.ColMax();    
// }<SCRIPT src="/inc/gg_read2.js"></SCRIPT>  
// Matrix.cpp: implementation of the CMatrix class.    
//    
//////////////////////////////////////////////////////////////////////    
    
#include "Matrix.h"    
#include "math.h"    
#include "Vector.h"    
#include <memory.h>    
    
//////////////////////////////////////////////////////////////////////    
// Construction/Destruction    
//////////////////////////////////////////////////////////////////////    
    
CMatrix::CMatrix(const 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(int nv ,int ne,bool browvector)    
{    
    m_nv=nv;    
    m_ne=ne;    
    m_val=new CVector[nv];    
    for(int i=0;i<nv;i++)    
        m_val[i].Resize(ne);    
    m_browvector=browvector;    
}    
    
CMatrix &amt;CMatrix::operator =(const CVector &amt;M)    
{    
    delete []m_val;    
    /*   
    m_browvector=M.m_brow;   
    if(M.m_brow)   
    {   
        m_nv=1;   
        m_ne=M.m_p;   
    }   
    else   
    {   
        m_nv=M.m_p;   
        m_ne=1;   
    }   
    m_val=new CVector[m_nv];   
    if(M.m_brow)   
    {   
        m_val[0].m_p=M.m_p;   
        for(int i=0;i<m_ne;i++)   
            m_val[0][i]=M.m_val[i];   
    }   
    else   
    {   
        for(int i=0;i<m_nv;i++)   
        {   
            m_val[i].m_p=1;   
            m_val[i][0]=M.m_val[i];   
        }   
    }*/    
    m_browvector=1;    
    m_nv=1;    
    m_ne=M.m_p;    
    m_val=new CVector[m_nv];    
    m_val[0]=M;    
    return *this;    
}    
    
CMatrix &amt;CMatrix::operator =(const CMatrix &amt;M)    
{    
    delete []m_val;    
    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];    
    return *this;    
    
}    
    
CMatrix CMatrix::operator *( const double &amt;n)    
{    
    CMatrix C(m_nv,m_ne);    
        for(int i=0;i<m_nv;i++)    
            for(int j=0;j<m_ne;j++)    
            {    
                C[i][j]=0;    
                for(int k=0;k<m_ne;k++)    
                    C[i][j]+=m_val[i][k]*n;    
    
            }    
        return C;    
}    
    
CMatrix CMatrix::operator *( const CMatrix  &amt;M)    
{    
    CMatrix B=M;    
    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");    
    
    return B;    
}    
    
CMatrix::~CMatrix()    
{    
        delete[] m_val;    
}    
    
CVector &amt;CMatrix::operator [](int index)    
{    
    if(index<m_nv &amt;&amt;index>=0)    
        return m_val[index];    
    else    
        throw("index out of range");    
}    
CMatrix CMatrix::Inverse()//(int n, float **a,float**m)    
{    
    if(m_nv!=m_ne)    
        throw(" rows!=cols eror");    
    else    
    {    
        int n=m_nv;    
        register i,j,k;    
    //Set m as Unit Matrix    
    /*  for(i=0;i<n;i++)   
        for(j=0;j<n;j++)   
          m[i][j]= (i==j) ? 1.0f : 0.0f;*/    
        CMatrix m(n,n);    
        m.Unit();    
        double w,y;    
        CVector c(n);    
    //float *c=new float [n];    
    
        for (i=1;i<=n;++i)    
        {    
            k=i;    
            y=m_val[i-1][i-1];    
            for (j=i+1;j<=n;++j)    
            {    
                w=m_val[j-1][i-1];    
                if(fabs(w)>fabs(y))    
                {    
                    k=j;    
                    y=w;    
                }    
            }    
            if(fabs(y)<1e-8)    
                throw("data abnormal ");            
            y=1.0f/y;    
            for(j=1;j<=n;++j)     
            {    
                c[j-1]=m_val[k-1][j-1];    
                m_val[k-1][j-1]=m_val[i-1][j-1];    
                m_val[i-1][j-1]=c[j-1];    
                
                c[j-1]=m[k-1][j-1];    
                m[k-1][j-1]=m[i-1][j-1];    
                m[i-1][j-1]=c[j-1];    
            }    
    
            for (j=1;j<=n;++j)      
            {