%% Script file to show regression using the EDM:
% First, construct some data to do the regression on.  In this case,
%  the domain will be 30 points in the plane (randomly selected).  The
%  range will be the function $e^{(-x^2-y^2)/4}$

X=randn(30,2);
Y=exp(-(X(:,1).^2+X(:,2).^2)./4);

% Be sure you have downloaded the file: edm.m from the class website
A=edm(X,X);
alpha=inv(A)*Y;

Z=A*alpha;
plot3(X(:,1),X(:,2),Y,'r*',X(:,1),X(:,2),Z,'b^');
hold on

% Create a 3-d surface to help visualize the function

% This function creates an array of points in the XY plane for the surface
[X1,Y1]=meshgrid(-3:0.125:3);
[m,n]=size(X1);

% Evaluate the function at x=X1(j,k) and y=Y1(j,k) created earlier
for j=1:m
    for k=1:n
        Z(j,k)=0;
        for t=1:30
            Z(j,k)=Z(j,k)+alpha(t)*norm([X1(j,k) Y1(j,k)]-X(t,:));
        end
    end
end
Z2=exp((-X1.^2-Y1.^2)/4);

% Surface command:
meshc(X1,Y1,Z);
hold off
figure(2)
meshc(X1,Y1,Z2-Z);
title('Error Between Surfaces');