function cluster_analysis(wfb, i)
  % What are we clustering
  label=dimension_label(i);
  disp(sprintf('Clustering dimension %s ',label));
  
  % Prepare data: remove NaN's
  wfbi = wfb(:,i);
  num_countries = length(wfbi);
  wfbi(isnan(wfbi)) = [];
  number_of_nans = num_countries - length(wfbi);
  disp(sprintf('Threw away %d NaNs', number_of_nans));
 
  % Plot data
  %% Histogram
  subplot(2,2,1);
  hist (wfb(:,i),20);
  title(label);
  %% Sorted Data with various descriptive statistic
  wfbi = sort(wfbi);
  mean_wfbi = mean (wfbi);
  min = min (wfbi);
  max = max (wfbi);
  num_datapoints = length(wfbi);
  std_deviation = std (wfbi);
  %% Draw lines for mean +/-1 STD dev
  subplot (2,2,2);
  plot (wfbi);
  line ([0 num_datapoints], [mean_wfbi mean_wfbi], 'Color', 'r', 'LineWidth', 2);
  line ([0 num_datapoints], [mean_wfbi+std_deviation mean_wfbi+std_deviation]);
  line ([0 num_datapoints], [mean_wfbi-std_deviation mean_wfbi-std_deviation]);
 
  % Do K-means  
  %% Three Clusters
  [cluster_id, centroids, sumdistances] = kmeans (wfbi, 3, 'rep', 5, 'EmptyAction', 'drop');
  subplot (2,2,3);
  hold on;
  for i=1:length(wfbi)
   switch cluster_id(i)
   case 1
      plot (i,wfbi(i), 'b.');
   case 2
      plot (i,wfbi(i), 'g.');
   case 3
      plot (i,wfbi(i), 'r.');
   % draw percentile lines
    for i=1:2
      percentile_x = round(i*num_datapoints/3);
      line([0 num_datapoints], [wfbi(percentile_x) wfbi(percentile_x)]);
    end
   end
  end
  hold off;


  %% Five Clusters
  wfbi = sort(wfbi);
  [cluster_id, centroids, sumdistances] = kmeans (wfbi, 5, 'rep', 5, 'EmptyAction', 'drop');
  subplot (2,2,4);
  hold on;
  for i=1:length(wfbi)
   switch cluster_id(i)
   case 1
      plot (i,wfbi(i), 'b.');
   case 2
      plot (i,wfbi(i), 'g.');
   case 3
      plot (i,wfbi(i), 'r.');
   case 4
      plot (i,wfbi(i), 'm.');
   case 5
      plot (i,wfbi(i), 'c.');
   % draw percentile lines
    for i=1:4
      percentile_x = round(i*num_datapoints/5);
      line([0 num_datapoints], [wfbi(percentile_x) wfbi(percentile_x)]);
    end
   end
  end
  hold off;


function label=dimension_label(i)
%all the labels
dimension_labels = strvcat('total area', 'land area', 'total borders', 'arable land', 'population', 'population growth rate', 'birth rate', 'death rate', 'net migration rate', 'infant mortality rate', 'life expectancy at birth total', 'literacy total population', 'total fertility rate', 'total labor force', 'national product', 'national product real growth rate', 'national product per capita', 'unemployment rate', 'budget expenditures', 'exports total', 'imports total', 'external debt', 'electricity capacity', 'electricity production', 'electricity consumption per capita', 'highways total', 'airports total', 'television broadcast stations', 'defense expenditures number');
label=dimension_labels(i,:);
