Probabilistic Clustering

The procedure of the previous lesson wrote a point into a cluster, placed it on the border, or left it outside. All three cases were definite, and only the result stood in the ledger of decisions. A subscriber who falls right between two dense regions and a subscriber who stands at the heart of a cluster were recorded in the same way; how hard the decision was to make was written nowhere. The same held for centroid-based and hierarchical partitioning: the output was a cluster number, not a degree of confidence.

This lesson turns the assignment into a number. Every point is attached to every cluster with a responsibility, the shares of a point sum to one, and the closeness of the shares directly measures the undecidedness of that point. For this, the clusters are no longer written as lists of centroids but as components with their own center and their own variability; the dataset as a whole is a mixture model of these components. There are two things to measure: how much the model beats a random assignment, and how much of the information carried by the soft assignment is lost when one steps down to hard assignment.

• CL37. The dataset is the training portion of 756 rows of the fictional subscriber table from the Introduction to Machine Learning, Data Preparation and Feature Engineering, and Supervised Learning courses, seed 20260218. The label column stays aside in this lesson as well.
• CL38. Three continuous columns are used: mean consumption, volatility and consumption per person. Columns with integer values are left out, because a component can collapse onto a single integer value and reduce the variability in that column to zero.
• CL39. Every component is defined by a weight, one center per column and one variability per column. Dependence between columns is not modeled; this is a model decision and it requires the components to be axis-aligned.
• CL40. A lower bound of 0.02 is placed on the variability. Without a lower bound a component can collapse onto a single point, the density at that point grows without limit, and the fit score loses its meaning. The bound is one fiftieth of unit variability on the scaled columns.
• CL41. The fit score is the per-point mean of the natural logarithm of the density that the model gives to a point. Higher is better, and the score is comparable only on the same set of columns.
• CL42. The baseline is random assignment: the points are distributed randomly into three groups and the parameters are computed from that distribution. This baseline is chosen because it uses the same model form and the same number of parameters; the difference between them comes only from the assignment, not from the model family.
• CL43. The number of clusters is three and is a hyperparameter; the justification for the choice and how arbitrary the choice is are measured in the closing lesson of the topic.
• CL44. The structureless reference is built, as in the previous lesson, by within-column shuffling; the tool belongs to the Statistical Analysis and Exploratory Data Analysis course and the testing theory is not repeated.

A Mixture of Components

The model is this: every point in the dataset is taken to have come from one of three components, but which one it came from is not recorded. Every component $c$ carries a weight $a_c$ (the share of the points that come from that component) and keeps, for each column, a center $m_{cd}$ and a variability $v_{cd}$. The density that the component gives to a point fades quickly with distance from the center:

$$f_c(x) = \prod_{d} \frac{1}{\sqrt{2\pi v_{cd}}}\, \exp\left(-\frac{(x_d - m_{cd})^2}{2\,v_{cd}}\right)$$

The share by which a point is attached to a component is the ratio of the contribution of that component to the total contribution:

$$r_{ic} = \frac{a_c\, f_c(x_i)}{\sum_j a_j\, f_j(x_i)}$$

Since the parameters are unknown the shares cannot be computed, and since the shares are unknown the parameters cannot be computed. Iterative estimation breaks this loop: it temporarily fixes one side, computes the other, then swaps. Every round raises the fit score or leaves it where it is; where it arrives depends on the start.

# mixture.py — MODEL. The FICTIONAL subscriber table of the Introduction to Machine Learning, Data Preparation and Feature Engineering and Supervised Learning courses is rebuilt with the same seed;
# the label column stays aside in this course as well.
import math
import statistics

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [(0.28, 21), (0.22, 17), (0.18, 26), (0.14, 14), (0.18, 23)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


SUBSCRIBERS, DATA = [], []
for i in range(RAW):
    r = rng(SEED + i)
    region = REGION[pick(r(), [x[0] for x in REGION])]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() >= 0.046:
        SUBSCRIBERS.append((10001 + i, region[1], household, rating))
for sub_id, base, household, rating in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + sub_id), []
    for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            base * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - d * 0.05) * 100 + 0.5) / 100)
        r()
    if v:
        mean_value = round(sum(v) / len(v), 2)
        DATA.append({"mean_consumption": mean_value, "volatility": round(max(v) - min(v), 2),
                     "per_person": round(mean_value / household, 3)})

CHOSEN = ["mean_consumption", "volatility", "per_person"]


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


TRAIN = shuffle(DATA, SEED + 90000)[:756]
SCALE = {a: (statistics.fmean(x[a] for x in TRAIN),
             statistics.pstdev([x[a] for x in TRAIN])) for a in CHOSEN}
P = [[(x[a] - SCALE[a][0]) / SCALE[a][1] for a in CHOSEN] for x in TRAIN]
K, FLOOR = 3, 0.02                               # number of clusters and variability floor


def parameters(P, member, k):                    # MODEL: weight, center, variability
    A, M, V = [], [], []
    for c in range(k):
        w = [u[c] for u in member]
        t = sum(w) + 1e-9
        m = [sum(w[i] * P[i][d] for i in range(len(P))) / t for d in range(len(P[0]))]
        V.append([max(sum(w[i] * (P[i][d] - m[d]) ** 2 for i in range(len(P))) / t, FLOOR)
                  for d in range(len(P[0]))])
        A.append(t / len(P))
        M.append(m)
    return A, M, V


def responsibility(P, A, M, V):                  # soft assignment and fit score
    R, total = [], 0.0
    for x in P:
        g = [math.log(A[c]) + sum(-0.5 * (math.log(2 * math.pi * V[c][d])
                                          + (x[d] - M[c][d]) ** 2 / V[c][d])
                                  for d in range(len(x))) for c in range(len(A))]
        b = max(g)
        s = sum(math.exp(u - b) for u in g)
        total += b + math.log(s)
        R.append([math.exp(u - b) / s for u in g])
    return R, total / len(P)


def random_membership(n, k, seed):               # BASELINE: random assignment
    r, U = rng(seed), []
    for _ in range(n):
        c = int(r() * k)
        U.append([1.0 if j == c else 0.0 for j in range(k)])
    return U


def estimate(P, k, seed, rounds=80):             # MODEL: iterative estimation
    U = random_membership(len(P), k, seed)
    for _ in range(rounds):
        U, score = responsibility(P, *parameters(P, U, k))
    return U, score


def column_shuffle(P, seed):
    Y = [list(x) for x in P]
    for s in range(len(P[0])):
        for i, j in enumerate(shuffle(list(range(len(P))), seed + s * 17)):
            Y[i][s] = P[j][s]
    return Y


def baseline_score(P, k, seed):
    return responsibility(P, *parameters(P, random_membership(len(P), k, seed), k))[1]


U, SCORE = estimate(P, K, SEED + 120000)
BASELINE = baseline_score(P, K, SEED + 120000)
SHUFFLED = column_shuffle(P, SEED + 610000)      # FICTIONAL: within-column reshuffling
_, SHUF_SCORE = estimate(SHUFFLED, K, SEED + 120000)
SHUF_BASELINE = baseline_score(SHUFFLED, K, SEED + 120000)
print(f"points {len(P)}, columns {len(CHOSEN)}, components {K}")
print(f"{'dataset':<16}{'baseline':>15}{'mixture model':>16}{'difference':>12}")
print(f"{'data':<16}{BASELINE:>15.4f}{SCORE:>16.4f}{SCORE - BASELINE:>+12.4f}")
print(f"{'shuffled':<16}{SHUF_BASELINE:>15.4f}{SHUF_SCORE:>16.4f}{SHUF_SCORE - SHUF_BASELINE:>+12.4f}")
print(f"share from the data {(SCORE - BASELINE) - (SHUF_SCORE - SHUF_BASELINE):+.4f}")

points 756, columns 3, components 3
dataset                baseline   mixture model  difference
data                    -4.2424         -3.7755     +0.4669
shuffled                -4.2443         -3.9085     +0.3358
share from the data +0.1311

The random assignment baseline is -4.2424, the mixture model -3.7755, the difference +0.4669. The baseline has been chosen deliberately hard here: random assignment also builds three components, computes the same number of parameters and uses the same density form. Its only lack is that it does not learn which point belongs to which component. In this way the difference comes not from the model family but only from the assignment.

The second row says how much of this difference comes from the data. On the dataset whose columns are shuffled within themselves the same procedure gains +0.3358. That is, about three quarters of the gain stands even when no relationship between columns is left. The reason is clear: three components wrap a skewed, long-tailed distribution axis by axis better than a single component could; this says nothing about the cluster structure in the data, it says something about the flexibility of the model. The share that falls to the structure between columns is +0.1311, about 28% of the total difference. The question of the course finds its answer exactly here: the answer to whether the difference comes from the measure or from the data is often “from both”, and without separating the shares the number cannot be read.

The Uncertainty That Soft Assignment Carries

The responsibility keeps the information that hard assignment has to throw away. If the shares of a point are $0.98 / 0.01 / 0.01$ the assignment is definite; if they are $0.40 / 0.38 / 0.22$ the assignment is a choice, and how fine the choice is stands written in the shares. To reduce the uncertainty to a single number, the sum of the shares weighted by their natural logarithm is used; this number is zero when all the weight is given to a single share and reaches its largest value when the shares are distributed equally.

def uncertainty(r):
    return -sum(p * math.log(p) for p in r if p > 1e-12)


HARD = [max(range(K), key=lambda c: r[c]) for r in U]
uncertain = [i for i, r in enumerate(U) if max(r) < 0.60]
A, _, _ = parameters(P, U, K)
print(f"points with largest responsibility below 0.60: {len(uncertain)}; "
      f"for these points the mean largest share is {sum(max(U[i]) for i in uncertain) / len(uncertain):.4f}")
print(f"points with the top two shares closer than 0.10: {sum(1 for r in U if sorted(r)[-1] - sorted(r)[-2] < 0.10)}"
      f"; largest share 0.99 and above: {sum(1 for r in U if max(r) >= 0.99)}")
print(f"mean uncertainty: soft assignment "
      f"{sum(uncertainty(r) for r in U) / len(U):.4f}, hard assignment 0.0000, "
      f"upper bound {math.log(K):.4f}")
print(f"{'cluster':<8}{'soft mass':>15}{'hard count':>12}{'difference':>12}")
for c in range(K):
    print(f"{c:<8}{A[c] * len(P):>15.1f}{HARD.count(c):>12}"
          f"{HARD.count(c) - A[c] * len(P):>+12.1f}")

points with largest responsibility below 0.60: 74; for these points the mean largest share is 0.5418
points with the top two shares closer than 0.10: 28; largest share 0.99 and above: 94
mean uncertainty: soft assignment 0.3797, hard assignment 0.0000, upper bound 1.0986
cluster       soft mass  hard count  difference
0                 217.8         189       -28.8
1                 200.3         213       +12.7
2                 337.9         354       +16.1

The threshold of 0.60 is chosen because an equal distribution over three components comes to 0.3333, and a largest share below 0.60 means that the majority of the weight of the point is not gathered in a single cluster. 74 of the 756 points are in this situation, and for these points the largest share is 0.5418 on average. For twenty-eight points the top two shares differ from each other by less than 0.10 — for these points the cluster chosen by hard assignment and the cluster in second place are almost equally justified. At the other end, 94 points are attached to a single cluster with a share of 0.99 and above.

The mean uncertainty is 0.3797 in soft assignment, against an upper bound of 1.0986. On moving to hard assignment this number falls to zero by definition: every point is written to a single cluster with full weight, and the information that 74 points are undecided does not remain in the table. The concrete counterpart of the loss is in the cluster sizes. In soft assignment the mass of the first cluster comes to 217.8 points, while hard assignment writes 189 points to that cluster; the deviation of 28.8 points comes from all the points on the border being rounded to one side. Every decision that rests on the relative size of the clusters inherits this rounding.

Testing the Shares

Up to here it has been assumed that the responsibility is meaningful. This is a testable claim: if the share really measures how hard the decision is, undecided points should change cluster when the conditions change a little, and certain points should not. The condition to change is ready — iterative estimation is sensitive to the start, and the second lesson of the topic measured this for centroid-based partitioning. The model run from a second start is compared with the first after the clusters are matched through their centers.

U2, SCORE2 = estimate(P, K, SEED + 555000)
_, M, _ = parameters(P, U, K)
_, M2, _ = parameters(P, U2, K)
MATCH = [min(range(K), key=lambda c: math.dist(M2[j], M[c])) for j in range(K)]
HARD2 = [MATCH[max(range(K), key=lambda c: r[c])] for r in U2]
changed = [i for i in range(len(P)) if HARD[i] != HARD2[i]]
certain = [i for i, r in enumerate(U) if max(r) >= 0.90]
changed_uncertain = sum(1 for i in changed if max(U[i]) < 0.60)
changed_certain = sum(1 for i in changed if max(U[i]) >= 0.90)
print(f"fit score of the second start {SCORE2:.4f} (first {SCORE:.4f}); "
      f"cluster matching {MATCH}")
print(f"points that changed cluster {len(changed)} / {len(P)}")
print(f"change among uncertain points {changed_uncertain}/{len(uncertain)} = {changed_uncertain / len(uncertain):.4f}")
print(f"change among certain points (share >= 0.90) {changed_certain}/{len(certain)} = {changed_certain / len(certain):.4f}")

fit score of the second start -3.7827 (first -3.7755); cluster matching [1, 0, 2]
points that changed cluster 125 / 756
change among uncertain points 32/74 = 0.4324
change among certain points (share >= 0.90) 29/403 = 0.0720

The two runs stand in different places: the fit scores are -3.7755 and -3.7827. The difference between them is small, but the difference between the partitions is not small — 125 of the 756 points change cluster. This shows that looking only at the fit score and saying “the model converged” is not enough.

The real number is in the disagreement. 32 of the 74 undecided points, that is 43.24%, change cluster. Of the 403 points with a share of 0.90 and above, only 29, that is 7.20%. The ratio is more than sixfold. The responsibility correctly estimates something outside the model it was computed from: how fragile the assignment of a point is. Hard assignment does not carry this information — a cluster number shows a number that will change with a 43% probability and a number that will change with a 7% probability in the same way.

Summary

• The mixture model writes the clusters as components defined by a center and a variability; every point is attached to every component with a responsibility, and the shares sum to one.
• The random assignment baseline is -4.2424, the fit score of the model is -3.7755, the difference is +0.4669; because the baseline uses the same model form and the same number of parameters, the difference comes only from the assignment.
• The same difference comes out as +0.3358 on the dataset shuffled within its columns; the share that falls to the structure between columns is +0.1311, that is, about 28% of the total difference. The rest is the flexibility of the model.
• For 74 of the 756 points the largest share is below 0.60 (mean 0.5418), and for 28 points the top two shares are closer than 0.10. On moving to hard assignment the mean uncertainty falls from 0.3797 to zero, and the count of one cluster deviates by 28.8 points from its soft mass.
• The responsibility also measures something outside the model: at a second start, 43.24% of the undecided points and 7.20% of the certain points change cluster.

Next Step

In this lesson the number of clusters was taken as three and this choice was not justified. Similar choices were made in the previous three lessons as well: how many clusters, which cut height, which neighborhood radius. Each time it was shown that the choice moved the result, but how to make the choice itself was postponed. The closing lesson of the topic takes up this question and shows with numbers that there is no answer: on the same data three different internal validity measures choose three different numbers of clusters, none of them more correct than another, and on a structureless dataset generated entirely at random the same three measures still choose a number of clusters. After that, the label is opened, for the first and only time in this course, for external validity, and it is counted how far the number of clusters that the label points to disagrees with those chosen by the internal measures.