|
Hello,
I had the same problem and did not find anything in matlab.
Below is the code, which I wrote up for this purpose. It
contains the reference to the source, which I used. Perhaps
I will consider uploading it to file exchange, but for now
I will just copy-paste it. Hope you will cope with text
wrapping and figure it out.
Yuri
------------------------
%*******************************
function [sesl, cisl, sein, ciin] = secisi(X,Y,pr);
% Returns Standard Errors and Confidence Intervals
% for the Slope and Intercept of a regression.
% Author: Yuri Geshelin, 2002.
% Source: N. R. Draper and H. Smith, 1966:
% Applied Regression Analysis, 407 p. Publ. by
% John Wiley and Sons, Inc. New York - London - Sydney.
% Input:
% 1-D arrays X and Y, for which the regression
% Y = B0 + B1 * X is computed
% pr - probability (e.g. if we want 95%
% confidence level, then pr = 95.
% Output:
% sesl - standard error for slope B1
% clsl - pr%-confindence interval for slope B1
% sein - standard error for intercept B0
% clin - pr%-confindence interval for intercept B0
if pr < 10 | pr > 99.9
error('pr should be between 10 and 99.9!')
end
n = length(X);
if length(Y) ~= n
error('X, Y should be the same lengths!')
end
n2 = n - 2; % n - 2 degrees of freedom
if n2 < 121 % if n2 = 120, then we already would
% have to use the next-to-last line
% of the table on p. 305: this code is only
% for the last line.
error('This code only works with n > 120!')
end
% The following line is the result of
% prob = 100 * ...
% (1 - [0.9:-0.2:0.3 0.2 0.1 0.05 0.02 0.01 0.001]) ...
% (see p. 305):
prob = [10 30 50 70 80 90 95 98 99 99.9]; % percentage
pp = [0.126 0.385 0.674 1.036 1.282 ...
1.645 1.960 2.326 2.576 3.291]; % last line of the
table,
% p. 305.
% Compute total corrected sum of squares (SS) of the Y's
% (see p.14 and table on p. 16):
SStot = sum((Y - mean(Y)).^2);
% Compute SS due to regression (see p. 16):
SX2 = sum(X .^ 2);
SSrgr = ...
(sum(X .* Y) - sum(X) * sum(Y) / n) ^2 / ...
(SX2 - sum(X) ^2 / n);
% Compute SS due to residual (see p. 16):
SSres = SStot - SSrgr;
% Compute mean square about regression
% (see tables on pp. 15, 16):
s2 = SSres / n2;
% Compute estimate of variance of slope B1 (p. 19)
sesl = sqrt(s2 / sum((X - mean(X)) .^2 ));
% Found the interpolated value from the table
% (page 305):
val_int = interp1(prob,pp,pr);
% Compute confidence limit for slope:
cisl = val_int * sesl;
% Compute estimate of variance of intercept B0 (p. 21)
sein = sesl * sqrt(SX2 / n);
% Compute confidence limit for intercept:
ciin = val_int * sein;
%***********************************
|
Feed for this Thread