% Example 3.5
a=sqrt(2); b=2;
X=[1 -1 -1 1 1 1;
   1 1 -1 1 -1 1;
   1 1 1 1 1 1;
   1 -1 1 1 -1 1;
   1 a 0 b 0 0;
   1 -a 0 b 0 0;
   1 0 a 0 0 b;
   1 0 -a 0 0 b;
   %1 0 0 0 0 0;
   %1 0 0 0 0 0;
   %1 0 0 0 0 0;
   %1 0 0 0 0 0;
   1 0 0 0 0 0];
XTX=X'*X;
XTXi=inv(XTX);
[X, Y]=meshgrid(-1:.1:1, -1:.1:1);
Z=zeros(size(X));
[n, m]=size(X);
for i=1:n
 for j=1:m
  xi=[1 X(i,j) Y(i,j) X(i,j).^2 X(i,j).*Y(i,j) Y(i,j).^2];
  Z(i,j)=sqrt(xi*XTXi*xi');
 end
end
v=linspace(min(min(Z)),max(max(Z)),10);
[C,h]=contour(X,Y,Z,v);
clabel(C,h)
Zmax=max(max(Z))
Zmin=min(min(Z))
stability=Zmax/Zmin
%-----------------------------------------------------------------------------
%
% Figure 3.10
dCC = ccdesign(2,'type','circumscribed');
plot(dCC(:,1),dCC(:,2),'ro','MarkerFaceColor','b')
X = [1 -1 -1 -1; 1 1 1 -1];
Y = [-1 -1 1 -1; 1 -1 1 1];
line(X,Y,'Color','b')
axis square equal off
%-----------------------------------------------------------------------------
%
% Example 3.8
ny=6; %With 6 samples:
[dce,X]=cordexch(2,ny,'quadratic');
scatter(dce(:,1),dce(:,2),200,'filled')
D6=det(X'*X)
%
ny=12; %With 12 samples:
[dce,X]=cordexch(2,ny,'quadratic');
scatter(dce(:,1),dce(:,2),200,'filled')
D12=det(X'*X)
%-----------------------------------------------------------------------------
%
% Example 3.10
X=[0 0;1 0;1 1]; % D-optimal design
%X=[0 0;1 0;0 1]; % A-optimal design
%X=[0 0;1 0;0.5 1]; % G-optimal design
%
XTX=X'*X;
XTXi=inv(XTX);
[X, Y]=meshgrid(0:.1:1, 0:.1:1);
Z=zeros(size(X));
[n, m]=size(X);
for i=1:n
 for j=1:m
  xi=[X(i,j) Y(i,j)];
  Z(i,j)=sqrt(xi*XTXi*xi');
 end
end
v=linspace(min(min(Z)),max(max(Z)),10);
[C,h]=contour(X,Y,Z,v);
clabel(C,h)
Zmax=max(max(Z))
Zmin=min(min(Z))
%-----------------------------------------------------------------------------
%
% Example 3.12
rng default;
x=rand(20,2);
subplot(2,2,1); plot(x(:,1), x(:,2), 'o');
subplot(2,2,2); hist(x(:,2),20);
subplot(2,2,3); hist(x(:,1),20);
%-----------------------------------------------------------------------------
%
% Figure 3.17
x=lhsdesign(10,2); plot(x(:,1), x(:,2), 'o');
xr=lhsdesign(10,2,'criterion','correlation');
hold on; plot(xr(:,1), xr(:,2), 'r+');
r=corrcoef(x)
%r = 1.0000 -0.6999
% -0.6999 1.0000
r=corrcoef(xr)
%r = 1.0000 -0.0545
% -0.0545 1.0000
%-----------------------------------------------------------------------------