function estvar = gsj(x,y)
%function estvar = gsj(x,y)
%Gasser-Sroka-Jennen-Steinmetz nonparametric variance estimator
%Arguments:
%   x: Design points (n,m)
%   y: Function observations (n,m)
%   Here n is the sample size and m is the number of functions or datasets

[n,m] = size(x);
if n==1
    x = x';
    y = y';
    [n,m] = size(x);
end
if n<3, return, end

[x,ind] = sort(x);
for j = 1:m, y(:,j) = y(ind(:,j),j); end
a = (x(3:n,:)-x(2:n-1,:))./(x(3:n,:)-x(1:n-2,:));
b = (x(2:n-1,:)-x(1:n-2,:))./(x(3:n,:)-x(1:n-2,:));
c = a.^2 + b.^2 + 1;
e = a.*y(1:n-2,:) + b.*y(3:n,:) - y(2:n-1,:);
estvar = mean(e.^2./c);