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
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 _______________
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
==============================================
> 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 =====================
> > > 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) ?