351 views

Skip to first unread message

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

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 _______________

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

==============================================

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 =====================

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

Search

Clear search

Close search

Google apps

Main menu