Weighted regression

351 views
Skip to first unread message

Gerard Middleton

unread,
Dec 10, 1996, 3:00:00 AM12/10/96
to

Is there anything written in MATLAB that does weighted regression?
--
Gerry Middleton
Department of Geology, McMaster University
Tel: (905) 525-9140 ext 24187 FAX 522-3141

W Tych

unread,
Dec 11, 1996, 3:00:00 AM12/11/96
to Gerard Middleton

Gerard Middleton wrote:
>
> Is there anything written in MATLAB that does weighted regression?


Do you mean the classic weighted LS solution:

a=(U'*W*U)\(U'*W*y);

where U is the regressors matrix, and W is a square, positive
defined matrix of weights (say with measurement weights on the
diagonal) ?

Example:
a=[1;2;3]; % assume a 3 parameters model
U=[ones(100,1), randn(100,1) randn(100,1)];
y=U*a; % and simulate it
W=eye(100); % attaching equal unit weights to all measurements
% I remember there are some addidtional conditions on W,
% you may have to look it up.

aa=(U'*W*U)\(U'*W*y) % re-estimate the parameters

aa1=U\y; % and compare with non-weighted LS:


Sorry if this is obvious and you are looking for something else.

Wlodek Tych

--
__________________Dr. Wlodzimierz (Wlodek) Tych___________________
| Control and Systems Group | phone: +44 1524 593973/593894 |
| IEBS/CRES,Lancaster University | fax: +44 1524 593985 |
| Lancaster LA1 4YQ U.K. | e-mail: w.t...@lancaster.ac.uk |
|__________WWW page: http://cres1.lancs.ac.uk/staff/wtych__________|
____ Lottery is a tax on people who did not listen to what the ____
__________ maths teacher was saying _______________

Matt Silvia

unread,
Dec 11, 1996, 3:00:00 AM12/11/96
to

W Tych wrote:
>
> Gerard Middleton wrote:
> >
> > Is there anything written in MATLAB that does weighted regression?
>
> Do you mean the classic weighted LS solution:
>
> a=(U'*W*U)\(U'*W*y);
>
> where U is the regressors matrix, and W is a square, positive
> defined matrix of weights (say with measurement weights on the
> diagonal) ?
>
> Example:
> a=[1;2;3]; % assume a 3 parameters model
> U=[ones(100,1), randn(100,1) randn(100,1)];
> y=U*a; % and simulate it
> W=eye(100); % attaching equal unit weights to all measurements
> % I remember there are some addidtional conditions on W,
> % you may have to look it up.
>
> aa=(U'*W*U)\(U'*W*y) % re-estimate the parameters
>
> aa1=U\y; % and compare with non-weighted LS:
>

Similarly, You can get get POLYFIT to create a fit in a "weighted"
fashion by modifying the M-file.

What you would need to do is first modify the first line of POLYFIT.M so
that it accepts a weighting matrix W as a fourth input argument. This
matrix W should be m-by-m where m is the length of x. It should also be
nonsingular and diagonal.

Next you will need to place the following lines after the section of the
code commented as "% Construct Vandermonde matrix":

V=W*V;
y=W*y;

Hope this helps!

--
==============================================
Matt Silvia sil...@mathworks.com
The MathWorks, Inc. http://www.mathworks.com
==============================================

Bradley Jones

unread,
Dec 13, 1996, 3:00:00 AM12/13/96
to

In article <32AEB461...@lancaster.ac.uk>, W.T...@lancaster.ac.uk (W
Tych) wrote:

> Gerard Middleton wrote:
> >
> > Is there anything written in MATLAB that does weighted regression?
>
>
> Do you mean the classic weighted LS solution:
>
> a=(U'*W*U)\(U'*W*y);
>
> where U is the regressors matrix, and W is a square, positive
> defined matrix of weights (say with measurement weights on the
> diagonal) ?
>
> Example:
> a=[1;2;3]; % assume a 3 parameters model
> U=[ones(100,1), randn(100,1) randn(100,1)];
> y=U*a; % and simulate it
> W=eye(100); % attaching equal unit weights to all measurements
> % I remember there are some addidtional conditions on W,
> % you may have to look it up.
>
> aa=(U'*W*U)\(U'*W*y) % re-estimate the parameters
>
> aa1=U\y; % and compare with non-weighted LS:
>
>

> Sorry if this is obvious and you are looking for something else.
>
> Wlodek Tych

The answer above is correct, but I thought it might be useful to show
another method that involves somewhat fewer calculations and has better
numeric properties. Here is an illustrative (I hope!) example:

X = [ones(10,1) rand(10,3)];
W = diag(rand(10,1));

% Note: the y data generated below are best estimated using weighted regression
y = X*(1:4)' + W*normrnd(0,0.1,10,1);


beta1 = (X'*W*X)\X'*(W*y); % These are Tych's estimates


beta = (sqrt(W)*X)\(sqrt(W)*y); % These are the estimates I am proposing.

beta1-beta

ans =

1.0e-14 *

-0.5551
0.7105
0.0888
-0.2220


The difference is a little larger than machine epsilon because of the loss
of precision in beta1 due to the computation of (X'*W*X)

--
==== Brad Jones ========================= br...@mathworks.com ======
| The MathWorks, Inc. in...@mathworks.com |
| 24 Prime Park Way http://www.mathworks.com |
| Natick, MA 01760-1500 ftp.mathworks.com |
==== Tel: 508-653-1415 ==== Fax: 508-653-6971 =====================

W Tych

unread,
Dec 19, 1996, 3:00:00 AM12/19/96
to br...@mathworks.com

Bradley Jones wrote:

> > > Is there anything written in MATLAB that does weighted regression?
> >
> >
> > Do you mean the classic weighted LS solution:
> >
> > a=(U'*W*U)\(U'*W*y);

[...]


> The answer above is correct, but I thought it might be useful to show
> another method that involves somewhat fewer calculations and has better
> numeric properties. Here is an illustrative (I hope!) example:
>
> X = [ones(10,1) rand(10,3)];
> W = diag(rand(10,1));
>
> % Note: the y data generated below are best estimated using weighted regression
> y = X*(1:4)' + W*normrnd(0,0.1,10,1);
>
> beta1 = (X'*W*X)\X'*(W*y); % These are Tych's estimates
>
> beta = (sqrt(W)*X)\(sqrt(W)*y); % These are the estimates I am proposing.

For completeness we should also mention that the difference between
the two is exactly as the one between the 'textbook' standard regression
(X'*X)\(X'*y),
the weighted equivalent of which I showed, and the Matlab shortened
X\y

Obviously, the shorter one should be used, the longer notation
although more 'educational', because it refers to the textbook solution
(that's why I used it), is also numerically inferior, I have clearly
missed the nifty one in my response. It is also very 'Matlabish'.

I believe that an even better, factorised solution can be suggested.
Any givers anyone ?

Thanks for the continuation of the thread.

Best wishes,

Wlodek Tych

P.S: (1) Brad, Tych is my surname, friends and colleagues call me
Wlodek.
No offence, just for info.
(2) Can you please generate the posted examples without using
toolboxes,
where it is possible ? It is a bit unfair, I can't replicate
your result. Is normrnd(0,0.1,10,1) the same as
0.1*randn(10,1),
or is it sqrt(0.1)*randn(10,1) ?

Reply all
Reply to author
Forward
0 new messages