Week 10: Unsupervised Learning — Clustering and PCA

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs, load_iris
from sklearn.preprocessing import StandardScaler

print("="*50)
print("UNSUPERVISED LEARNING FUNDAMENTALS")
print("="*50)

print("""
SUPERVISED vs UNSUPERVISED:

SUPERVISED (Week 9):
- We have labeled training data (X, y)
- Learn mapping from X to y
- Examples: emails→spam/not-spam, images→cat/dog
- Easy to evaluate: compare predictions to labels

UNSUPERVISED (Week 10):
- Only have data X, no labels y
- Discover structure and patterns
- Examples: group customers, compress images
- Harder to evaluate: no ground truth

TYPES OF UNSUPERVISED LEARNING:

1. CLUSTERING
- Group similar items
- Questions: How many groups? What's similar?
- Algorithms: K-means, Hierarchical, DBSCAN, GMM

2. DIMENSIONALITY REDUCTION
- Reduce number of features
- Questions: Which features matter? How to visualize?
- Algorithms: PCA, t-SNE, Autoencoders

3. ANOMALY DETECTION
- Find unusual items
- Questions: What's "normal"? What's outlier?
- Algorithms: Isolation Forest, One-class SVM

4. ASSOCIATION RULES
- Find relationships
- Questions: What items bought together?
- Algorithms: Apriori, Eclat

THIS WEEK: CLUSTERING AND DIMENSIONALITY REDUCTION
""")

# Create sample data
print("\nEXAMPLE: CUSTOMER SEGMENTATION")
print("-"*50)

# Generate synthetic customer data
np.random.seed(42)
n_samples = 300

# Create clusters
cluster1 = np.random.normal([2, 2], 0.5, (100, 2))
cluster2 = np.random.normal([8, 1], 0.5, (100, 2))
cluster3 = np.random.normal([5, 8], 0.5, (100, 2))

X = np.vstack([cluster1, cluster2, cluster3])
y_true = np.hstack([np.zeros(100), np.ones(100), np.ones(100)*2])

print(f"Number of customers: {len(X)}")
print(f"Features: 2 (e.g., spending, frequency)")
print(f"True underlying clusters: 3")
print(f"First few samples:")
print(X[:5])

# Visualize
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# True labels (for reference - we won't use in unsupervised learning)
scatter1 = axes[0].scatter(X[:, 0], X[:, 1], c=y_true, cmap='viridis', s=50, alpha=0.6)
axes[0].set_xlabel('Feature 1')
axes[0].set_ylabel('Feature 2')
axes[0].set_title('True Clusters (for reference only)')
axes[0].grid(True, alpha=0.3)
plt.colorbar(scatter1, ax=axes[0])

# Unlabeled data (what we actually have)
axes[1].scatter(X[:, 0], X[:, 1], c='gray', s=50, alpha=0.6)
axes[1].set_xlabel('Feature 1')
axes[1].set_ylabel('Feature 2')
axes[1].set_title('Unlabeled Data (What We Have)')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Key considerations
print("\n" + "="*50)
print("KEY CHALLENGES IN UNSUPERVISED LEARNING")
print("="*50)

print("""
1. NO GROUND TRUTH
- Can't directly evaluate correctness
- Need validation metrics (silhouette, inertia)
- Require domain knowledge

2. CHOOSING NUMBER OF CLUSTERS
- How many groups are there?
- Too few: Oversimplify, miss patterns
- Too many: Overcomplicate, find noise
- Use elbow method, silhouette score

3. INTERPRETING RESULTS
- What do clusters mean?
- Are they meaningful or just mathematical artifacts?
- Require domain expertise

4. SCALABILITY
- High dimensions, many samples
- Distance metrics become problematic
- Computational complexity important

5. PARAMETER SENSITIVITY
- Results depend on parameters
- Different parameters → different clusters
- Need to explore multiple settings
""")

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

print("="*50)
print("K-MEANS CLUSTERING")
print("="*50)

print("""
HOW K-MEANS WORKS:

1. INITIALIZE
- Randomly choose k points as initial centroids
- k = number of clusters (we choose this)

2. ASSIGN
- Assign each point to nearest centroid
- "Nearest" = smallest Euclidean distance

3. UPDATE
- Update centroids to mean of assigned points
- Centroid moves toward cluster center

4. REPEAT
- Repeat assign-update until convergence
- Stops when centroids don't move much

PROPERTIES:
- Local optimization (may not find global optimum)
- Sensitive to initialization
- Assumes spherical clusters
- Requires specifying k in advance
- Fast: O(n*k*i) where i = iterations
""")

# Generate data
np.random.seed(42)
X, y_true = make_blobs(n_samples=300, centers=3, n_features=2,
cluster_std=0.6, random_state=42)

# Standardize features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("\nEXAMPLE: CLUSTERING SYNTHETIC DATA")
print("-"*50)

print(f"Number of samples: {X_scaled.shape[0]}")
print(f"Number of features: {X_scaled.shape[1]}")

# Fit K-means
kmeans = KMeans(n_clusters=3, random_state=42, n_init=10)
clusters = kmeans.fit_predict(X_scaled)

print(f"\nK-means Results (k=3):")
print(f"Centroids shape: {kmeans.cluster_centers_.shape}")
print(f"Inertia (sum of squared distances): {kmeans.inertia_:.2f}")
print(f"\nCluster sizes:")
unique, counts = np.unique(clusters, return_counts=True)
for cluster_id, count in zip(unique, counts):
print(f" Cluster {cluster_id}: {count} samples")

# Predictions for new points
new_point = np.array([[0.5, 0.5]])
new_point_scaled = scaler.transform(new_point)
prediction = kmeans.predict(new_point_scaled)
distance_to_centroid = np.min(np.linalg.norm(new_point_scaled - kmeans.cluster_centers_, axis=1))

print(f"\nNew point prediction:")
print(f" Point (scaled): {new_point_scaled}")
print(f" Assigned to cluster: {prediction[0]}")
print(f" Distance to centroid: {distance_to_centroid:.3f}")

# Choosing k: Elbow Method
print("\n" + "="*50)
print("CHOOSING K: ELBOW METHOD")
print("="*50)

print("""
ELBOW METHOD:
1. Fit K-means for different k values
2. Plot inertia (within-cluster sum of squares) vs k
3. Look for "elbow" - point where decrease slows
4. Choose k at elbow
""")

inertias = []
silhouette_scores = []
K_range = range(1, 10)

from sklearn.metrics import silhouette_score

for k in K_range:
kmeans_k = KMeans(n_clusters=k, random_state=42, n_init=10)
clusters_k = kmeans_k.fit_predict(X_scaled)
inertias.append(kmeans_k.inertia_)

if k > 1: # Silhouette needs at least 2 clusters
silhouette_scores.append(silhouette_score(X_scaled, clusters_k))
else:
silhouette_scores.append(None)

print("\nInertia for different k values:")
for k, inertia in zip(K_range, inertias):
print(f" k={k}: {inertia:.2f}")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# Plot 1: Clustered data
scatter = axes[0, 0].scatter(X_scaled[:, 0], X_scaled[:, 1], c=clusters,
cmap='viridis', s=50, alpha=0.6)
axes[0, 0].scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],
c='red', marker='X', s=200, edgecolors='black', linewidths=2,
label='Centroids')
axes[0, 0].set_xlabel('Feature 1')
axes[0, 0].set_ylabel('Feature 2')
axes[0, 0].set_title('K-means Clustering (k=3)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
plt.colorbar(scatter, ax=axes[0, 0])

# Plot 2: Elbow curve
axes[0, 1].plot(K_range, inertias, 'b-o', linewidth=2, markersize=8)
axes[0, 1].axvline(x=3, color='red', linestyle='--', linewidth=2, label='Elbow at k=3')
axes[0, 1].set_xlabel('Number of Clusters (k)')
axes[0, 1].set_ylabel('Inertia')
axes[0, 1].set_title('Elbow Method')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: Silhouette scores
valid_k = [k for k, s in zip(K_range, silhouette_scores) if s is not None]
valid_scores = [s for s in silhouette_scores if s is not None]
axes[1, 0].plot(valid_k, valid_scores, 'g-o', linewidth=2, markersize=8)
axes[1, 0].axvline(x=3, color='red', linestyle='--', linewidth=2)
axes[1, 0].set_xlabel('Number of Clusters (k)')
axes[1, 0].set_ylabel('Silhouette Score')
axes[1, 0].set_title('Silhouette Score vs k')
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: Comparison of k values
for k in [2, 3, 4]:
kmeans_k = KMeans(n_clusters=k, random_state=42, n_init=10)
clusters_k = kmeans_k.fit_predict(X_scaled)

axes[1, 1].scatter(X_scaled[:, 0], X_scaled[:, 1], c=clusters_k,
alpha=0.3, s=30)

axes[1, 1].set_xlabel('Feature 1')
axes[1, 1].set_ylabel('Feature 2')
axes[1, 1].set_title('Comparison: Different k values')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Advantages and disadvantages
print("\n" + "="*50)
print("K-MEANS ADVANTAGES AND DISADVANTAGES")
print("="*50)

print("""
ADVANTAGES:
✓ Fast: O(nki) where n=samples, k=clusters, i=iterations
✓ Scalable: Works well with large datasets
✓ Simple to understand and implement
✓ Easy to interpret: cluster centers are meaningful
✓ Works well with roughly spherical clusters

DISADVANTAGES:
✗ Must specify k in advance
✗ Sensitive to initialization (may find local minima)
✗ Assumes spherical, similarly-sized clusters
✗ Affected by outliers
✗ Doesn't scale well with high dimensions
✗ All points must be assigned (no outlier detection)

WHEN TO USE:
- Large datasets
- Spherical cluster shapes
- Need speed
- Know or can estimate k

WHEN NOT TO USE:
- Complex cluster shapes
- Unknown number of clusters
- Many outliers
- Very high dimensions
""")

import numpy as np
from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

print("="*50)
print("HIERARCHICAL CLUSTERING")
print("="*50)

print("""
HIERARCHICAL CLUSTERING:
- Build tree of clusters (dendrogram)
- Two approaches: Agglomerative and Divisive

AGGLOMERATIVE (Bottom-Up):
1. Start with each point as own cluster
2. Merge closest clusters iteratively
3. Continue until one big cluster
4. Result: Tree showing merge sequence

DIVISIVE (Top-Down):
1. Start with all points in one cluster
2. Recursively split clusters
3. Continue until each point separate
4. Less common, more expensive

LINKAGE METHODS (how to measure cluster distance):
- Single: Distance between closest points
- Complete: Distance between farthest points
- Average: Average distance between all pairs
- Ward: Minimizes within-cluster variance
""")

# Generate data
np.random.seed(42)
X, y_true = make_blobs(n_samples=100, centers=4, n_features=2,
cluster_std=0.7, random_state=42)

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("\nEXAMPLE: HIERARCHICAL CLUSTERING")
print("-"*50)

print(f"Number of samples: {len(X_scaled)}")
print(f"Number of features: {X_scaled.shape[1]}")

# Hierarchical clustering with different linkages
linkages = ['ward', 'complete', 'average', 'single']

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

for idx, linkage_method in enumerate(linkages):
Z = linkage(X_scaled, method=linkage_method)

dendrogram(Z, ax=axes[idx], no_labels=True)
axes[idx].set_title(f'Dendrogram ({linkage_method.capitalize()} Linkage)')
axes[idx].set_xlabel('Sample Index')
axes[idx].set_ylabel('Distance')
axes[idx].axhline(y=5, color='red', linestyle='--', linewidth=2, label='Cut line')
axes[idx].legend()

plt.suptitle('Hierarchical Clustering with Different Linkages',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Use Ward linkage (most common)
print("\nUsing Ward linkage:")
Z_ward = linkage(X_scaled, method='ward')

# Agglomerative clustering
n_clusters = 4
hierarchical = AgglomerativeClustering(n_clusters=n_clusters, linkage='ward')
clusters = hierarchical.fit_predict(X_scaled)

print(f"\nResults (k={n_clusters}):")
unique, counts = np.unique(clusters, return_counts=True)
for cluster_id, count in zip(unique, counts):
print(f" Cluster {cluster_id}: {count} samples")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Dendrogram
dendrogram(Z_ward, ax=axes[0], no_labels=True)
axes[0].axhline(y=10, color='red', linestyle='--', linewidth=2, label='Cut for 4 clusters')
axes[0].set_title('Dendrogram (Ward Linkage)')
axes[0].set_xlabel('Sample Index')
axes[0].set_ylabel('Distance')
axes[0].legend()
axes[0].grid(True, alpha=0.3, axis='y')

# Clustered data
scatter = axes[1].scatter(X_scaled[:, 0], X_scaled[:, 1], c=clusters,
cmap='viridis', s=50, alpha=0.6)
axes[1].set_xlabel('Feature 1')
axes[1].set_ylabel('Feature 2')
axes[1].set_title(f'Hierarchical Clustering (k={n_clusters})')
axes[1].grid(True, alpha=0.3)
plt.colorbar(scatter, ax=axes[1])

plt.tight_layout()
plt.show()

# Advantages and disadvantages
print("\n" + "="*50)
print("HIERARCHICAL CLUSTERING ADVANTAGES AND DISADVANTAGES")
print("="*50)

print("""
ADVANTAGES:
✓ No need to specify k in advance (read dendrogram)
✓ Dendrogram provides interpretable tree structure
✓ Works with different cluster shapes
✓ Deterministic (no random initialization)
✓ Can use any distance metric
✓ Good for hierarchical data

DISADVANTAGES:
✗ Computationally expensive: O(n²) to O(n³)
✗ Can't undo merges (greedy algorithm)
✗ Sensitive to noise and outliers
✗ Different linkages give different results
✗ Not scalable to very large datasets
✗ All points must be assigned (no outlier handling)

WHEN TO USE:
- Small to medium datasets
- Hierarchical relationships matter
- Want to explore multiple k values
- Need interpretable structure

WHEN NOT TO USE:
- Very large datasets
- Need speed/scalability
- Have no hierarchical structure
- Working with high dimensions
""")

import numpy as np
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import NearestNeighbors
import matplotlib.pyplot as plt

print("="*50)
print("DBSCAN (Density-Based Spatial Clustering)")
print("="*50)

print("""
HOW DBSCAN WORKS:

1. CORE POINTS
- Points with >= min_samples neighbors within eps distance
- "Dense" regions

2. BORDER POINTS
- Non-core points within eps of core point
- On cluster boundary

3. NOISE POINTS
- Points not core or border
- Outliers/anomalies

4. CLUSTERING
- Group core points that are close
- Include border points
- Mark noise points as outliers (-1)

PARAMETERS:
- eps: Maximum distance between points (radius)
- min_samples: Minimum points in neighborhood for core point
Usually: 2*n_features or larger

PROPERTIES:
- Arbitrary cluster shapes
- Automatic outlier detection
- No need to specify number of clusters
- Good with varying density (if eps set right)
""")

# Generate data with non-convex shapes
print("\nEXAMPLE: CLUSTERING NON-CONVEX DATA")
print("-"*50)

# Create moons dataset (two crescents)
X, y_true = make_moons(n_samples=300, noise=0.05, random_state=42)
X_scaled = StandardScaler().fit_transform(X)

print(f"Generated moons dataset with {len(X)} samples")
print(f"True structure: Two crescents (K-means would fail here!)")

# Find optimal eps using k-distance graph
neighbors = NearestNeighbors(n_neighbors=5)
neighbors_fit = neighbors.fit(X_scaled)
distances, indices = neighbors_fit.kneighbors(X_scaled)
distances = np.sort(distances[:, -1], axis=0)

print("\nFinding optimal eps:")
# Optimal eps usually at "knee" of distance curve
print(f"Distance at 90th percentile: {np.percentile(distances, 90):.3f}")

# Fit DBSCAN
eps = 0.2
min_samples = 5

dbscan = DBSCAN(eps=eps, min_samples=min_samples)
clusters = dbscan.fit_predict(X_scaled)

n_clusters = len(set(clusters)) - (1 if -1 in clusters else 0)
n_outliers = list(clusters).count(-1)

print(f"\nDBSCAN Results (eps={eps}, min_samples={min_samples}):")
print(f"Number of clusters: {n_clusters}")
print(f"Number of outliers: {n_outliers}")
print(f"\nCluster sizes:")
for cluster_id in sorted(set(clusters)):
if cluster_id == -1:
print(f" Outliers: {list(clusters).count(cluster_id)}")
else:
print(f" Cluster {cluster_id}: {list(clusters).count(cluster_id)}")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# Plot 1: K-distance graph
axes[0, 0].plot(distances)
axes[0, 0].axhline(y=eps, color='red', linestyle='--', linewidth=2, label=f'eps={eps}')
axes[0, 0].set_xlabel('Points sorted by distance')
axes[0, 0].set_ylabel('5-nearest neighbor distance')
axes[0, 0].set_title('K-distance Graph (for choosing eps)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: DBSCAN clustering
scatter = axes[0, 1].scatter(X_scaled[:, 0], X_scaled[:, 1], c=clusters,
cmap='viridis', s=50, alpha=0.6)
# Highlight outliers
outlier_mask = clusters == -1
axes[0, 1].scatter(X_scaled[outlier_mask, 0], X_scaled[outlier_mask, 1],
c='red', marker='X', s=200, edgecolors='black',
linewidths=2, label='Outliers')
axes[0, 1].set_xlabel('Feature 1')
axes[0, 1].set_ylabel('Feature 2')
axes[0, 1].set_title(f'DBSCAN Clustering')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
plt.colorbar(scatter, ax=axes[0, 1])

# Plot 3: Different eps values
eps_values = [0.1, 0.2, 0.3]
for idx, eps_val in enumerate(eps_values):
ax = axes[1, idx // 2] if idx < 2 else axes[1, 1]

dbscan_eps = DBSCAN(eps=eps_val, min_samples=5)
clusters_eps = dbscan_eps.fit_predict(X_scaled)

ax.scatter(X_scaled[:, 0], X_scaled[:, 1], c=clusters_eps,
cmap='viridis', s=50, alpha=0.6)
n_clust = len(set(clusters_eps)) - (1 if -1 in clusters_eps else 0)
n_out = list(clusters_eps).count(-1)

ax.set_title(f'eps={eps_val} (k={n_clust}, outliers={n_out})')
ax.grid(True, alpha=0.3)

# Remove extra subplot
if len(eps_values) < 3:
fig.delaxes(axes[1, 1])

plt.suptitle('DBSCAN: Effect of eps Parameter', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Advantages and disadvantages
print("\n" + "="*50)
print("DBSCAN ADVANTAGES AND DISADVANTAGES")
print("="*50)

print("""
ADVANTAGES:
✓ Finds arbitrary cluster shapes
✓ Automatic outlier detection
✓ No need to specify number of clusters
✓ Works well with varying density (if eps chosen well)
✓ Deterministic (no random initialization)

DISADVANTAGES:
✗ Difficult to choose eps and min_samples
✗ Bad with clusters of varying density
✗ Sensitivity to parameters
✗ O(n²) to O(n log n) complexity
✗ High-dimensional data problematic

WHEN TO USE:
- Non-convex cluster shapes
- Need outlier detection
- Don't know number of clusters
- Varying cluster density acceptable

WHEN NOT TO USE:
- Uniform cluster density required
- Need to specify exact clusters
- Very high dimensions
- Large datasets (slow)
""")

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

print("="*50)
print("DIMENSIONALITY REDUCTION")
print("="*50)

print("""
WHY REDUCE DIMENSIONS?

1. VISUALIZATION
- Reduce to 2D or 3D for plotting
- Understand data structure

2. CURSE OF DIMENSIONALITY
- Too many features, not enough samples
- Model overfitting, poor generalization
- Distance metrics become meaningless
- Computational cost explodes

3. NOISE REDUCTION
- Remove irrelevant features
- Keep only important information
- Improve model performance

4. FEATURE EXTRACTION
- Create new features from old ones
- May be more interpretable

METHODS:

LINEAR:
- PCA: Unsupervised, finds variance directions
- Eigenvectors capture principal components
- Fast, interpretable

NONLINEAR:
- t-SNE: Excellent for visualization
- UMAP: Fast t-SNE alternative
- Autoencoders: Deep learning approach
""")

# Load data
iris = load_iris()
X = iris.data
y = iris.target
target_names = iris.target_names
feature_names = iris.feature_names

print(f"\nEXAMPLE: IRIS DATASET")
print("-"*50)

print(f"Original data shape: {X.shape}")
print(f"Features: {feature_names}")
print(f"Classes: {target_names}")

# High-dimensional visualization problem
print("\nChallenge:")
print("- 4 features (4D data)")
print("- Can't visualize 4D directly")
print("- Solution: Reduce to 2D for visualization")

# Standardize
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCA with different components
print("\n" + "="*50)
print("PCA: PRINCIPAL COMPONENT ANALYSIS")
print("="*50)

print("""
HOW PCA WORKS:

1. STANDARDIZE DATA
- Mean = 0, variance = 1

2. COMPUTE COVARIANCE MATRIX
- Show relationships between features

3. FIND EIGENVECTORS AND EIGENVALUES
- Eigenvectors: principal component directions
- Eigenvalues: variance explained per component

4. ORDER BY VARIANCE
- Largest eigenvalue = first principal component
- Explains most variance

5. SELECT COMPONENTS
- Keep top k components
- Usually retain 90-95% variance
""")

# Fit full PCA
pca_full = PCA()
X_pca_full = pca_full.fit_transform(X_scaled)

# Explained variance
explained_var = pca_full.explained_variance_ratio_
cumsum_var = np.cumsum(explained_var)

print(f"\nExplained Variance by Component:")
for i, (var, cumsum) in enumerate(zip(explained_var, cumsum_var)):
print(f" PC{i+1}: {var:.4f} ({var*100:.2f}%) [Cumulative: {cumsum*100:.2f}%]")

print(f"\nTo retain 90% variance: need {np.argmax(cumsum_var >= 0.9) + 1} components")
print(f"To retain 95% variance: need {np.argmax(cumsum_var >= 0.95) + 1} components")

# Fit 2-component PCA
pca_2 = PCA(n_components=2)
X_pca_2 = pca_2.fit_transform(X_scaled)

print(f"\n2-component PCA:")
print(f" Explained variance: {pca_2.explained_variance_ratio_.sum()*100:.2f}%")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# Plot 1: Scree plot (variance explained)
axes[0, 0].plot(range(1, len(explained_var)+1), explained_var, 'b-o', linewidth=2)
axes[0, 0].set_xlabel('Principal Component')
axes[0, 0].set_ylabel('Explained Variance Ratio')
axes[0, 0].set_title('Scree Plot')
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Cumulative explained variance
axes[0, 1].plot(range(1, len(cumsum_var)+1), cumsum_var, 'g-o', linewidth=2)
axes[0, 1].axhline(y=0.9, color='red', linestyle='--', label='90% variance')
axes[0, 1].axhline(y=0.95, color='orange', linestyle='--', label='95% variance')
axes[0, 1].set_xlabel('Number of Components')
axes[0, 1].set_ylabel('Cumulative Explained Variance')
axes[0, 1].set_title('Cumulative Explained Variance')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: 2D PCA projection
for i, target in enumerate(np.unique(y)):
axes[1, 0].scatter(X_pca_2[y == target, 0], X_pca_2[y == target, 1],
label=target_names[target], s=50, alpha=0.6)
axes[1, 0].set_xlabel(f'PC1 ({pca_2.explained_variance_ratio_[0]*100:.1f}%)')
axes[1, 0].set_ylabel(f'PC2 ({pca_2.explained_variance_ratio_[1]*100:.1f}%)')
axes[1, 0].set_title('Iris Data Projected to 2 Principal Components')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: Component loadings (which original features matter)
loadings = pca_2.components_.T * np.sqrt(pca_2.explained_variance_)
loading_df = pd.DataFrame(loadings, columns=['PC1', 'PC2'], index=feature_names)

loading_df.plot(kind='barh', ax=axes[1, 1])
axes[1, 1].set_xlabel('Loading')
axes[1, 1].set_title('Feature Contributions to Principal Components')
axes[1, 1].grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.show()

# Interpretation
print("\n" + "="*50)
print("INTERPRETING PCA")
print("="*50)

print("\nComponent Loadings (which features matter):")
print(loading_df)

print("\nInterpretation:")
print("- High loading: Feature important for that component")
print("- Positive loading: Feature increases with component")
print("- Negative loading: Feature decreases with component")
print("- Loadings help understand what components represent")

import numpy as np
from sklearn.datasets import load_iris, load_digits
from sklearn.preprocessing import StandardScaler
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt

print("="*50)
print("ADVANCED DIMENSIONALITY REDUCTION")
print("="*50)

print("""
PCA LIMITATIONS:
- Linear transformation only
- Can miss nonlinear structures
- Not ideal for visualization
- Example: Swiss roll data (PCA fails)

NONLINEAR ALTERNATIVES:

1. t-SNE (t-Distributed Stochastic Neighbor Embedding)
PROS:
- Excellent for visualization
- Preserves local structure well
- Great for exploratory analysis

CONS:
- Slow: O(n²) or O(n log n)
- Non-deterministic (use random_state)
- Hyperparameters matter: perplexity, learning_rate
- Not good for prediction (can't transform new data easily)
- Can create misleading clusters

2. UMAP (Uniform Manifold Approximation and Projection)
PROS:
- Faster than t-SNE
- Better preserves global structure
- Can transform new data
- Works well for high dimensions

CONS:
- Newer, less established
- Still hyperparameter sensitive
- Interpretation can be tricky

3. AUTOENCODERS
PROS:
- Deep learning approach
- Can learn complex transformations
- Good for feature extraction

CONS:
- Requires neural network knowledge
- More parameters to tune
- Can overfit
""")

# Load a more interesting dataset
digits = load_digits()
X = digits.data
y = digits.target

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print(f"\nEXAMPLE: HANDWRITTEN DIGITS")
print("-"*50)

print(f"Original shape: {X.shape} (64 features = 8x8 pixels)")
print(f"Classes: 10 (digits 0-9)")

# t-SNE
print("\nFitting t-SNE (this may take a moment)...")

tsne = TSNE(n_components=2, random_state=42, perplexity=30, n_iter=1000)
X_tsne = tsne.fit_transform(X_scaled[:500]) # Use subset for speed
y_subset = y[:500]

print("t-SNE fitting complete!")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# PCA
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled[:500])

scatter1 = axes[0].scatter(X_pca[:, 0], X_pca[:, 1], c=y_subset,
cmap='tab10', s=50, alpha=0.6)
axes[0].set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)')
axes[0].set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)')
axes[0].set_title('PCA Projection')
axes[0].grid(True, alpha=0.3)
plt.colorbar(scatter1, ax=axes[0])

# t-SNE
scatter2 = axes[1].scatter(X_tsne[:, 0], X_tsne[:, 1], c=y_subset,
cmap='tab10', s=50, alpha=0.6)
axes[1].set_xlabel('t-SNE 1')
axes[1].set_ylabel('t-SNE 2')
axes[1].set_title('t-SNE Projection')
axes[1].grid(True, alpha=0.3)
plt.colorbar(scatter2, ax=axes[1])

plt.suptitle('PCA vs t-SNE: Digits Dataset', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("\nComparison:")
print("- PCA: Fast, global structure, linear")
print("- t-SNE: Slow, local structure, nonlinear")
print("- t-SNE shows clusters better but distorts global distances")

print("\n" + "="*50)
print("CHOOSING DIMENSIONALITY REDUCTION METHOD")
print("="*50)

print("""
USE PCA IF:
✓ Need speed and interpretability
✓ Want explained variance metric
✓ Need to transform new data
✓ Linear relationships sufficient

USE t-SNE IF:
✓ Need visualization only
✓ Have time for computation
✓ Want to see local clusters
✓ Non-linear structure important

USE UMAP IF:
✓ Need balance of speed and quality
✓ Want global + local structure
✓ Need to transform new data
✓ Have large dataset

USE AUTOENCODERS IF:
✓ Expert in deep learning
✓ Complex nonlinear relationships
✓ Feature learning important
✓ Have lots of data and compute
""")

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import (silhouette_score, davies_bouldin_score,
calinski_harabasz_score)
import pandas as pd

print("="*50)
print("EVALUATING CLUSTERING QUALITY")
print("="*50)

print("""
CHALLENGE:
- No ground truth labels in unsupervised learning
- Can't directly measure correctness
- Need internal validation metrics

CLUSTERING EVALUATION METRICS:

1. SILHOUETTE SCORE (most common)
- Range: -1 to 1
- 1: Perfect clustering
- 0: Overlapping clusters
- -1: Wrong clustering
- Measures: How similar point is to own cluster vs others

2. DAVIES-BOULDIN INDEX
- Range: 0 to ∞
- Lower is better
- Measures: Ratio of distances

3. CALINSKI-HARABASZ INDEX
- Range: 0 to ∞
- Higher is better
- Measures: Ratio of between to within cluster variance

4. INERTIA (for K-means)
- Sum of squared distances to nearest centroid
- Lower is better
- Can't compare across different k
""")

# Generate data
np.random.seed(42)
X, y_true = make_blobs(n_samples=300, centers=3, cluster_std=0.6, random_state=42)
X_scaled = StandardScaler().fit_transform(X)

print("\nEXAMPLE: EVALUATE K-MEANS WITH DIFFERENT K")
print("-"*50)

results = []

for k in range(2, 10):
kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
clusters = kmeans.fit_predict(X_scaled)

silhouette = silhouette_score(X_scaled, clusters)
davies_bouldin = davies_bouldin_score(X_scaled, clusters)
calinski = calinski_harabasz_score(X_scaled, clusters)
inertia = kmeans.inertia_

results.append({
'k': k,
'Silhouette': silhouette,
'Davies-Bouldin': davies_bouldin,
'Calinski-Harabasz': calinski,
'Inertia': inertia
})

results_df = pd.DataFrame(results)
print(results_df.to_string(index=False))

# Find optimal k
optimal_k_silhouette = results_df.loc[results_df['Silhouette'].idxmax(), 'k']
optimal_k_davies = results_df.loc[results_df['Davies-Bouldin'].idxmin(), 'k']
optimal_k_calinski = results_df.loc[results_df['Calinski-Harabasz'].idxmax(), 'k']

print(f"\nOptimal k by:")
print(f" Silhouette Score: k={int(optimal_k_silhouette)}")
print(f" Davies-Bouldin Index: k={int(optimal_k_davies)}")
print(f" Calinski-Harabasz: k={int(optimal_k_calinski)}")

# Visualization
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Silhouette (higher is better)
axes[0, 0].plot(results_df['k'], results_df['Silhouette'], 'b-o', linewidth=2)
axes[0, 0].axvline(x=optimal_k_silhouette, color='red', linestyle='--', alpha=0.7)
axes[0, 0].set_xlabel('Number of Clusters (k)')
axes[0, 0].set_ylabel('Silhouette Score')
axes[0, 0].set_title('Silhouette Score (Higher is Better)')
axes[0, 0].grid(True, alpha=0.3)

# Davies-Bouldin (lower is better)
axes[0, 1].plot(results_df['k'], results_df['Davies-Bouldin'], 'g-o', linewidth=2)
axes[0, 1].axvline(x=optimal_k_davies, color='red', linestyle='--', alpha=0.7)
axes[0, 1].set_xlabel('Number of Clusters (k)')
axes[0, 1].set_ylabel('Davies-Bouldin Index')
axes[0, 1].set_title('Davies-Bouldin Index (Lower is Better)')
axes[0, 1].grid(True, alpha=0.3)

# Calinski-Harabasz (higher is better)
axes[1, 0].plot(results_df['k'], results_df['Calinski-Harabasz'], 'm-o', linewidth=2)
axes[1, 0].axvline(x=optimal_k_calinski, color='red', linestyle='--', alpha=0.7)
axes[1, 0].set_xlabel('Number of Clusters (k)')
axes[1, 0].set_ylabel('Calinski-Harabasz Score')
axes[1, 0].set_title('Calinski-Harabasz Index (Higher is Better)')
axes[1, 0].grid(True, alpha=0.3)

# Inertia (for reference)
axes[1, 1].plot(results_df['k'], results_df['Inertia'], 'c-o', linewidth=2)
axes[1, 1].set_xlabel('Number of Clusters (k)')
axes[1, 1].set_ylabel('Inertia')
axes[1, 1].set_title('Inertia (Lower is Better, but use Elbow)')
axes[1, 1].grid(True, alpha=0.3)

plt.suptitle('Clustering Evaluation Metrics', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("\n" + "="*50)
print("BEST PRACTICES FOR UNSUPERVISED LEARNING")
print("="*50)

print("""
1. DATA PREPARATION
✓ Standardize/normalize features
✓ Remove highly correlated features
✓ Handle missing values
✓ Outlier detection/removal appropriate
✗ Don't use target variable

2. ALGORITHM SELECTION
✓ Start with simple methods (K-means)
✓ Try multiple algorithms
✓ Consider interpretability
✗ Don't assume one algorithm works for all

3. PARAMETER TUNING
✓ Use multiple evaluation metrics
✓ Cross-validate results
✓ Try range of parameters
✓ Visualization helps
✗ Don't rely on single metric

4. VALIDATION
✓ Use internal validation metrics
✓ Domain knowledge assessment
✓ Stability across runs
✓ Interpretability check
✗ Don't over-interpret results

5. COMMON PITFALLS
✗ Expecting ground truth accuracy
✗ Using only one evaluation metric
✗ Not scaling features
✗ Over-interpreting clusters
✗ Ignoring outliers
✗ Not validating with domain experts

6. DIMENSIONALITY REDUCTION
✓ Choose method based on goal
✓ Retain sufficient variance (90-95%)
✓ Interpret components
✓ Validate on downstream task
✗ Don't remove features blindly
""")