Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: standard errors of fminsearch parameters
Replies: 6   Last Post: May 29, 2014 9:59 AM

 Messages: [ Previous | Next ]
 Alan Weiss Posts: 1,430 Registered: 11/27/08
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)*2-1;
>>
>> % 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 t-statistic 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:
>>
>> objectId=8553&objectType=FILE
>>
>> (Beware, the newsreader tends to wrap these URLs. Paste
>> in BOTH lines!)
>>
>> You can find jacobianest here:
>>
>> 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: Wiley-Interscience, 1998.

Seber, G.A.F., and Wild, C.J. /Nonlinear Regression/, New York:
Wiley-Interscience, 2003.

Alan Weiss
MATLAB mathematical toolbox documentation

Date Subject Author