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
Spatial Data Structures - KD-trees for clustering
advanced_filters_comparison - Gaussian sum filters