syms a b c m y x real;
f=a*exp(-b*x)+c*exp(-m*x);
Jsym=jacobian(f,[a b c m]);
data_1=[30 45 75 105 150 200 300 375 475 600 675 750 900];
obs_1=[4.8 4.3 3.45 2.95 2.67 2.6 2.67 2.81 3.16 3.39 3.5 4 4.25];

a0=2.5699; b0=0.0176;c0=2.7675;m0=-0.0003776;
y_init = a0*exp(-b0*data_1)+c0*exp(-m0*data_1);
Ndata=length(obs_1);
Nparams=4;
n_iters=50;
lamda=0.01;

updateJ=1;
a_est=a0;
b_est=b0;
c_est=c0;
m_est=m0;

for it=1:n_iters
    if updateJ==1
     
        J=zeros(Ndata,Nparams);
        for i=1:length(data_1)
            J(i,:)=[exp(-b_est*data_1(i)) -a_est*data_1(i)*exp(-b_est*data_1(i))+exp(-m_est*data_1(i))-c_est*data_1(i)*exp(-m_est*data_1(i))];
        end
        y_est = a_est*exp(-b_est*data_1)+c_est*exp(-m_est*data_1);
        d=obs_1-y_est;
        H=J'*J;
        if it==1
            e=dot(d,d);
        end
    end

    H_lm=H+(lamda*eye(Nparams,Nparams));
    dp=pinv(J'*d(:))*(H_lm);
    g = J'*d(:);
    a_lm=a_est+dp(1);
    b_lm=b_est+dp(2);
    c_lm=c_est+dp(3);
    m_lm=m_est+dp(4);
    y_est_lm = a_lm*exp(-b_lm*data_1)+c_lm*exp(-m_lm*data_1);
    d_lm=obs_1-y_est_lm;
    e_lm=dot(d_lm,d_lm);
    if e_lm<e
        lamda=lamda/10;
        a_est=a_lm;
        b_est=b_lm;
        c_est=c_lm;
        m_est=m_lm;
        e=e_lm;
        disp(e);
        updateJ=1;
    else
        updateJ=0;
        lamda=lamda*10;
    end
end
a_est;
b_est;
c_est;
m_est;

3Â¥: Originally posted by pdl9527 at 2017-05-14 21:39:39
J(i,=;

ÕâÀ︳ֵÓÐÎÊÌ⣬Ä㼺¾­¶¨ÒåÁËJµÄά¶È

J=zeros(Ndata,Nparams);


ÕâÀïºóÃæ¼ÆËãµÄÖµºÍÄ㶨ÒåµÄά¶È²»Æ¥Å䣬Çë¼ì²é¼ÆËãµÄʱºòÊÇ·ñ´æÔÚÎÊÌ⣬±ÈÈç¾ØÕóÔËËãÊÇ·ñÓõã³ËµÈµÈ¡£...

J(i,:)=[exp(-b_est*data_1(i)) -a_est*data_1(i)*exp(-b_est*data_1(i))+exp(-m_est*data_1(i))-c_est*data_1(i)*exp(-m_est*data_1(i))];

J=zeros(Ndata,Nparams);