function [out, e1] = testcv(x, resi)

% a code of SCV method in Xia (2006) for checking parametric and
% semiparametric models
%
% input 
%       x (nxp): the covariates of the fitted model
%       resi (nx1): the fitted residuals
% output
%       out: 1 indicates the model is not adequate
%            0 indicatss the model is adequate
%        e1: SCV value
%
% Example (see Stute et al 1999. If a = 0, a linear regression model is adequate)
%   n = 100; p = 3; a = 25;
%   x = rand(n,p);
%   y = 2+5*x(:,1) - x(:,2) + a*x(:,1).*x(:,2) + 1.0*randn(n,1);
%   x1 = [ones(n,1) x];
%   b = inv(x1'*x1)*x1'*y;   
%   resi = y-x1*b;
%   out = SCV(x, resi)
% 

   [n,p] = size(x);
   x=(x-repmat(mean(x),n,1))./repmat(std(x),n,1);
   v = inv(x'*x/n);
   x = x*v^0.5;
   resi = (resi-mean(resi))/std(resi);
  
   e0 = mean( (resi - mean(resi)).^2 );
   e1 = 0; e2 = 0;
   for i = 1:n
       Ii = [1:(i-1) (i+1):n];
       xi = x(Ii, :);
       yi = resi(Ii);
       ye = CVsim2c(xi, yi, x(i,:));   
       e1 = e1 + mean((ye-resi(i)).^2);
   end
   e1 = e1/n;
   out = 0 + (e0 > e1);

return
   

function ye = CVsim2c(x, y, xp);
[n,p]=size(x);

onep = ones(p,1);
onen = ones(n,1);
B = eye(p);

h2 = 2*n^(-2/(p+4));
p1 = 2:(p+1);
onexi = ones(n, p+1);
e = eye(p+1)/n;

onex = [ones(n*n, 1) repmat(x, n, 1) - kron(x, ones(n,1))];
for iter = 1:p+5;
	ab = ones(p,n);
    kernel = exp(-sum((onex(:,p1)*B).^2,2)/h2);
    xk = (onex.*repmat(kernel, 1, p+1))';
    for i = 1:n;
        I = (i-1)*n+1:i*n;
        abi = inv(xk(:,I)*onex(I, :) + e)*(xk(:,I)*y);
        ab(:,i) = abi(p1);
    end;
	[B0 D] = eig(ab*ab');
	[D I] = sort(diag(D));
	B = B0;
    B0 = B(:,I);
    B = B0(:, p+1-(1:p));
    nd = max(1, p-iter);
    B = B(:,1:nd);
    h = mean(std(x*B))/n^(1/(nd+4));
    h2 = 2*h*h;
end;

xB = x*B;
xp = xp*B;
[I, i] = sort(abs(xB-xp));
h = max(I(5), mean(std(xB))/n^0.225);
h2 = 2*h*h;

xi = xB - xp;
ker = exp(-sum(xi.^2,2)/h2);
ker1 = xi.*ker;
s1 = sum(ker1); s2 = xi'*ker1;
ker = (ker.*s2 - ker1.*s1);
ye = ker'*y./(sum(ker)+1/n/n);