% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% A modified Example 5.7 - but using the 
%
load lungB
% Take the transpose, because the
% colums are the observations.
X = lungB';
[n,p] = size(X);
% Standardize the columns.
for i = 1:p
    X(:,i) = X(:,i)/std(X(:,i));
end

K = 10; 
B = 100; 
link_method = 'average'; 
pdist_method = 'euclidean'; 
[Z, khat, gap, Wobs, muWb] = gap_uniform(X,K,B,link_method,pdist_method); 
fprintf('khat = %10d\n',khat); 

figure; plot(1:K,Wobs,'o-',1:K,muWb,'x-')
legend({'Observed';'Expected'})
xlabel('Number of Clusters k')
ylabel('Observed and Expected log(W_k)')
saveas( gcf, '../../WriteUp/Graphics/Chapter5/prob_5_17_pca_oe_clust', 'epsc' ); 

figure,plot(1:K,gap,'o-'),title('Gap')
xlabel('Number of Clusters k'); ylabel('Gap')
saveas( gcf, '../../WriteUp/Graphics/Chapter5/prob_5_17_pca_gap_clust', 'epsc' );