"Yuri Geshelin" wrote in message <frbt1b$9nn$1@fred.mathworks.com>...
> 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 copypaste it. Hope you will cope with text
> wrapping and figure it out.
>
> Yuri
> 
I took the liberty to improve on this by including a ttable, allowing for any degrees of freedom:
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:
% 1D 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 <50  pr > 99.9
error('pr should be between 50 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
% 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);
val_int = tvalue(pr,n2,2);
% 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;
end
function [ tval ] = tvalue( confint, df, sides )
%tvalue returns interpolated values from a ttable
one = [75 80 85 90 95 97.5 99 99.5 99.75 99.9 99.95];
two = [50 60 70 80 90 95 98 99 99.5 99.8 99.9];
tabledf = [1:30 40 50 60 80 100 120 10^40];
t_table = [ 1 1.3760 1.9630 3.0780 6.3140 12.710 31.820 63.660 127.30 318.30 636.60;
0.81600 1.06100 1.3860 1.8860 2.9200 4.3030 6.9650 9.9250 14.090 22.330 31.600;
0.76500 0.97800 1.2500 1.6380 2.3530 3.1820 4.5410 5.8410 7.4530 10.210 12.920;
0.74100 0.94100 1.1900 1.5330 2.1320 2.7760 3.7470 4.6040 5.5980 7.1730 8.6100;
0.72700 0.92000 1.1560 1.4760 2.0150 2.5710 3.3650 4.0320 4.7730 5.8930 6.8690;
0.71800 0.90600 1.1340 1.4400 1.9430 2.4470 3.1430 3.7070 4.3170 5.2080 5.9590;
0.71100 0.89600 1.1190 1.4150 1.8950 2.3650 2.9980 3.4990 4.0290 4.7850 5.4080;
0.70600 0.88900 1.1080 1.3970 1.8600 2.3060 2.8960 3.3550 3.8330 4.5010 5.0410;
0.70300 0.88300 1.1000 1.3830 1.8330 2.2620 2.8210 3.2500 3.6900 4.2970 4.7810;
0.70000 0.87900 1.0930 1.3720 1.8120 2.2280 2.7640 3.1690 3.5810 4.1440 4.5870;
0.69700 0.87600 1.0880 1.3630 1.7960 2.2010 2.7180 3.1060 3.4970 4.0250 4.4370;
0.69500 0.87300 1.0830 1.3560 1.7820 2.1790 2.6810 3.0550 3.4280 3.9300 4.3180;
0.69400 0.87000 1.0790 1.3500 1.7710 2.1600 2.6500 3.0120 3.3720 3.8520 4.2210;
0.69200 0.86800 1.0760 1.3450 1.7610 2.1450 2.6240 2.9770 3.3260 3.7870 4.1400;
0.69100 0.86600 1.0740 1.3410 1.7530 2.1310 2.6020 2.9470 3.2860 3.7330 4.0730;
0.69000 0.86500 1.0710 1.3370 1.7460 2.1200 2.5830 2.9210 3.2520 3.6860 4.0150;
0.68900 0.86300 1.0690 1.3330 1.7400 2.1100 2.5670 2.8980 3.2220 3.6460 3.9650;
0.68800 0.86200 1.0670 1.3300 1.7340 2.1010 2.5520 2.8780 3.1970 3.6100 3.9220;
0.68800 0.86100 1.0660 1.3280 1.7290 2.0930 2.5390 2.8610 3.1740 3.5790 3.8830;
0.68700 0.86000 1.0640 1.3250 1.7250 2.0860 2.5280 2.8450 3.1530 3.5520 3.8500;
0.68600 0.85900 1.0630 1.3230 1.7210 2.0800 2.5180 2.8310 3.1350 3.5270 3.8190;
0.68600 0.85800 1.0610 1.3210 1.7170 2.0740 2.5080 2.8190 3.1190 3.5050 3.7920;
0.68500 0.85800 1.0600 1.3190 1.7140 2.0690 2.5000 2.8070 3.1040 3.4850 3.7670;
0.68500 0.85700 1.0590 1.3180 1.7110 2.0640 2.4920 2.7970 3.0910 3.4670 3.7450;
0.68400 0.85600 1.0580 1.3160 1.7080 2.0600 2.4850 2.7870 3.0780 3.4500 3.7250;
0.68400 0.85600 1.0580 1.3150 1.7060 2.0560 2.4790 2.7790 3.0670 3.4350 3.7070;
0.68400 0.85500 1.0570 1.3140 1.7030 2.0520 2.4730 2.7710 3.0570 3.4210 3.6900;
0.68300 0.85500 1.0560 1.3130 1.7010 2.0480 2.4670 2.7630 3.0470 3.4080 3.6740;
0.68300 0.85400 1.0550 1.3110 1.6990 2.0450 2.4620 2.7560 3.0380 3.3960 3.6590;
0.68300 0.85400 1.0550 1.3100 1.6970 2.0420 2.4570 2.7500 3.0300 3.3850 3.6460;
0.68100 0.85100 1.0500 1.3030 1.6840 2.0210 2.4230 2.7040 2.9710 3.3070 3.5510;
0.67900 0.84900 1.0470 1.2990 1.6760 2.0090 2.4030 2.6780 2.9370 3.2610 3.4960;
0.67900 0.84800 1.0450 1.2960 1.6710 2.0000 2.3900 2.6600 2.9150 3.2320 3.4600;
0.67800 0.84600 1.0430 1.2920 1.6640 1.9900 2.3740 2.6390 2.8870 3.1950 3.4160;
0.67700 0.84500 1.0420 1.2900 1.6600 1.9840 2.3640 2.6260 2.8710 3.1740 3.3900;
0.67700 0.84500 1.0410 1.2890 1.6580 1.9800 2.3580 2.6170 2.8600 3.1600 3.3730;
0.67400 0.84200 1.0360 1.2820 1.6450 1.9600 2.3260 2.5760 2.8070 3.0900 3.2910;];
if sides == 1
tval = interp2(one,tabledf, t_table, confint, df);
elseif sides == 2
tval = interp2(two,tabledf, t_table, confint, df);
else
error('must be either 1 or 2sided')
end
end
