Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
vectorized rolling regression?

Subject: vectorized rolling regression?

From: Michael Robbins

Date: 16 May, 2003 08:06:38

Message: 1 of 8

Has anyone come up with a vectorized rolling linear regression or
robust linear regression?


WINDOW=20;
[L,W]=size(y);
for i=L:-1:20
   X=[ones(size(x(i-20:i))) X(i-20:i)];
   Y=y(i-20:i);
   [b{i},bint{i},r{i},rint{i},stats{i}] = ...
   regress(Y,X);
   [br{i},statsr{i}] = robustfit(X,Y);
end;


Or even better one that will vectorize this


WINDOW=20;
[L,W]=size(y);
for i=1:L
   for j=1:W
       X=[ones(size(x(i-20:i,j))) X(i-20:i,j)];
       Y=y(i-20:i,j);
       [b{i,j},bint{i,j},r{i,j},rint{i,j},stats{i,j}] = ...
   regress(Y,X);
   [br{i,j},statsr{i,j}] = robustfit(X,Y);
end;


What about just for either regress or robustfit if you don't know how
to do both?

Subject: vectorized rolling regression?

From: John D'Errico

Date: 16 May, 2003 09:11:47

Message: 2 of 8

In article <eebdf1d.-1@WebX.raydaftYaTP>,
 "Michael Robbins" <michael.robbins@us.cibc.com> wrote:

> Has anyone come up with a vectorized rolling linear regression or
> robust linear regression?

This would be easy to write as a filter, at least for the
linear regression case. Something like this:

X = [ones(21,1),(-20:0)'];

filtercoef=pinv(X);

intercept_coef=filter(filtercoef(1,:),1,y);
slope_coef=filter(filtercoef(2,:),1,y);

HTH,
John D'Errico

Subject: vectorized rolling regression?

From: Michael Robbins

Date: 16 May, 2003 11:24:12

Message: 3 of 8

> This would be easy to write as a filter, at least for the
> linear regression case. Something like this:
>
> X = [ones(21,1),(-20:0)'];
>
> filtercoef=pinv(X);
>
> intercept_coef=filter(filtercoef(1,:),1,y);
> slope_coef=filter(filtercoef(2,:),1,y);
>
> HTH,
> John D'Errico


Thanks again for your help. Maybe I coded it wrong, but the results
are quite different, though much faster;


%-----------------------------------
function [is,st,t]=RollingRegress(x,y,WINDOW,TEST)


if nargin<2
   x=[1:1000]';
   y=x+rand(1000,1);
end;
if nargin<3
   WINDOW=20;
end;
if nargin<4
   TEST=1;
end;


if TEST tic; end;
[L,W]=size(y);
for i=L:-1:(WINDOW+1)
   X=[ones(size(x(i-20:i))) x(i-20:i)];
   Y=y(i-20:i);
   [b{i},bint{i},r{i},rint{i},stats{i}] = ...
   regress(Y,X);
   %[br{i},statsr{i}] = robustfit(X,Y);
end;
if TEST t.loop=toc; end;
is.loop=[b{:}].';
st.stats=stats;
st.bint=bint;
st.r=r;
st.rint=rint;


if TEST tic; end;
X = [ones(21,1),(-20:0)'];
filtercoef=pinv(X);
intercept_coef=filter(filtercoef(1,:),1,y);
slope_coef=filter(filtercoef(2,:),1,y);
if TEST t.filter=toc; end;
is.filter=[intercept_coef slope_coef];

Subject: vectorized rolling regression?

From: Denis Gilbert

Date: 16 May, 2003 16:41:55

Message: 4 of 8

On Fri, 16 May 2003 11:24:12 -0400, "Michael Robbins"
<michael.robbins@us.cibc.com> wrote:


>Thanks again for your help. Maybe I coded it wrong, but the results
>are quite different, though much faster;

Michael,

I tried John's code and obtained good results by changing

 X = [ones(21,1), (-20:0)'];

to

 X = [ones(21,1),(0 : -1 : -20)'];

You will also need to decide whether you want to have your slope
estimates centered or not about your rolling window width. By that I
mean that if your time unit is years and your last data point
corresponds to year 2002, do you want your last slope estimate to be
in year 2002 or centered at 1992, the middle of your latest rolling
window? In the latter case you will need to shift back the results by
windowWidth/2 indices and set to NaN the last windowWidth/2 indices of
your slope estimate vector.

Hope this helps, Denis.

Subject: vectorized rolling regression?

From: John D'Errico

Date: 16 May, 2003 13:09:35

Message: 5 of 8

In article <vg4acvsea1tin33d2aomchq2j97e54aski@4ax.com>,
 Denis Gilbert <nospam@thank.you> wrote:

> On Fri, 16 May 2003 11:24:12 -0400, "Michael Robbins"
> <michael.robbins@us.cibc.com> wrote:
>
>
> >Thanks again for your help. Maybe I coded it wrong, but the results
> >are quite different, though much faster;
>
> Michael,
>
> I tried John's code and obtained good results by changing
>
> X = [ones(21,1), (-20:0)'];
>
> to
>
> X = [ones(21,1),(0 : -1 : -20)'];
>
> You will also need to decide whether you want to have your slope
> estimates centered or not about your rolling window width. By that I
> mean that if your time unit is years and your last data point
> corresponds to year 2002, do you want your last slope estimate to be
> in year 2002 or centered at 1992, the middle of your latest rolling
> window? In the latter case you will need to shift back the results by
> windowWidth/2 indices and set to NaN the last windowWidth/2 indices of
> your slope estimate vector.
>
> Hope this helps, Denis.
>

Thanks Denis.

I was sloppy in creating my filter coefficients.

John

Subject: vectorized rolling regression?

From: Michael Robbins

Date: 16 May, 2003 13:48:21

Message: 6 of 8

>>>Thanks again for your help. Maybe I coded it wrong, but the
>>> results are quite different, though much faster;


>> I tried John's code and obtained good results


Thanks to both of you. I'm not sure what you mean by good. Look at
what I computed:


K>> [is,st,t]=RollingRegress;
fprintf(' Loop Beta Intcpt Filter Beta Intcpt\n');
format short;
i=1:50:900;
[is.loop(i,:) is.filter(i,:)]


   Loop Beta Intcpt Filter Beta Intcpt
    0.6174 0.9981 0.3461 0.0253
    0.1782 1.0063 51.4125 0.9920
    0.8881 0.9963 101.3951 0.9870
   -1.7105 1.0126 151.4827 0.9959
   -1.1928 1.0078 201.4332 0.9977
   -0.2148 1.0030 251.4869 0.9971
   -3.1599 1.0120 301.5340 0.9958
    0.1821 1.0010 351.4433 0.9950
   -8.8355 1.0228 401.4653 1.0037
    4.8185 0.9906 451.4497 0.9902
    0.7916 0.9994 501.3675 0.9906
    9.9444 0.9832 551.7281 1.0264
    2.6864 0.9964 601.5338 1.0042
    2.5644 0.9967 651.6866 1.0142
    1.0942 0.9992 701.5832 1.0025
   22.8601 0.9706 751.6297 1.0016
    5.9169 0.9934 801.7259 1.0193
  -19.7636 1.0236 851.3268 0.9885


%__-__-__-__-__-__-__-__-__-__-__-__-__-


function [is,st,t]=RollingRegress(x,y,WINDOW,TEST)


if nargin<2
   x=[1:1000]';
   y=x+rand(1000,1);
end;
if nargin<3
   WINDOW=20;
end;
if nargin<4
   TEST=1;
end;


if TEST tic; end;
[L,W]=size(y);
for i=L:-1:(WINDOW+1)
   X=[ones(size(x(i-20:i))) x(i-20:i)];
   Y=y(i-20:i);
   [b{i},bint{i},r{i},rint{i},stats{i}] = ...
   regress(Y,X);
   %[br{i},statsr{i}] = robustfit(X,Y);
end;
if TEST t.loop=toc; end;
is.loop=[b{:}].';
st.stats=stats;
st.bint=bint;
st.r=r;
st.rint=rint;


if TEST tic; end;
X = [ones(21,1),(0:-1:-20)'];
filtercoef=pinv(X);
intercept_coef=filter(filtercoef(1,:),1,y);
slope_coef=filter(filtercoef(2,:),1,y);
if TEST t.filter=toc; end;
is.filter=[intercept_coef slope_coef];

Subject: vectorized rolling regression?

From: Michael Robbins

Date: 16 May, 2003 13:50:12

Message: 7 of 8

> i=1:50:900;


That was i=1:20:500 because of the offset used (WINDOW)


I cut-and-pasted it incorrectly.

Subject: Help refine code

From: Ashish Ojha

Date: 19 Dec, 2004 21:58:39

Message: 8 of 8

Hi People below is my code.
I am trying to output my struct to a file along with the fieldnames as the
headers. I first try to read the 'class' of each field and their
dimension. If the dimesions are more than 1 then I modify that field name
by increasing the field-names along with element number appended to
it..for e.g
if the fieldname is 'go' and it has 2 elements then the header reads
go(1) go(2).
and then I build up the format by the classtype and #dimensions and then
print that out....so I have it completed can any1 of y'll refine it for me
and help vectorize it as well.Thx for your time.
OZ
==========================================================================
clear;clc;rand('seed',0)
data=struct('anum',{rand},'bfournum',{rand(1,4)},'cchar',{'ozo'},'dtwonum',{rand(2,1)});
data(2).anum=rand;
data(2).bfournum=rand(1,4);
data(2).cchar='bozo';
data(2).dtwonum=rand(2,1)
cdata=struct2cell(data);
datafn=num2cell(fieldnames(data));
cdatadouble=find(cellfun('isclass',cdata(:,:,1),'double'));
cdatachar=find(cellfun('isclass',cdata(:,:,1),'char'));
numericcdata=cdata(cdatadouble);
numericdatafn=cat(1,datafn{cdatadouble});
cdatadim=cellfun('length',numericcdata);
vfn=numericdatafn(cdatadim>1);
vcdatadim=cdatadim(cdatadim>1);
datafnm=datafn;
for i=1:length(vcdatadim)
    si(i)=strmatch(vfn(i),cat(1,datafn{:}));
    datafnm{si(i)}=strcat(vfn(i),'(',num2str([1:vcdatadim(i)]'),')');
end
headers=cat(1,datafnm{:});tmp=headers';
fprintf([repmat('%s\t',1,length(headers)-1) '%s\n'],tmp{:})
M=NaN*zeros(length(cdatadouble),max(cdatadouble));
for k=1:length(cdatadouble)
    M(k,1:cdatadim(k))=cdatadouble(k);
end
Mm=sort(M(:)');
MmM=Mm(~isnan(Mm));
fm=sort([MmM,cdatachar]);
fm1=[{'%f'} repmat({'\t%f'},1,length(fm)-1) {'\n'}];
fm1(fm==cdatachar)={'\t%s'};
fstr=[num2str(cat(2,fm1{:}))];
fprintf(fstr,cdata{:});
====================================================================

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us