function [a,b,C,m]=quadreg(Y,X,z)
% fits a quadratic regression model by least squares of the form 
% Y=a+(x-z)b+0.5*(x-z)C(x-z)'  in the multivariate case, X is n by p, C is p by p etc.
% 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); X=X-Z; w=sum(X.^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
%w=ones(1,length(w));  % replace by vector of ones-it wasn't working well
%D=design matrix,  number of parameters =1+p+p(p+1)/2=
D=ones(n,1+p+p*(p+1)/2);
D(:,2:(p+1))=X;   colnum=zeros(p,p);
m=p+2;
for i=1:p
    D(:,m)=X(:,i).^2;
    colnum(i,i)=m;
    m=m+1;
    for j=(i+1):p
    D(:,m)=2*X(:,i).*X(:,j);
    colnum(i,j)=m;
    m=m+1;
end
end
par=inv(D'*diag(w)*D)*D'*diag(w)*Y;
a=par(1);
b=par(2:p+1);
for i=1:p
    for j=i:p
        C(i,j)=par(colnum(i,j));
        C(j,i)=C(i,j);
    end
end
m=-.5*inv(C)*b;  % this is the max or min of the quadratic
C=2*C;