leasqr problem in the partial derivative

5 messages
Open this post in threaded view
|

leasqr problem in the partial derivative

 I went to the examples of leasqr and copied example #2, then adjusted it for my data.% Define functions leasqrfunc = @(x, p) p(1) * exp (-p(2) * x ) +p(3); leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)*x), -p(1)*x.*exp(-p(2)*x),0*x+1];what I did was added +p(3) to the original function. I then calculated the partial diff with respect to p(3). this is " 1 ".when I put a 1 in the 3rd position in leasqrdfdp, then the program does not work. but when I make the 3rd term to by 0*x + 1  then the program runs and gives the correct answer.should just +1 work???here is the program that works % Define functions leasqrfunc = @(x, p) p(1) * exp (-p(2) * x ) +p(3); leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)*x), -p(1)*x.*exp(-p(2)*x),0*x+1]; % generate test data t = [0 2 4 6 13]'; #p = [1; 0.1]; #data = leasqrfunc (t, p); data=[1520   1457   1426   1408   1387]'; #rnd = [0.352509; -0.040607; -1.867061; -1.561283; 1.473191; ... #       0.580767;  0.841805;  1.632203; -0.179254; 0.345208]; % add noise % wt1 = 1 /sqrt of variances of data % 1 / wt1 = sqrt of var = standard deviation wt1 = (1 + 0.01 * t) ./ sqrt (data);  #data = data + 0.05 * rnd ./ wt1; #wt1=[1 1 1 1 1]' *1; % Note by Thomas Walter : % % Using a step size of 1 to calculate the derivative is WRONG !!!! % See numerical mathbooks why. % A derivative calculated from central differences need: s  %     step = 0.001...1.0e-8 % And onesided derivative needs: %     step = 1.0e-5...1.0e-8 and may be still wrong F = leasqrfunc; dFdp = leasqrdfdp; % exact derivative % dFdp = dfdp;     % estimated derivative dp = [01; 0.1;1]; pin = [160; 1.0; 1150];  stol=0.000001; niter=50; minstep = [0.01; 0.00001;.1]; maxstep = [010; 0.2; 1000]; options = [minstep, maxstep]; global verbose; verbose = 1; [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...    leasqr (t, data, pin, F, stol, niter, wt1, [], dFdp, options);        f1    p1    xx=0:.1:30; yy=F(xx,p1); hold on plot(xx,yy)    hold off    -- DAS _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|

Fwd: leasqr problem in the partial derivative

 added Olaf Till---------- Forwarded message ----------From: Doug Stewart Date: Thu, Sep 21, 2017 at 8:08 PMSubject: leasqr problem in the partial derivativeTo: Help GNU Octave <[hidden email]>I went to the examples of leasqr and copied example #2, then adjusted it for my data.% Define functions leasqrfunc = @(x, p) p(1) * exp (-p(2) * x ) +p(3); leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)*x), -p(1)*x.*exp(-p(2)*x),0*x+1];what I did was added +p(3) to the original function. I then calculated the partial diff with respect to p(3). this is " 1 ".when I put a 1 in the 3rd position in leasqrdfdp, then the program does not work. but when I make the 3rd term to by 0*x + 1  then the program runs and gives the correct answer.should just +1 work???here is the program that works % Define functions leasqrfunc = @(x, p) p(1) * exp (-p(2) * x ) +p(3); leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)*x), -p(1)*x.*exp(-p(2)*x),0*x+1]; % generate test data t = [0 2 4 6 13]'; #p = [1; 0.1]; #data = leasqrfunc (t, p); data=[1520   1457   1426   1408   1387]'; #rnd = [0.352509; -0.040607; -1.867061; -1.561283; 1.473191; ... #       0.580767;  0.841805;  1.632203; -0.179254; 0.345208]; % add noise % wt1 = 1 /sqrt of variances of data % 1 / wt1 = sqrt of var = standard deviation wt1 = (1 + 0.01 * t) ./ sqrt (data);  #data = data + 0.05 * rnd ./ wt1; #wt1=[1 1 1 1 1]' *1; % Note by Thomas Walter : % % Using a step size of 1 to calculate the derivative is WRONG !!!! % See numerical mathbooks why. % A derivative calculated from central differences need: s  %     step = 0.001...1.0e-8 % And onesided derivative needs: %     step = 1.0e-5...1.0e-8 and may be still wrong F = leasqrfunc; dFdp = leasqrdfdp; % exact derivative % dFdp = dfdp;     % estimated derivative dp = [01; 0.1;1]; pin = [160; 1.0; 1150];  stol=0.000001; niter=50; minstep = [0.01; 0.00001;.1]; maxstep = [010; 0.2; 1000]; options = [minstep, maxstep]; global verbose; verbose = 1; [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...    leasqr (t, data, pin, F, stol, niter, wt1, [], dFdp, options);        f1    p1    xx=0:.1:30; yy=F(xx,p1); hold on plot(xx,yy)    hold off    -- DAS -- DAS _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|

Re: Fwd: leasqr problem in the partial derivative

 On Thu, Sep 21, 2017 at 10:12:31PM -0400, Doug Stewart wrote: > added Olaf Till > ---------- Forwarded message ---------- > From: Doug Stewart <[hidden email]> > Date: Thu, Sep 21, 2017 at 8:08 PM > Subject: leasqr problem in the partial derivative > To: Help GNU Octave <[hidden email]> > > > I went to the examples of leasqr and copied example #2, then adjusted it > for my data. > > % Define functions >  leasqrfunc = @(x, p) p(1) * exp (-p(2) * x ) +p(3); >  leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)*x), > -p(1)*x.*exp(-p(2)*x),0*x+1]; > > what I did was added +p(3) to the original function. I then calculated the > partial diff with respect to p(3). this is " 1 ". > when I put a 1 in the 3rd position in leasqrdfdp, then the program does not > work. but when I make the 3rd term to by 0*x + 1  then the program runs and > gives the correct answer. > > should just +1 work??? 'x' is a vector, no scalar, and 'leasqrdfdp' returns a (numel(x) by 3) matrix, not a vector... So '1' is wrong (error: horizontal dimensions mismatch (10x2 vs 1x1)), 'ones(numel(x), 1)' is ok (no spaces between 'ones' or 'numel' and '('), which is the same as '0*x+1'. Off topic: I recommend to use nonlin_curvefit instead of leasqr, the default algorithm is the same. Olaf -- public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave signature.asc (836 bytes) Download Attachment