"Liqing " <liqingj@gmail.com> wrote in message
<g7atfm$a9$1@fred.mathworks.com>...
> Hi,
>
> I have two columns of data: x and y. I fitted y against x
> with exponential regression. The resulting equation is:
> y = a*exp(b*x).
>
> My question is how can I estimate the standard deviation
> of a and b?
>
> Your help is greatly appreciated!
I describe how to do this in my optimization
tip and tricks file. Find it here:
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
objectId=8553&objectType=FILE
Make up some data as an example...
x = rand(100,1)*42;
y = 1.5*exp(.5*x) + randn(size(x))/10;
plot(x,y,'o')
% Now, fit it using the model y=a*exp(b*x)
abfun = @(ab) ab(1)*exp(ab(2)*x)y;
abstart = [1 1]';
abEst = fminsearch(@(ab) sum(abfun(ab).^2),abstart)
% abEst =
% 1.493
% 0.50739
% (reasonably close)
% Compute the Jacobian matrix, evaluated at
% the final estimate. (If you don't know how
% to do this yourself, download my derivest
% tools from the file exchange.) Thus, either
jac = [exp(abEst(2)*x),abEst(1)*exp(abEst(2)*x).*x];
% or ...
jac = jacobianest(abfun,abEst);
% Compute the residuals
res = abfun(abEst);
% Compute the standard error of the residuals.
% Note that I've divided by n2, since there are
% 2 parameters that we estimated.
ser = sqrt(sum(res.^2)/(length(x)2));
% Compute a covariance matrix of the parameters
% (yeah, I know I used inv.)
Cab = ser^2*inv(jac'*jac);
% The parameter standard deviations...
Sab = sqrt(diag(Cab))
Sab =
0.011591
0.0060296
% +/ 2*sigma bounds on ab
[abEst2*Sab,abEst,abEst+2*Sab]
ans =
1.4699 1.493 1.5162
0.49533 0.50739 0.51945
You can find derivest here:
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
objectId=13490&objectType=FILE#
HTH,
John
