function [a,b,tk]=linreg(Y,X,z)
% like quadrag only linear
% Y=a+Xb  in the multivariate case, X is n by p,
% z is a point (row vector) about which e wish to weight observations. weights decrease
% with distance from z
p=size(X,2); n=size(X,1);
Z=repmat(z,n,1); w=sum((X-Z).^2,2);  me=median(w); % choose tuning constant tk so that farthest away 50% receive weight<10%.
tk=9/me;   
w=1./(1+tk*w);  % w is the weight vector
%D=design matrix,  number of parameters =1+p
D=ones(n,1+p);
D(:,2:(p+1))=X;  
par=inv(D'*diag(w)*D)*D'*diag(w)*Y;
a=par(1);
b=par(2:p+1);