Estimation Problem in Newton Method
in [Matlab]
|
From: Antony zerrar on 14 Jun 2010 06:22 Hi guys, I got a problem in fitting a distribution in more than 5 observations by using newton method. Here is the code. Thanks! d = [-1 5 2 -2 4 3 -3 6 3 1 -2]; % Data x0=[0 1]'; % Starting values x(:,1)=x0; N=1; for n=N:15 m = sum(d-x(1,n))/x(2,n); s= -3/(2*(x(2,n)))+(1/(2*(x(2,n)^2)))*(sum((d-x(1,n)).^2)); G = [m ; s]; % Gradient Vector mm = (-3/x(2,n)); ms = (-sum(d-x(1,n))/(2*x(2,n)^2)); ss = (3/(2*(x(2,n)^2)))-(sum((d-x(1,n)).^2))/(x(2,n)^3); H = [ mm ms ; ms ss]; % Hessian Vector x(:,n+1) = x(:,n) - inv(H)*G; end
From: Roger Stafford on 14 Jun 2010 19:59 "Antony zerrar" <zerrar(a)yahoo.com> wrote in message <hv4vse$9h1$1(a)fred.mathworks.com>... > Hi guys, > > I got a problem in fitting a distribution in more than 5 observations by using newton method. Here is the code. Thanks! > > d = [-1 5 2 -2 4 3 -3 6 3 1 -2]; % Data > x0=[0 1]'; % Starting values > x(:,1)=x0; > N=1; > > for n=N:15 > m = sum(d-x(1,n))/x(2,n); > s= -3/(2*(x(2,n)))+(1/(2*(x(2,n)^2)))*(sum((d-x(1,n)).^2)); > G = [m ; s]; % Gradient Vector > mm = (-3/x(2,n)); > ms = (-sum(d-x(1,n))/(2*x(2,n)^2)); > ss = (3/(2*(x(2,n)^2)))-(sum((d-x(1,n)).^2))/(x(2,n)^3); > H = [ mm ms ; ms ss]; % Hessian Vector > x(:,n+1) = x(:,n) - inv(H)*G; > end Apparently the function you are trying to minimize or maximize is: f(x1,x2) = -3/2*log(x2) - sum((d-x1)^2)/(2*x2) If so, I see what looks like a couple of errors in your second partial derivatives. You should have mm = -length(d)/x2 = -11/x2 instead of -3/x2 and ms = -sum(d-x1)/x2^2 instead of -sum(d-x1)/(2*x2^2) Roger Stafford
From: Antony zerrar on 15 Jun 2010 11:25 Thanks a lot Roger! Silly mistake there
From: Antony zerrar on 17 Jun 2010 06:22 "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hv6foa$3ho$1(a)fred.mathworks.com>... > "Antony zerrar" <zerrar(a)yahoo.com> wrote in message <hv4vse$9h1$1(a)fred.mathworks.com>... > > Hi guys, > > > > I got a problem in fitting a distribution in more than 5 observations by using newton method. Here is the code. Thanks! > > > > d = [-1 5 2 -2 4 3 -3 6 3 1 -2]; % Data > > x0=[0 1]'; % Starting values > > x(:,1)=x0; > > N=1; > > > > for n=N:15 > > m = sum(d-x(1,n))/x(2,n); > > s= -3/(2*(x(2,n)))+(1/(2*(x(2,n)^2)))*(sum((d-x(1,n)).^2)); > > G = [m ; s]; % Gradient Vector > > mm = (-3/x(2,n)); > > ms = (-sum(d-x(1,n))/(2*x(2,n)^2)); > > ss = (3/(2*(x(2,n)^2)))-(sum((d-x(1,n)).^2))/(x(2,n)^3); > > H = [ mm ms ; ms ss]; % Hessian Vector > > x(:,n+1) = x(:,n) - inv(H)*G; > > end > > Apparently the function you are trying to minimize or maximize is: > > f(x1,x2) = -3/2*log(x2) - sum((d-x1)^2)/(2*x2) > > If so, I see what looks like a couple of errors in your second partial derivatives. You should have > > mm = -length(d)/x2 = -11/x2 instead of -3/x2 > > and > > ms = -sum(d-x1)/x2^2 instead of -sum(d-x1)/(2*x2^2) > > Roger Stafford Hi Roger, Can I ask you how would it be possible to change this model to fsolve or fmincon? Thanks!
From: Roger Stafford on 17 Jun 2010 13:30
"Antony zerrar" <zerrar(a)yahoo.com> wrote in message <hvct0c$n90$1(a)fred.mathworks.com>... > Can I ask you how would it be possible to change this model to fsolve or fmincon? I don't see why you are using matlab for this problem at all, Antony. Using elementary calculus you can easily determine that in the upper half of the x1,x2 plane (x2 > 0) your function, f(x1,x2), attains a local maximum at: x1 = 1/11*sum(d) = mean(d) x2 = 1/3*sum((d-mean(d)).^2) = 10/3*var(d) That is where both components of the gradient, m and s, become zero. A 'surf' plot in the neighborhood of this point should make this evident. In the lower half of the x1,x2 plane (x2 < 0) there is no local maximum or minimum as the function varies from minus infinity to plus infinity throughout that half-plane. Newton's method (properly done,) fsolve, and fmincon should all find that solution easily if given an appropriate initial estimate, but why do it that way when the solution is so easily found by hand? I am assuming here, based on your earlier statement, that your function is: f(x1,x2) = -3/2*log(abs(x2)) - sum((d-x1)^2)/(2*x2) where d = [-1 5 2 -2 4 3 -3 6 3 1 -2]. Roger Stafford |