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

15 次查看（过去 30 天）
Julio
Julio 2012-7-24
评论： Carlo Grillenzoni 2024-11-25
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 个评论
bym
bym 2012-7-24
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
Julio 2012-7-25
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;

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Andrei Bobrov
Andrei Bobrov 2012-7-25
编辑：Andrei Bobrov 2012-7-25
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;
  3 个评论
显示 1更早的评论
Image Analyst
Image Analyst 2024-11-24
@Carlo Grillenzoni, what is "Nadaraya-Watson regression" and why do you want that specific denoising algorithm? There are lots of denoising algorithms for 2-D matrices/images. How about a simple averaging via convolution, or a median filter via stdfilt? Or there are lots of other algorithms like nonlocal means (there is a MATLAB function for that) or BM3D. See "A tour of modern image filtering": https://users.soe.ucsc.edu/~milanfar/publications/journal/ModernTour.pdf
Carlo Grillenzoni
Carlo Grillenzoni 2024-11-25
I just answered to the post of Andrei Bobrov (assuming he could receive it)
I need a simple MatLab implementation of Bivariate Local Linear Regression (LLR)
for sparse spatial points, and for my further methodological developments, CG

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear Regression 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by