function r=ksrmvAICBIC(x,y,criteria, z)
% KSRMV   Multivariate kernel smoothing regression
% criteria : AIC  (criteria=1) or BIC (criteria=2)
% r=ksrmv(x,y) returns the Gaussian kernel regression in structure r such that
%   r.f(r.x) = y(x) + e
% The bandwidth and number of samples are also stored in r.h and r.n
% respectively.
%
% r=ksrmv(x,y,h) performs the regression using the specified bandwidth, h.
%
% r=ksrmv(x,y,h,z) calculates the regression at location z (default z=x).
%
% Algorithm
% The kernel regression is a non-parametric approach to estimate the
% conditional expectation of a random variable:
%
% E(Y|X) = f(X)
%
% where f is a non-parametric function. Based on the kernel density
% estimation, this code implements the Nadaraya-Watson kernel regression
% using the Gaussian kernel as follows:
%
% f(x) = sum(kerf((x-X)/h).*Y)/sum(kerf((x-X)/h))
%
% See also gkde, ksdensity, ksr

% Check input and output
% Example 1:
% L = 40*membrane(1,25)+randn(51);
% [x,y]=meshgrid(0:50);
% x = [x(:) y(:)];
% y = L(:);

y=y(:);
if size(x,1)~=size(y,1)
    error('x and y have different rows.');
end
d=size(x,2);


r.x=x;
N=size(x,1);

% clean missing or invalid data points
inv=(y~=y);
x(inv,:)=[];
y(inv)=[];
r.n=numel(y);


% optimal bandwidth suggested by Bowman and Azzalini (1997) p.31
hy=median(abs(y-median(y)))/0.6745*(4/(d+2)/r.n)^(1/(d+4));
hx=median(abs(x-repmat(median(x),r.n,1)))/0.6745*(4/(d+2)/r.n)^(1/(d+4));
hx=sqrt(hy*hx);



seq = linspace(0.0001.*hx(1), hx(1), 100);
AIC=zeros(100,1);
BIC=zeros(100,1);
 for i= 1:100
    h = seq(i);
    hx = h*ones(1,d);
       
    % Scaling first
    H=diag(1./hx);
    newx=x*H;
    x1=r.x*H;

    % Gaussian kernel function
    kerf=@(x)exp(-sum(x.*x,2)/2);

    % allocate memory
    r.f=zeros(N,1);
    trace = 0;
    % Loop through each regression point
    for k=1:N
        % scaled deference from regression point
        xx=abs(newx-x1(k+zeros(r.n,1),:));
        % kernel function
        tr=kerf(xx);
        % regression
        r.f(k)=sum(tr.*y)/sum(tr);
        trace = trace + tr(k)/sum(tr);
    end
    RSS = sum((r.f-y).^2);
    AIC(i) = -2.*log(RSS)+2.*trace;
    BIC(i) = -2.*log(RSS)+trace.*log(N); 
    %[i h log(RSS) trace]
 end
 
     
if criteria==1
   % criteria is AIC
   value = AIC;    
elseif criteria==2
   % criteria is BIC 
   value = BIC;
end
%delete inf and NaN
idx = -(isinf(value)+isnan(value)-1);
idx = logical(idx);
%locate the minimum numeric value
idx = (value==min(value(idx)));
h = seq(idx);
hx = h*ones(1,d);
r.h = hx;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% testing smoothing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r.x=z;
N=size(z,1);
% Scaling first
H=diag(1./hx);
x=x*H;
x1=r.x*H;

% allocate memory
r.f=zeros(N,1);

% Loop through each regression point
for k=1:N
    % scaled deference from regression point
    xx=abs(x-x1(k+zeros(r.n,1),:));
    % kernel function
    z=kerf(xx);
    % regression
    r.f(k)=sum(z.*y)/sum(z);
end
