Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



Re: standard errors of fminsearch parameters
Posted:
Dec 11, 2013 10:31 AM


On 12/10/2013 5:54 PM, Derek wrote: > "John D'Errico" <woodchips@rochester.rr.com> wrote in message > <fele3u$jf2$1@fred.mathworks.com>... >> "Cecilia Persson" <cecper@passagen.se> wrote in message <fel5qr$716 >> $1@fred.mathworks.com>... >> > does someone know how I can get the standard errors of the > >> parameters found by fminsearch? (i.e. not the errors of > each data >> value compared to the calculated one but the > parameters fitted to >> my nonlinear data) >> >> n = 200; >> % n data points, in the interval [1,1]. >> x = rand(n,1)*21; >> >> % true parameters >> ab0 = [2 .5]; >> >> % a function that predicts y, given x and the parameters >> predfun = @(ab) ab(1)*exp(ab(2)*x); >> >> % our data, to be used for estimation >> y = randn(size(x))*.1 + predfun(ab0); >> >> plot(x,y,'o') >> >> % the estimation step, using fminsearch >> sumofsquares = @(ab) sum((y  predfun(ab)).^2); >> abstart = [1 1]; >> abfinal = fminsearch(sumofsquares,abstart) >> >> % degrees of freedom in the problem >> dof = n  2; >> >> % standard deviation of the residuals >> sdr = sqrt(sum((y  predfun(abfinal)).^2)/dof) >> >> % jacobian matrix >> J = jacobianest(predfun,abfinal) >> >> % I'll be lazy here, and use inv. Please, no flames, >> % if you want a better approach, look in my tips and >> % tricks doc. >> Sigma = sdr^2*inv(J'*J); >> >> % Parameter standard errors >> se = sqrt(diag(Sigma))' >> >> % which suggest rough confidence intervalues around >> % the parameters might be... >> abupper = abfinal + 2*se >> ablower = abfinal  2*se >> >> >> I'll let you decide what the confidence level is, or if >> you would prefer to use a tstatistic for the above >> computation. >> >> The output from the code above is... >> >> >> abfinal = >> 1.98856241997929 0.503604377341979 >> >> sdr = >> 0.102056089110678 >> >> se = >> 0.00786391742840254 0.0063796968802395 >> >> abupper = >> 2.00429025483609 0.516363771102458 >> >> ablower = >> 1.97283458512248 0.4908449835815 >> >> >> You can find my optimization tips and tricks document at: >> >> http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do? >> objectId=8553&objectType=FILE >> >> (Beware, the newsreader tends to wrap these URLs. Paste >> in BOTH lines!) >> >> You can find jacobianest here: >> >> http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do? >> objectId=13490&objectType=FILE# >> >> CAVEATS: >> The computations above assume many things about >> your problem. These assumptions include a normal >> (Gaussian) iid error structure, that the optimizaton >> has converged to a minimizer of the error surface, >> no lack of fit in the model, and probably a few other >> things I've forgotten to list here. >> >> HTH, >> John D'Errico > > Dear John, > This example has worked very well for me but now I am trying to > reference this method. Where does Sigma = sdr^2*inv(J'*J) come from? > Any help/references would be appreciated. Thanks! Derek You can try either of these, as well as many other textbooks: Examples That Use Standard Algorithms :: Tutorial (Optimization Toolbox?)
Draper, Norman R., and Smith, Harry. /Applied Regression Analysis/, Third Edition. New York: WileyInterscience, 1998.
Seber, G.A.F., and Wild, C.J. /Nonlinear Regression/, New York: WileyInterscience, 2003.
Alan Weiss MATLAB mathematical toolbox documentation



