clear,clc
x1 = xlsread ('data.xls','B2:K2');
x2 = xlsread ('data.xls','A3:A7');
y = xlsread ('data.xls','B3:K7');
n1=length(x1);n2=length(x2);
x1=ones(n2,1)*x1;x1=x1(:);
x2=x2*ones(1,n1);x2=x2(:);
y=y(:); X=[x1,x2]; n=length(y);
str=num2str([1:n]');
fx1=@(c,x1,x2)(c(1)+c(2)*x1+c(3)*x2+c(4)*x1.^2+c(5)*x1.*x2+c(6)*x2.^2+c(7).*x1.^3+c(8)*x2.*x1.^2+c(9)*x1.*x2.^2+c(10)*x2.^3);
fx2=@(c,x1,x2)(c(1)+c(2)*X(:,1)+c(3)*X(:,2)+c(4)*X(:,1).^2+c(5)*X(:,1).*X(:,2)+c(6)*X(:,2).^2+c(7)*X(:,1).^3+c(8)*X(:,2).*X(:,1).^2+c(9)*X(:,1).*X(:,2).^2+c(10)*X(:,2).^3);
c=[8.29927747717989E+03 2.82383157394932E+02 -5.89576900114284E+01 4.10191749116755E+00 -1.50873655516070E+00 -4.11074321839960E-02 2.14355866374489E-02 -2.76297884478486E-02 -6.66250572642417E-03 -8.28569814231031E-04];
for l=1:5
c=lsqcurvefit(fx2,c,X,y);
c=nlinfit(X,y,fx2,c);
end
c
figure(1),clf
plot3(x1,x2,y,'o')
stem3(x1,x2,y,'filled')
text(x1,x2,y+.01,str)
hold on
[x11,x22]=meshgrid(min(x1):range(x1)/80:max(x1),min(x2):range(x2)/80:max(x2));
yhat=fx1(c,x11,x22);
surf(x11,x22,yhat)
shading interp
alpha(.8)
axis tight
xlabel('X1');ylabel('X2'),zlabel('Y')
SSy=var(y)*(n-1)
y1=fx1(c,x1,x2);
RSS=(y-y1)'*(y-y1)
rsquare=(SSy-RSS)/SSy
MSe=RSS/(n-length(c))
c=vpa(c,10)