Méthode de Newton Raphson avec matlab

laribi · April 2019

Bonjour, je voudrais estimer les paramètres a, b, c et d par la méthode de Newton Raphson sous matlab, donc j'ai calculé la hessienne et le grad de la fonction vraisemblance et j'ai l'écrit sous matlab comme suit au dessous, quand je l’exécute, les SK1,SK2,SK3,SK21,... sont tous des vecteurs infinis
par exemple

SK2 =
   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf   Inf

je ne comprends pas pourquoi ils sont tous des vecteurs infinis

syms a b c d;
a0=0.8;
b0=2;
d0=0.5;
c0=1.2
N=15;
x=linspace(0.001,20,N)
K1=[];
K2=[];
SK2=[];
SK3=[];
fL=[];
G=[];
SK21=[];
SK22=[];
SK23=[];
i=1:N
G(i)= (c0-1)/c0.*gammainc(a0.*x(i),b0);
fL(i) = (1/c0).*(a0^(b0+1).*x(i).^(b0-1).*(b0+d0.*x(i)).*exp(-a0.*x(i)))/((d0+a0).*gamma(b0+1));
K1(i)=G(i)/(G(i)+fL(i));
K2(i)=1-K1(i)
SK3(i)=K2(i)/(b+d.*x(i)).^2;
SK2(i)=K2(i)/(b+d.*x(i));
SK21(i)=(K2(i).*x(i))/(b+d.*x(i))
SK22(i)=(K2(i).*x(i))/(b+x(i)*d).^2;
SK23(i)=(K2(i).*x(i).^2)/(b+d.*x(i)).^2;
S=sum(x(i));
S1=sum(log(x(i)));
k1=sum (K1(i));
k2=sum (K2(i));
sk3=sum (SK3(i));
sk2=sum(SK2(i));
sk21=sum(SK21(i));
sk22=sum(SK22(i));
sk23=sum(SK23(i));
K=a0+1;
F1=a^2*S+a*(k2-N*(b+1)+d*S)-N*(b+1)*d;
F2=N.*log(a)+S1-N.*psi(b+1)+sk2;
F3=sk21-(N/(a+d));
F=size(3,1);
F=[F1 ; F2 ; F3 ];
H=size(3,3);
H11=k2/(d+a)^2-(b+1)/a^2;
H12= N/a;
H13=N/(d+a)^2;
H21=-N*(a+d);
H22=-N*psi(1,b+1)-sk3;
H23=-sk22;
H31=a.*S-N*(b+1);
H32=-sk22;
H33=-sk23+1/(a+d)^2;

marsup · April 2019

Qu'est ce que tu veux dire avec l'instruction suivante ?

i = 1:N
G(i) = ...

Tu ne voudrais pas parler d'une boucle, par hasard ?

for i = 1:N
  G(i) = ...
  ...
  ...
end

marsup · April 2019

Ah ok.
Je pense que le problème est simplement que tu as oublié les `./`
Au lieu de

SK3(i)=K2(i)/(b+d.*x(i)).^2;

Tu veux dire

SK3(i) = K2(i) ./ (b+d.*x(i)).^2;

laribi · April 2019

merci pour votre réponse, j'ai essayé de faire ça et j'ai obtenu le même résultat, ça ne marche pas

laribi · April 2019

Inutile d'ouvrir une nouvelle discussion pour ton problème. AD

Bonjour, je veux estimer les paramètres a b et d par la méthode de Newton Raphson sous matlab, voila mon code, après l’exécution j'ai trouvé que a et b et d sont des nombres complexes je ne comprends pas pourquoi, toutes mes donnes initiales sont des réels

syms a b d ;
a0=0.8;
b0=0.1;
d0=0.5;
c0=1.2;
N=15;
K1=[];
K2=[];
SK2=[];
SK3=[];
fL=[];
G=[];
SK21=[];
SK22=[];
SK23=[];
n=10
aold=a0;
bold=b0;
dold=d0;
for k = 0:n
x=linspace(0.001,20,N)
i=1:N
G(i)= (c0-1)/c0.*gammainc(aold.*x(i),bold);
fL(i) = (1/c0).*(aold^(bold+1).*x(i).^(bold-1).*(bold+dold.*x(i)).*exp(-aold.*x(i)))/((dold+aold).*gamma(bold+1));
K1(i)=G(i)/(G(i)+fL(i))
K2(i)=1-K1(i)
SK3(i)=K2(i)/(bold+dold.*x(i)).^2
SK2(i)=K2(i)/(bold+dold.*x(i))
SK21(i)=(K2(i).*x(i))/(bold+dold.*x(i))
SK22(i)=(K2(i).*x(i))/(bold+x(i)*dold).^2
SK23(i)=(K2(i).*x(i).^2)/(bold+dold.*x(i)).^2
S=sum(x(i))
S1=sum(log(x(i)))
k1=sum (K1(i))
k2=sum (K2(i))
sk3=sum (SK3(i))
sk2=sum(SK2(i))
sk21=sum(SK21(i))
sk22=sum(SK22(i))
sk23=sum(SK23(i))
F1=aold^2*S+aold*(k2-N*(bold+1)+dold*S)-N*(bold+1)*dold
F2=N.*log(aold)+S1-N.*psi(bold+1)+sk2
F3=sk21-(N/(aold+dold))
F=size(3,1)
F=[F1 ; F2 ; F3 ]
H=size(3,3)
H11=k2/(dold+aold)^2-(bold+1)/aold^2
H12= N/aold
H13=N/(dold+aold)^2
H21=-N*(aold+dold)
H22=-N*psi(1,bold+1)-sk3
H23=-sk22
H31=aold.*S-N*(bold+1)
H32=-sk22
H33=-sk23+1/(aold+dold)^2
H=[H11,H12,H13;H21,H22,H23;H31,H32,H33]

    A=inv(H)*F;
    a=aold-A(1,1)
    b=bold-A(2,1)
    d=dold-A(3,1)
    aold=a;
    bold=b;
    dold=d;
    
    if norm(F)<0.0001
        break;
    end
end

marsup · April 2019

C'est peut-être que $F_2 = N \cdot \log(\dots) \dots$ reçoit un argument $a_{old}$ négatif à un moment donné ?

laribi · April 2019

ahhhh oui ta raison, avez-vous une idée pour choisir les donnes initiales pour obtenir a, b et d des réelles positifs?

marsup · April 2019

Aucune idée, d'ailleurs je ne reconnais même pas la méthode de Newton dans tout ce que tu as écrit.

À quel problème t'intéresses-tu au juste ?

laribi · April 2019

je ne peux pas te dire parce que c'est mon sujet de thèse, merci beaucoup pour ta réponse :-)

Méthode de Newton Raphson avec matlab

Réponses

Bonjour!

Catégories

In this Discussion