Gaussian Mixtures and Clustering

This example demonstrates Gaussian mixture operations and clustering algorithms.

Overview

Gaussian mixtures and clustering are essential for:

  • Multi-hypothesis tracking: Representing multiple target hypotheses

  • Mixture reduction: Limiting computational complexity

  • Data clustering: Grouping measurements or tracks

Clustering Algorithms

K-Means
  • Partitions data into k clusters

  • Minimizes within-cluster variance

  • Fast and scalable

DBSCAN
  • Density-based clustering

  • Handles arbitrary cluster shapes

  • Identifies outliers automatically

Hierarchical
  • Builds cluster tree (dendrogram)

  • Multiple linkage options

  • Flexible number of clusters

Gaussian Mixture Operations

  • Moment matching: Reducing mixture to single Gaussian

  • Runnalls’ algorithm: Optimal mixture reduction

  • West’s algorithm: Alternative reduction method

  • Splitting/merging: Dynamic mixture management

Code Highlights

The example demonstrates:

  • K-means clustering with kmeans()

  • DBSCAN with dbscan()

  • Hierarchical clustering with hierarchical_cluster()

  • Gaussian mixture reduction with reduce_mixture_runnalls()

Source Code

  1"""
  2Gaussian Mixtures Example
  3=========================
  4
  5This example demonstrates Gaussian mixture operations and clustering
  6algorithms in PyTCL:
  7
  8Gaussian Mixture Operations:
  9- Component representation and manipulation
 10- Moment matching (computing mean and covariance)
 11- Mixture merging and reduction
 12- Runnalls' and West's reduction algorithms
 13
 14Clustering Algorithms:
 15- K-means with K-means++ initialization
 16- DBSCAN (density-based clustering)
 17- Hierarchical/agglomerative clustering
 18- Elbow method for K selection
 19
 20These algorithms are essential for multi-target tracking (PHD filters),
 21hypothesis reduction in MHT, and general density estimation.
 22"""
 23
 24from pathlib import Path
 25
 26import numpy as np
 27import plotly.graph_objects as go
 28
 29# Output directory for generated plots
 30OUTPUT_DIR = Path(__file__).parent.parent / "docs" / "_static" / "images" / "examples"
 31OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
 32from plotly.subplots import make_subplots
 33
 34# Global flag to control plotting
 35SHOW_PLOTS = True
 36
 37
 38def create_ellipse_trace(mean, cov, color="blue", opacity=0.3, n_std=2, name=None):
 39    """Create a plotly trace for an ellipse representing a 2D Gaussian covariance."""
 40    eigvals, eigvecs = np.linalg.eigh(cov)
 41    # Angle in radians
 42    angle = np.arctan2(eigvecs[1, 0], eigvecs[0, 0])
 43
 44    # Semi-axes lengths
 45    a = n_std * np.sqrt(eigvals[0])
 46    b = n_std * np.sqrt(eigvals[1])
 47
 48    # Parametric ellipse
 49    t = np.linspace(0, 2 * np.pi, 100)
 50    x = a * np.cos(t)
 51    y = b * np.sin(t)
 52
 53    # Rotate
 54    x_rot = x * np.cos(angle) - y * np.sin(angle) + mean[0]
 55    y_rot = x * np.sin(angle) + y * np.cos(angle) + mean[1]
 56
 57    return go.Scatter(
 58        x=x_rot,
 59        y=y_rot,
 60        mode="lines",
 61        fill="toself",
 62        fillcolor=color,
 63        opacity=opacity,
 64        line=dict(color=color, width=2),
 65        name=name,
 66        showlegend=name is not None,
 67    )
 68
 69
 70from pytcl.clustering import (  # Gaussian mixture operations; K-means; DBSCAN; Hierarchical clustering
 71    DBSCANResult,
 72    GaussianComponent,
 73    GaussianMixture,
 74    HierarchicalResult,
 75    KMeansResult,
 76    agglomerative_clustering,
 77    compute_distance_matrix,
 78    cut_dendrogram,
 79    dbscan,
 80    dbscan_predict,
 81    kmeans,
 82    kmeans_elbow,
 83    kmeans_plusplus_init,
 84    merge_gaussians,
 85    moment_match,
 86    prune_mixture,
 87    reduce_mixture_runnalls,
 88    reduce_mixture_west,
 89    runnalls_merge_cost,
 90    west_merge_cost,
 91)
 92
 93
 94def demo_gaussian_components():
 95    """Demonstrate Gaussian component operations."""
 96    print("=" * 70)
 97    print("Gaussian Component Operations Demo")
 98    print("=" * 70)
 99
100    # Create individual Gaussian components
101    comp1 = GaussianComponent(
102        weight=0.4,
103        mean=np.array([0.0, 0.0]),
104        covariance=np.array([[1.0, 0.0], [0.0, 1.0]]),
105    )
106
107    comp2 = GaussianComponent(
108        weight=0.6,
109        mean=np.array([3.0, 3.0]),
110        covariance=np.array([[2.0, 0.5], [0.5, 2.0]]),
111    )
112
113    print("\nComponent 1:")
114    print(f"  Weight: {comp1.weight}")
115    print(f"  Mean: {comp1.mean}")
116    print(f"  Covariance diagonal: {np.diag(comp1.covariance)}")
117
118    print("\nComponent 2:")
119    print(f"  Weight: {comp2.weight}")
120    print(f"  Mean: {comp2.mean}")
121    print(f"  Covariance diagonal: {np.diag(comp2.covariance)}")
122
123    # Create a mixture
124    mixture = GaussianMixture([comp1, comp2])
125    print(f"\nMixture has {len(mixture)} components")
126    print(f"Total weight: {sum(c.weight for c in mixture.components)}")
127
128
129def demo_moment_matching():
130    """Demonstrate moment matching for mixture approximation."""
131    print("\n" + "=" * 70)
132    print("Moment Matching Demo")
133    print("=" * 70)
134
135    # Create a mixture of 3 components
136    weights = np.array([0.3, 0.5, 0.2])
137    means = [np.array([0.0, 0.0]), np.array([2.0, 1.0]), np.array([1.0, 3.0])]
138    covariances = [np.eye(2) * 0.5, np.eye(2) * 0.8, np.eye(2) * 0.3]
139
140    print("\nOriginal mixture (3 components):")
141    for i, (w, m) in enumerate(zip(weights, means)):
142        print(f"  Component {i+1}: weight={w:.1f}, mean={m}")
143
144    # Moment match to single Gaussian - takes (weights, means, covariances)
145    mean, cov = moment_match(weights, means, covariances)
146
147    print("\nMoment-matched single Gaussian:")
148    print(f"  Mean: ({mean[0]:.3f}, {mean[1]:.3f})")
149    print(f"  Covariance:\n{cov}")
150
151    # The mean should be the weighted average
152    weighted_mean = sum(w * m for w, m in zip(weights, means))
153    print(
154        f"\nVerification - weighted mean: ({weighted_mean[0]:.3f}, "
155        f"{weighted_mean[1]:.3f})"
156    )
157
158
159def demo_mixture_merging():
160    """Demonstrate merging Gaussian components."""
161    print("\n" + "=" * 70)
162    print("Mixture Merging Demo")
163    print("=" * 70)
164
165    # Two nearby components
166    comp1 = GaussianComponent(0.4, np.array([0.0, 0.0]), np.eye(2) * 1.0)
167    comp2 = GaussianComponent(0.3, np.array([0.5, 0.5]), np.eye(2) * 1.0)
168
169    # Two far apart components
170    comp3 = GaussianComponent(0.3, np.array([5.0, 5.0]), np.eye(2) * 1.0)
171
172    print("\nThree components:")
173    print(f"  1: mean={comp1.mean}, weight={comp1.weight}")
174    print(f"  2: mean={comp2.mean}, weight={comp2.weight}")
175    print(f"  3: mean={comp3.mean}, weight={comp3.weight}")
176
177    # Merge costs
178    cost_12 = runnalls_merge_cost(comp1, comp2)
179    cost_13 = runnalls_merge_cost(comp1, comp3)
180    cost_23 = runnalls_merge_cost(comp2, comp3)
181
182    print("\nRunnalls merge costs (lower = better candidates for merging):")
183    print(f"  Cost(1,2): {cost_12:.4f}")
184    print(f"  Cost(1,3): {cost_13:.4f}")
185    print(f"  Cost(2,3): {cost_23:.4f}")
186
187    # Merge the closest pair
188    merge_result = merge_gaussians(comp1, comp2)
189    merged = merge_result.component
190
191    print("\nMerged component (1+2):")
192    print(f"  Weight: {merged.weight:.2f}")
193    print(f"  Mean: ({merged.mean[0]:.3f}, {merged.mean[1]:.3f})")
194    print(f"  Merge cost: {merge_result.cost:.4f}")
195
196
197def demo_mixture_reduction():
198    """Demonstrate mixture reduction algorithms."""
199    print("\n" + "=" * 70)
200    print("Mixture Reduction Demo")
201    print("=" * 70)
202
203    np.random.seed(42)
204
205    # Create a large mixture (e.g., from PHD filter output)
206    n_components = 20
207    components = []
208
209    # Cluster 1: around (0, 0)
210    for _ in range(8):
211        mean = np.array([0, 0]) + np.random.randn(2) * 0.5
212        cov = np.eye(2) * (0.3 + np.random.rand() * 0.2)
213        weight = 0.3 + np.random.rand() * 0.4
214        components.append(GaussianComponent(weight, mean, cov))
215
216    # Cluster 2: around (5, 3)
217    for _ in range(7):
218        mean = np.array([5, 3]) + np.random.randn(2) * 0.5
219        cov = np.eye(2) * (0.3 + np.random.rand() * 0.2)
220        weight = 0.2 + np.random.rand() * 0.3
221        components.append(GaussianComponent(weight, mean, cov))
222
223    # Scattered components
224    for _ in range(5):
225        mean = np.random.rand(2) * 10
226        cov = np.eye(2) * 0.5
227        weight = 0.05 + np.random.rand() * 0.1
228        components.append(GaussianComponent(weight, mean, cov))
229
230    print(f"\nOriginal mixture: {len(components)} components")
231    print(f"Total weight: {sum(c.weight for c in components):.2f}")
232
233    # Prune low-weight components
234    pruned = prune_mixture(components, weight_threshold=0.1)
235    print(f"\nAfter pruning (weight_threshold=0.1): {len(pruned)} components")
236
237    # Reduce using Runnalls' algorithm
238    n_target = 5
239    result_runnalls = reduce_mixture_runnalls(components, n_target)
240    print(f"\nRunnalls reduction to {n_target} components:")
241    print(f"  Final components: {result_runnalls.n_reduced}")
242    print(f"  Total merge cost: {result_runnalls.total_cost:.4f}")
243    for i, c in enumerate(result_runnalls.components):
244        print(
245            f"    {i+1}: weight={c.weight:.3f}, "
246            f"mean=({c.mean[0]:.2f}, {c.mean[1]:.2f})"
247        )
248
249    # Reduce using West's algorithm
250    result_west = reduce_mixture_west(components, n_target)
251    print(f"\nWest reduction to {n_target} components:")
252    print(f"  Final components: {result_west.n_reduced}")
253    print(f"  Total merge cost: {result_west.total_cost:.4f}")
254
255
256def demo_kmeans():
257    """Demonstrate K-means clustering."""
258    print("\n" + "=" * 70)
259    print("K-Means Clustering Demo")
260    print("=" * 70)
261
262    np.random.seed(42)
263
264    # Generate clustered data
265    n_per_cluster = 50
266    centers_true = np.array(
267        [
268            [0, 0],
269            [5, 5],
270            [0, 5],
271        ]
272    )
273
274    data = []
275    for center in centers_true:
276        cluster = center + np.random.randn(n_per_cluster, 2) * 0.8
277        data.append(cluster)
278    data = np.vstack(data)
279
280    print(f"\nGenerated {len(data)} points in 3 clusters")
281    print(f"True centers:\n{centers_true}")
282
283    # K-means clustering
284    result = kmeans(data, n_clusters=3, max_iter=100)
285
286    print(f"\nK-means result:")
287    print(f"  Iterations: {result.n_iter}")
288    print(f"  Inertia (within-cluster sum of squares): {result.inertia:.2f}")
289    print(f"  Found centers:\n{result.centers}")
290
291    # Cluster sizes
292    unique, counts = np.unique(result.labels, return_counts=True)
293    print("\n  Cluster sizes:", dict(zip(unique, counts)))
294
295    # Plot K-means result
296    if SHOW_PLOTS:
297        fig = make_subplots(
298            rows=1,
299            cols=2,
300            subplot_titles=["Ground Truth Clusters", "K-means Clustering Result"],
301        )
302
303        # True clusters
304        colors = ["blue", "green", "orange"]
305        for i, center in enumerate(centers_true):
306            mask = np.arange(len(data)) // n_per_cluster == i
307            fig.add_trace(
308                go.Scatter(
309                    x=data[mask, 0],
310                    y=data[mask, 1],
311                    mode="markers",
312                    marker=dict(color=colors[i], size=6, opacity=0.6),
313                    name=f"True cluster {i}",
314                ),
315                row=1,
316                col=1,
317            )
318
319        fig.add_trace(
320            go.Scatter(
321                x=centers_true[:, 0],
322                y=centers_true[:, 1],
323                mode="markers",
324                marker=dict(color="black", size=15, symbol="x", line=dict(width=3)),
325                name="True centers",
326            ),
327            row=1,
328            col=1,
329        )
330
331        # K-means result
332        colors_km = ["red", "green", "blue"]
333        for i in range(3):
334            mask = result.labels == i
335            fig.add_trace(
336                go.Scatter(
337                    x=data[mask, 0],
338                    y=data[mask, 1],
339                    mode="markers",
340                    marker=dict(color=colors_km[i], size=6, opacity=0.6),
341                    name=f"Cluster {i}",
342                    showlegend=True,
343                ),
344                row=1,
345                col=2,
346            )
347
348        fig.add_trace(
349            go.Scatter(
350                x=result.centers[:, 0],
351                y=result.centers[:, 1],
352                mode="markers",
353                marker=dict(color="black", size=15, symbol="x", line=dict(width=3)),
354                name="K-means centers",
355            ),
356            row=1,
357            col=2,
358        )
359
360        fig.update_xaxes(title_text="x")
361        fig.update_yaxes(title_text="y")
362        fig.update_layout(height=500, width=1000, showlegend=True)
363        fig.write_html(str(OUTPUT_DIR / "gaussian_kmeans.html"))
364        print("\n  [Plot saved to gaussian_kmeans.html]")
365
366
367def demo_kmeans_plusplus():
368    """Demonstrate K-means++ initialization."""
369    print("\n" + "=" * 70)
370    print("K-Means++ Initialization Demo")
371    print("=" * 70)
372
373    np.random.seed(42)
374
375    # Generate data with 4 well-separated clusters
376    data = np.vstack(
377        [
378            np.random.randn(30, 2) + [0, 0],
379            np.random.randn(30, 2) + [10, 0],
380            np.random.randn(30, 2) + [0, 10],
381            np.random.randn(30, 2) + [10, 10],
382        ]
383    )
384
385    n_clusters = 4
386
387    # Random initialization
388    random_centers = data[np.random.choice(len(data), n_clusters, replace=False)]
389
390    # K-means++ initialization
391    plusplus_centers = kmeans_plusplus_init(data, n_clusters)
392
393    print(f"\nComparing initialization methods for n_clusters={n_clusters}:")
394
395    # Run K-means with each initialization
396    # For random, run multiple times (use init='random' for random init)
397    results_random = []
398    for _ in range(10):
399        idx = np.random.choice(len(data), n_clusters, replace=False)
400        result = kmeans(data, n_clusters=n_clusters, init=data[idx], n_init=1)
401        results_random.append(result.inertia)
402
403    result_plusplus = kmeans(
404        data, n_clusters=n_clusters, init=plusplus_centers, n_init=1
405    )
406
407    print(f"\n  Random initialization (10 runs):")
408    print(f"    Mean inertia: {np.mean(results_random):.2f}")
409    print(f"    Std inertia: {np.std(results_random):.2f}")
410    print(f"    Best inertia: {np.min(results_random):.2f}")
411
412    print(f"\n  K-means++ initialization:")
413    print(f"    Inertia: {result_plusplus.inertia:.2f}")
414
415    print("\nNote: K-means++ typically provides better, more consistent results.")
416
417
418def demo_elbow_method():
419    """Demonstrate elbow method for K selection."""
420    print("\n" + "=" * 70)
421    print("Elbow Method Demo")
422    print("=" * 70)
423
424    np.random.seed(42)
425
426    # Generate data with 3 true clusters
427    data = np.vstack(
428        [
429            np.random.randn(40, 2) + [0, 0],
430            np.random.randn(40, 2) + [4, 4],
431            np.random.randn(40, 2) + [8, 0],
432        ]
433    )
434
435    print(f"\nData: 120 points from 3 true clusters")
436    print("\nInertia for different K values:")
437    print("-" * 40)
438
439    elbow_result = kmeans_elbow(data, k_range=range(1, 8))
440    k_values = elbow_result["k_values"]
441    inertias = elbow_result["inertias"]
442
443    for k, inertia in zip(k_values, inertias):
444        bar = "#" * int(inertia / max(inertias) * 30)
445        print(f"  K={k}: {inertia:>8.1f} {bar}")
446
447    print("\nThe 'elbow' should appear around K=3")
448    print("(where adding more clusters gives diminishing returns)")
449
450    # Plot elbow method
451    if SHOW_PLOTS:
452        fig = go.Figure()
453
454        fig.add_trace(
455            go.Scatter(
456                x=list(k_values),
457                y=inertias,
458                mode="lines+markers",
459                line=dict(color="blue", width=2),
460                marker=dict(size=10),
461                name="Inertia",
462            )
463        )
464
465        fig.add_vline(
466            x=3, line_dash="dash", line_color="red", annotation_text="True K=3"
467        )
468
469        fig.update_layout(
470            title="Elbow Method for K Selection",
471            xaxis_title="Number of Clusters (K)",
472            yaxis_title="Inertia (Within-cluster Sum of Squares)",
473            height=500,
474            width=700,
475            showlegend=True,
476        )
477        fig.write_html(str(OUTPUT_DIR / "gaussian_elbow.html"))
478        print("\n  [Plot saved to gaussian_elbow.html]")
479
480
481def demo_dbscan():
482    """Demonstrate DBSCAN clustering."""
483    print("\n" + "=" * 70)
484    print("DBSCAN Clustering Demo")
485    print("=" * 70)
486
487    np.random.seed(42)
488
489    # Generate data: two dense clusters + noise
490    cluster1 = np.random.randn(50, 2) * 0.5 + [0, 0]
491    cluster2 = np.random.randn(50, 2) * 0.5 + [4, 4]
492    noise = np.random.uniform(-2, 8, (20, 2))  # Scattered noise points
493
494    data = np.vstack([cluster1, cluster2, noise])
495
496    print(f"\nData: 100 cluster points + 20 noise points")
497
498    # DBSCAN clustering
499    result = dbscan(data, eps=0.8, min_samples=5)
500
501    print(f"\nDBSCAN result (eps=0.8, min_samples=5):")
502    print(f"  Clusters found: {result.n_clusters}")
503    print(f"  Core sample indices: {len(result.core_sample_indices)}")
504    print(f"  Noise points: {result.n_noise}")
505
506    # Cluster sizes
507    unique_labels = np.unique(result.labels)
508    for label in unique_labels:
509        count = np.sum(result.labels == label)
510        if label == -1:
511            print(f"  Noise points: {count}")
512        else:
513            print(f"  Cluster {label}: {count} points")
514
515    print("\nNote: DBSCAN identifies noise points (label=-1)")
516    print("and doesn't require specifying the number of clusters.")
517
518    # Plot DBSCAN result
519    if SHOW_PLOTS:
520        fig = make_subplots(
521            rows=1,
522            cols=2,
523            subplot_titles=[
524                "Ground Truth",
525                f"DBSCAN Result ({result.n_clusters} clusters)",
526            ],
527        )
528
529        # Ground truth
530        fig.add_trace(
531            go.Scatter(
532                x=cluster1[:, 0],
533                y=cluster1[:, 1],
534                mode="markers",
535                marker=dict(color="blue", size=8, opacity=0.6),
536                name="Cluster 1",
537            ),
538            row=1,
539            col=1,
540        )
541
542        fig.add_trace(
543            go.Scatter(
544                x=cluster2[:, 0],
545                y=cluster2[:, 1],
546                mode="markers",
547                marker=dict(color="green", size=8, opacity=0.6),
548                name="Cluster 2",
549            ),
550            row=1,
551            col=1,
552        )
553
554        fig.add_trace(
555            go.Scatter(
556                x=noise[:, 0],
557                y=noise[:, 1],
558                mode="markers",
559                marker=dict(color="red", size=8, opacity=0.6),
560                name="Noise",
561            ),
562            row=1,
563            col=1,
564        )
565
566        # DBSCAN result
567        colors = ["blue", "green", "purple", "orange"]
568        for label in unique_labels:
569            mask = result.labels == label
570            if label == -1:
571                fig.add_trace(
572                    go.Scatter(
573                        x=data[mask, 0],
574                        y=data[mask, 1],
575                        mode="markers",
576                        marker=dict(color="red", size=8, opacity=0.6, symbol="x"),
577                        name="Noise",
578                        showlegend=True,
579                    ),
580                    row=1,
581                    col=2,
582                )
583            else:
584                fig.add_trace(
585                    go.Scatter(
586                        x=data[mask, 0],
587                        y=data[mask, 1],
588                        mode="markers",
589                        marker=dict(
590                            color=colors[label % len(colors)], size=8, opacity=0.6
591                        ),
592                        name=f"Cluster {label}",
593                        showlegend=True,
594                    ),
595                    row=1,
596                    col=2,
597                )
598
599        # Mark core samples
600        core_mask = np.zeros(len(data), dtype=bool)
601        core_mask[result.core_sample_indices] = True
602        fig.add_trace(
603            go.Scatter(
604                x=data[core_mask, 0],
605                y=data[core_mask, 1],
606                mode="markers",
607                marker=dict(
608                    color="rgba(0,0,0,0)", size=15, line=dict(color="black", width=1)
609                ),
610                name="Core samples",
611                showlegend=True,
612            ),
613            row=1,
614            col=2,
615        )
616
617        fig.update_xaxes(title_text="x")
618        fig.update_yaxes(title_text="y")
619        fig.update_layout(height=500, width=1000, showlegend=True)
620        fig.write_html(str(OUTPUT_DIR / "gaussian_dbscan.html"))
621        print("\n  [Plot saved to gaussian_dbscan.html]")
622
623
624def demo_hierarchical():
625    """Demonstrate hierarchical clustering."""
626    print("\n" + "=" * 70)
627    print("Hierarchical Clustering Demo")
628    print("=" * 70)
629
630    np.random.seed(42)
631
632    # Generate small dataset for visualization
633    data = np.array(
634        [
635            [0, 0],
636            [0.5, 0.5],
637            [1, 0],  # Cluster A
638            [5, 5],
639            [5.5, 5.5],
640            [5, 6],  # Cluster B
641            [2.5, 2.5],  # Between clusters
642        ]
643    )
644
645    print(f"\nData points:\n{data}")
646
647    # Compute distance matrix
648    dist_matrix = compute_distance_matrix(data)
649    print(f"\nDistance matrix:\n{np.round(dist_matrix, 2)}")
650
651    # Agglomerative clustering
652    result = agglomerative_clustering(data, linkage="average")
653
654    print("\nHierarchical clustering (average linkage):")
655    print(f"  Labels: {result.labels}")
656    print(f"  Number of clusters: {result.n_clusters}")
657
658    # Cut at different thresholds
659    n_samples = len(data)
660    for threshold in [1.0, 3.0, 5.0]:
661        labels = cut_dendrogram(
662            result.linkage_matrix, n_samples, distance_threshold=threshold
663        )
664        n_clusters = len(set(labels))
665        print(
666            f"  Threshold {threshold:.1f}: {n_clusters} clusters, labels={list(labels)}"
667        )
668
669
670def demo_tracking_application():
671    """Demonstrate mixture reduction in tracking context."""
672    print("\n" + "=" * 70)
673    print("Tracking Application Demo")
674    print("=" * 70)
675
676    np.random.seed(42)
677
678    # Simulated PHD filter output: mixture representing target density
679    # After several updates, the mixture can have many components
680
681    # True targets at these locations
682    true_targets = np.array(
683        [
684            [10.0, 20.0],
685            [30.0, 40.0],
686            [25.0, 25.0],
687        ]
688    )
689
690    print(f"\nTrue target positions: {len(true_targets)} targets")
691    for i, t in enumerate(true_targets):
692        print(f"  Target {i+1}: ({t[0]:.1f}, {t[1]:.1f})")
693
694    # Create mixture with components clustered around true targets
695    # Plus some spurious components (false alarms, etc.)
696    components = []
697
698    # Components near true targets (higher weight)
699    for target in true_targets:
700        for _ in range(4):
701            mean = target + np.random.randn(2) * 1.0
702            cov = np.eye(2) * (0.5 + np.random.rand() * 0.5)
703            weight = 0.6 + np.random.rand() * 0.4
704            components.append(GaussianComponent(weight, mean, cov))
705
706    # Spurious components (lower weight)
707    for _ in range(8):
708        mean = np.random.uniform(5, 45, 2)
709        cov = np.eye(2) * 2.0
710        weight = 0.05 + np.random.rand() * 0.1
711        components.append(GaussianComponent(weight, mean, cov))
712
713    print(f"\nPHD mixture: {len(components)} components")
714    total_weight = sum(c.weight for c in components)
715    print(f"  Total weight (expected target count): {total_weight:.2f}")
716
717    # Reduce to extract target estimates
718    n_expected = int(round(total_weight))
719    reduced = reduce_mixture_runnalls(components, n_expected)
720
721    print(f"\nAfter reduction to {n_expected} components:")
722    for i, c in enumerate(reduced.components):
723        # Find closest true target
724        dists = [np.linalg.norm(c.mean - t) for t in true_targets]
725        closest = np.argmin(dists)
726        error = min(dists)
727        print(
728            f"  Estimate {i+1}: ({c.mean[0]:.1f}, {c.mean[1]:.1f}), "
729            f"weight={c.weight:.2f}, error to target {closest+1}={error:.2f}"
730        )
731
732
733def main():
734    """Run all demonstrations."""
735    print("\n" + "#" * 70)
736    print("# PyTCL Gaussian Mixtures and Clustering Example")
737    print("#" * 70)
738
739    # Gaussian mixture operations
740    demo_gaussian_components()
741    demo_moment_matching()
742    demo_mixture_merging()
743    demo_mixture_reduction()
744
745    # Clustering algorithms
746    demo_kmeans()
747    demo_kmeans_plusplus()
748    demo_elbow_method()
749    demo_dbscan()
750    demo_hierarchical()
751
752    # Application
753    demo_tracking_application()
754
755    print("\n" + "=" * 70)
756    print("Example complete!")
757    if SHOW_PLOTS:
758        print(
759            "Plots saved: gaussian_kmeans.html, gaussian_elbow.html, gaussian_dbscan.html"
760        )
761    print("=" * 70)
762
763
764if __name__ == "__main__":
765    main()

Running the Example

python examples/gaussian_mixtures.py

See Also