% This script will visually illustrate how a couple of different methods of unsupervised clustering work.  Unsupervised clustering is a common task within data science.  The basic problem can be summed up as follows:
% Suppose you have a collection of N data points, each of which is described by some set of properties that can be put into the form of a numeric vector with m entries.  We will use x_i to denote the vector of properties for data point i.  Further suppose that each of these data points also belongs to one of K different fundamental classes of objects.  But, and this is important, you don't have any direct information about which of your data points correspond to which class of object, or even any other similar dataset that might tell you the type of properties that objects of a certain class possess.  This is what makes the problem unsupervised.  Despite this lack of information, you are tasked with determining the class that each of these data points belongs to; that is, you must classify each data point as one of the K types.  How can you do this?
% To get started, let's first just visualize an example of what we are really talking about here:

xa=mvnrnd([0,0],.2*eye(2),50); % example using data from three different multivariate normal distributions
xb=mvnrnd([1,1],.2*eye(2),50);
xc=mvnrnd([1,-1],.2*eye(2),50);
x=[xa;xb;xc];
scatter(xa(:,1),xa(:,2)) % visualize the data
hold on;
scatter(xb(:,1),xb(:,2))
scatter(xc(:,1),xc(:,2))
hold off

% In this example, we have N=150 data points, each of which is just a numeric vector with m=2 entries.  Each data point also belongs to one of three different classes, which are illustrated here with three different colors.  But it is important to remember that, while for visualization purposes here the colors of each point are shown, we have to pretend that we don't actually know the color of any of these datapoints, nor do we possess any other similar data set that would contain such colors.  If we didn't know these colors, could we nonetheless cluster our data into three distinct groups while retaining some semblence of accuracy?
% Before proceeding to discuss some techniques that can be used for this task, we note another potential difficulty we haven't point out yet.  While for some specific applications we might have knowledge of what the number of classes K is, in many applications we don't even know this!  And, even if we did know this, since we have no knowledge of the true classification of any of the data points, there is no good way to actually assess the accuracy of whatever classifications we end up making on our data.  These are some of the things that make unsupervised clustering a challenging problem, at least on a theoretical level.
%
% Having described the problem at hand, and the challenges associated with it, let us now discuss some popular methods that can be used.  The first method is called K-means clustering.  The fundamental justification for this method is roughly as follows: in many scenarios where unsupervised clustering might be relevant and plausible, we might expect that the different classes of objects inhabit somewhat distinct regions of the m dimensional space being used to describe them.  This is largely true of our numerical example above, where we can visually see that the data points with the same color are all naturally clustered together in some way in that two dimensional space, albeit with some overlapping regions.  Given this plausible assumption, we might also be tempted to believe that, in general, objects of the same class are "closer" together, on average, than objects of different classes.
% In order to proceed mathematically, let's stipulate that there exists some way to measure how "distant" two different data points are.  For certain types of data, a definition of distance might be obvious, while for others it might require some thinking.  But in any case, we will simply define the distance between two data points i and j as d(i,j), which we will insist upon being a non-negative scalar quantity.  We might also refer to d(i,j) as the "dissimilartiy" of two data points; we use this interchangably with distance from now on.
% Finally, let's use C_k to denote the set of all points belonging to category k, where k is any of the K categories, and use |C_k| to denote the number of objects in category k.  Then K-means takes the position that we should expect the average dissimilarity of objects within each category to be small, at least relative to the dissimilarity of objects in different categories.  One can compute the total average within-group dissimilarity as
% T_W=\sum_{k=1}^K \frac{1}{|C_k|} \sum_{i,j\in C_k} d(i,j)
% For simplicity, assume now that d(i,j) is just the squared Euclidean distance.  Then the formula for T_W simplifies to something with an interesting interpretation:
% T_W=2\sum_{k=1}^K \sum_{i\in C_k} (x_i -\bar{x}_k)^2
% where \bar{x}_k is the mean value of the data points within group k.  So, we can state then that K-means thinks each data point should be close to the mean of the group that it belongs to.  Let's compute these values for the groups in our example:

mna=mean(xa); % find the means of the three groups
mnb=mean(xb);
mnc=mean(xc);
diffa=xa-mna; % find the distance from each point to its group's mean
diffb=xb-mnb;
diffc=xc-mnc;
Twa=sum(sum(diffa.^2,2)); % find average within group distances
Twb=sum(sum(diffb.^2,2));
Twc=sum(sum(diffc.^2,2));
Tw=Twa+Twb+Twc; % total average within group distances

% As a point of comparison, let's also compute the total average between-group dissimilarity via
% T_B=\sum_{k=1}^K \frac{1}{N-|C_k|} \sum_{i\in C_k}\sum_{j\notin C_k} d(i,j)

alldists=squareform(pdist([xa;xb;xc],"squaredeuclidean")); % find all interpoint distances
distsa=alldists(1:50,51:150); % find distances between groups
distsb=alldists(51:100,[1:50,101:150]);
distsc=alldists(101:150,1:100);
Tba=mean(sum(distsa,2)); 
Tbb=mean(sum(distsb,2));
Tbc=mean(sum(distsc,2));
Tb=Tba+Tbb+Tbc;
fprintf('Total average between-group distance is %g, while within-group is %g.',Tb,Tw)

% We see, then, that for this particular dataset the K-means expectation of relatively low within-group distance holds true.
% So, how does K-means actually go about using this expected behavior to assign group membership to the datapoints?  Refer back to our formula above for T_W.  For a fixed value of K, if our algorithm assigns to each datapoint i a category C_k(i), the result will be some numerical value for the quantity T_W.  In other words, T_W is a function of the various category labels of the datapoints.  It is the assumption of K-means that, if data points are assigned to their correct categories, the value of this function should be small.  So, K-means turns this into an optimization problem: given our set of datapoints and a value for K, assign each datapoint i into a category C_k(i) in order to minimize T_W.
% Unfortunately, this minimization problem is a bit tricky to solve exactly.  So, typically K-means uses an approximate algorithm in order to try to solve it instead.  The algorithm works as follow:
% Begin by assigning each datapoint to some category in any way you like (including randomly), so long as each of the K categories has at least one point assigned to it.
% Based on these current assignments, compute the means of the current categories.
% Consider the category means as fixed for this step.  For each datapoint, update its category membership so that it belongs to the category with the closest mean.
% Go back to step 2.  Continue alternating steps 2 and 3 until no datapoints change their category in step 3.
% Here it is in code:

K=3; % number of categories
x=x(randperm(150),:); % initial random category assignment
assign=zeros(150,1);
for i=1:150
    assign(i)=mod(i-1,K)+1;
end
mns=zeros(K,2); % setting up variable to store the category means
assignold=assign; % setting up looping variables
assignold(1)=0;

while (sum(assignold~=assign)>0) % continue looping while some points are changing category
    for k=1:K % compute current means
        sub=x(assign==k,:);
        mns(k,:)=mean(sub);
    end
    assignold=assign; % store current categories for comparison
    for i=1:150 % update categories by assigning to closest mean
        diff=mns-x(i,:);
        diff=sum(diff.^2,2);
        [~,ind]=min(diff);
        assign(i)=ind;
    end
end

x1p=x(assign==1,:); % define K-means category sets
x2p=x(assign==2,:);
x3p=x(assign==3,:);
% end K-means.  Code below is for assessment and visualization only

ca1=sum(ismember(xa,x1p,'rows')); % assign each of these to the closest of the true categories for assessment purposes
xap=x1p;
camax=ca1;
ca2=sum(ismember(xa,x2p,'rows'));
if ca2>ca1
    xap=x2p;
    camax=ca2;
end
ca3=sum(ismember(xa,x3p,'rows'));
if ca3>camax
    xap=x3p;
    camax=ca3;
end
cb1=sum(ismember(xb,x1p,'rows'));
xbp=x1p;
cbmax=cb1;
cb2=sum(ismember(xb,x2p,'rows'));
if cb2>cb1
    xbp=x2p;
    cbmax=cb2;
end
cb3=sum(ismember(xb,x3p,'rows'));
if cb3>cbmax
    xbp=x3p;
    cbmax=cb3;
end
cc1=sum(ismember(xc,x1p,'rows'));
xcp=x1p;
ccmax=cc1;
cc2=sum(ismember(xc,x2p,'rows'));
if cc2>cc1
    xcp=x2p;
    ccmax=cc2;
end
cc3=sum(ismember(xc,x3p,'rows'));
if cc3>ccmax
    xcp=x3p;
    ccmax=cc3;
end

figure; % visualize, including plotting the means
scatter(xap(:,1),xap(:,2))
hold on;
scatter(xbp(:,1),xbp(:,2))
scatter(xcp(:,1),xcp(:,2))
scatter(mns(:,1),mns(:,2),[],'MarkerFaceColor',[0 0 0])
hold off

% Now, we can visually compare the K-means result to the true categories shown above. In general, it appears to have done well!
% We can be a little more precise by considering that accuracy of the K-means result, defined as the percentage of data points assigned to their correct true category:

acc=(camax+cbmax+ccmax)/150;
fprintf('The K-means accuracy was %g.',acc)

% We can also consider what the total average within-group distance was for our assigned categories:

mnap=mean(xap); % find the means of the three groups
mnbp=mean(xbp);
mncp=mean(xcp);
diffap=xap-mnap; % find the distance from each point to its group's mean
diffbp=xbp-mnbp;
diffcp=xcp-mncp;
Twap=sum(sum(diffap.^2,2)); % find average within group distances
Twbp=sum(sum(diffbp.^2,2));
Twcp=sum(sum(diffcp.^2,2));
Twp=Twap+Twbp+Twcp; % total average within group distances
fprintf('The K-means total average within-group distance was %g.',Twp)

% The result here will vary a little bit from one run of this code to the next (the original data is random, plus the K-means result may vary from one run to the next on a fixed dataset, based on the initial random group assignments), but in general you will likely find that the within-group distance for the K-means categories is smaller than the true distance.  This is of course because K-means is directly attempting to minimize this distance, so (assuming it succeeds in finding the global minimum) the result cannot be greater than the true within-group distance.
%
% In the implementation above, the correct value of K was used by K-means.  But, as noted previously, there is no guarantee that we will know the correct K value.  Let's see what happens when we provide a different K:

K=8;

%

allgroups=zeros(K,1);
while (sum(allgroups)<K) % make sure we end up with K groups
    x=x(randperm(150),:); % initial random category assignment
    assign=zeros(150,1);
    for i=1:150
        assign(i)=mod(i-1,K)+1;
    end
    mns=zeros(K,2); % setting up variable to store the category means
    assignold=assign; % setting up looping variables
    assignold(1)=0;

    while (sum(assignold~=assign)>0) % continue looping while some points are changing category
        for k=1:K % compute current means
            sub=x(assign==k,:);
            mns(k,:)=mean(sub);
        end
        assignold=assign; % store current categories for comparison
        for i=1:150 % update categories by assigning to closest mean
            diff=mns-x(i,:);
            diff=sum(diff.^2,2);
            [~,ind]=min(diff);
            assign(i)=ind;
        end
    end
    allgroups=ismember(1:K,assign); % see if all K groups are used
end

xp=x(assign==1,:); % visualize and diagnose, including plotting the means
diff=xp-mns(1,:);
Tw=sum(sum(diff.^2,2));
figure; 
scatter(xp(:,1),xp(:,2))
hold on;
for k=2:K
    xp=x(assign==k,:);
    diff=xp-mns(k,:);
    Tw=Tw+sum(sum(diff.^2,2));
    scatter(xp(:,1),xp(:,2));
end
scatter(mns(:,1),mns(:,2),[],'MarkerFaceColor',[0 0 0])
hold off
fprintf('The K-means total average within-group distance was %g.',Tw)

%
% Kmeans relies on a particular dissimilarity measure between points: Euclidean distance.  This simplifies the problem considerably, since the point that minimizes total Euclidean distance to a set of points is the mean of those points, which allows the formulation of Kmeans to simply use the mean of the cluster to assess within-group dissimilarity.  However, in some cases a different distance measure may be more appropriate, which may not have this property.  In these cases, we can alter Kmeans into Kmedoids.  Kmedoids still attempts to minimize total intra-group dissimilarity, measured via
% T_W=\sum_{k=1}^K  \sum_{i\in C_k} d(i,m_k).
% Here, m_k is the location of the medoid (some measure of the "center") of cluster k.  Ideally, this problem would be solved by optimizing both over cluster membership and medoid location simultaneously, but this is hard.  Instead, an approximate method is used (there are several such methods, but we present one of them here).  As a first step of this approximation, we will restrict the values of the various m_k so that they must each be located at one of the datapoints in our set.  Let i^*_k be the index of the datapoint that gives the location of m_k.  Then proceed as follows:
% Begin by assigning each datapoint to some category in any way you like (including randomly), so long as each of the K categories has at least one point assigned to it.
% Based on these current assignments, find i^*_k by locating the datapoint j that minimizes  \sum_{i\in C_k} d(i,j).
% Consider the category medoids as fixed for this step.  For each datapoint, update its category membership so that it belongs to the category with the closest medoid as measured via d(i,i^*_k).
% Go back to step 2.  Continue alternating steps 2 and 3 until no datapoints change their category in step 3.
% This is implemented below in code.  For the distance metric here, we use "cityblock" distance: for a vector (x_1,x_2,\ldots,x_p), the cityblock length is \sum_{i=1}^p |x_i|.

K=3; % number of categories

x=x(randperm(150),:); % initial random category assignment

distmat=squareform(pdist(x,"cityblock")); %get all point pair distances using cityblock measure

assign=zeros(150,1);
for i=1:150
    assign(i)=mod(i-1,K)+1;
end
mds=zeros(1,K); % setting up variable to store the category medoids.  This will just store the index of the medoid point
assignold=assign; % setting up looping variables
assignold(1)=0;

while (sum(assignold~=assign)>0) % continue looping while some points are changing category
    for k=1:K % find current medoids
        sub=distmat(assign==k,:);
        totd=sum(sub);
        [~,mds(k)]=min(totd);
    end
    assignold=assign; % store current categories for comparison
    for i=1:150 % update categories by assigning to closest medoid
        diff=distmat(i,mds);
        [~,ind]=min(diff);
        assign(i)=ind;
    end
end

x1p=x(assign==1,:); % define K-means category sets
x2p=x(assign==2,:);
x3p=x(assign==3,:);
% end K-medoids.  Code below is for assessment and visualization only

ca1=sum(ismember(xa,x1p,'rows')); % assign each of these to the closest of the true categories for assessment purposes
xap=x1p;
camax=ca1;
ca2=sum(ismember(xa,x2p,'rows'));
if ca2>ca1
    xap=x2p;
    camax=ca2;
end
ca3=sum(ismember(xa,x3p,'rows'));
if ca3>camax
    xap=x3p;
    camax=ca3;
end
cb1=sum(ismember(xb,x1p,'rows'));
xbp=x1p;
cbmax=cb1;
cb2=sum(ismember(xb,x2p,'rows'));
if cb2>cb1
    xbp=x2p;
    cbmax=cb2;
end
cb3=sum(ismember(xb,x3p,'rows'));
if cb3>cbmax
    xbp=x3p;
    cbmax=cb3;
end
cc1=sum(ismember(xc,x1p,'rows'));
xcp=x1p;
ccmax=cc1;
cc2=sum(ismember(xc,x2p,'rows'));
if cc2>cc1
    xcp=x2p;
    ccmax=cc2;
end
cc3=sum(ismember(xc,x3p,'rows'));
if cc3>ccmax
    xcp=x3p;
    ccmax=cc3;
end

figure; % visualize, including plotting the medoids
scatter(xap(:,1),xap(:,2))
hold on;
scatter(xbp(:,1),xbp(:,2))
scatter(xcp(:,1),xcp(:,2))
scatter(x(mds,1),x(mds,2),[],'MarkerFaceColor',[0 0 0])
hold off