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

How do I vectorize the loop for the local linear regression?

Asked by Julio on 24 Jul 2012

I am currently writing code for the project that my professor is thinking of, and this involves writing code for local linear regression and Nadaraya-Watson regression. The first code that I wrote uses a loop for the NW estimator and the local linear regression estimator.

% This code simulates the local linear estimator and NW kernel.
% Generating noisy data. (I use the same as the example in the Lo paper.)
 x=linspace(0,4*pi,100);
 y=sin(x)+0.5*randn(size(x));
% Set up the space to store the estimates.
 yhatnw=zeros(size(x));
 yhatll=zeros(size(x));
 n=length(x);
% Set up the bandwidth.
 hx=median(abs(x-median(x)))/0.6745*(4/3/n)^0.2;
 hy=median(abs(y-median(y)))/0.6745*(4/3/n)^0.2;
 h=sqrt(hy*hx);
% Find smooth at each value of x.
 for i=1:n
    w=wfun(h,x(i),x);
    xc=x-x(i);
    s2=mean(xc.^2.*w);
    s1=mean(xc.*w);
    s0=mean(w);
    yhatnw(i)=sum(w.*y)/sum(w); % Nadaraya-Watson kernel
    yhatll(i)=sum(((s2-s1*xc).*w.*y)/(s2*s0-s1^2))/n; % local linear estimator
 end
plot(x,y,'.',x,yhatnw,'-',x,yhatll,':')
The corresponding weighting function is:
function w=wfun(h,mu,x);
w=exp((-1/2)*(((x-mu)/h).^2))/sqrt(2*pi*h^2);

I tried to vectorize the code, and it works for the NW estimator (since I am getting the same graph), but not for the local linear estimator (the graph with the loop is better than that of the matrix). The relevant code (for the local linear estimator) is here:

% Evaluating the weighting matrix.
 xi=ones(n,1)*x; 
 data=x'*ones(1,n);
 w=wfun(h,xi,data);
% % Local linear regression
 xc=data-xi;
 s2=mean(xc.^2.*w);
 s1=mean(xc.*w);
 s0=mean(w);
 s2i=ones(n,1)*s2;
 s1i=ones(n,1)*s1;
 s0i=ones(n,1)*s0;
 yhatll=sum(((s2i-s1i*xc).*w.*yi)'./(s2i*s0i-(s1i^2)))/n;

I think there is something wrong with the way I am vectorizing the local linear regression, but I am unsure, so I'd like to ask all of you. Thanks!

2 Comments

bym on 24 Jul 2012

I think your post is missing the definition of yi, which sidetracked me for awhile. Please confirm that yi is defined something like

yi = repmat(y,100,1);
Julio on 25 Jul 2012

Hi, I was looking at it, and you are right. I was missing my definition of yi. I'll put it right here:

yi=ones(n,1)*y;
Julio

Products

No products are associated with this question.

1 Answer

Answer by Andrei Bobrov on 25 Jul 2012
Edited by Andrei Bobrov on 25 Jul 2012
xc = bsxfun(@minus,x',x);
w = exp(-.5*(xc/h).^2)/sqrt(2*pi*h^2);
xw = xc.*w;
s0 = mean(w);
s1 = mean(xw);
s2 = mean(xw.*xc);
yhatnw = sum(bsxfun(@times,w,y'))./sum(w);
yhatll = sum(...
          bsxfun(@rdivide,...
          bsxfun(@times,...
          bsxfun(@minus,s2,...
          bsxfun(@times,xc,s1)).*w,y'),s2.*s0 - s1.^2))/n;

or your variant

[ii,jj] = meshgrid(1:n);
xc = x(jj)-x(ii);
w=exp( -1/2*(xc/h).^2 )/sqrt(2*pi*h^2);
s2=mean(xc.^2.*w);
s1=mean(xc.*w);
s0=mean(w);
yhatll3=sum(((s2(ii)-s1(ii).*xc).*w.*y(jj))./(s2(ii).*s0(ii)-s1(ii).^2))/n;

0 Comments

Andrei Bobrov

Contact us