Path: news.mathworks.com!newsfeed-00.mathworks.com!kanaga.switch.ch!switch.ch!newsfeed00.sul.t-online.de!t-online.de!news.nask.pl!news.nask.org.pl!news.onet.pl!not-for-mail
From: Zebbik <zebik@op.pl>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Linear regression.
Date: Tue, 18 Mar 2008 14:01:56 +0000
Organization: Onet.pl
Lines: 110
Message-ID: <froi0n$e0s$1@news.onet.pl>
References: <fr9c9t$o1p$1@news.onet.pl> <frbt1b$9nn$1@fred.mathworks.com>
NNTP-Posting-Host: acca5a67.ipt.aol.com
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
X-Trace: news.onet.pl 1205848920 14364 172.202.90.103 (18 Mar 2008 14:02:00 GMT)
X-Complaints-To: usenet@news.onet.pl
NNTP-Posting-Date: Tue, 18 Mar 2008 14:02:00 +0000 (UTC)
X-Sender: Bp2tzS8CFS8FGDaD6qs7Hg==
In-Reply-To: <frbt1b$9nn$1@fred.mathworks.com>
User-Agent: Thunderbird 2.0.0.12 (Windows/20080213)
Xref: news.mathworks.com comp.soft-sys.matlab:457830



Yuri Geshelin wrote:
> 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;
> %***********************************
> 
thanks for that =))