sta5635_hw12.py

# -*- coding: utf-8 -*-
"""STA5635_HW12.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1oo9gw5Naw0njiiua-BjwcK6R37UjMELR
"""

#Imported Packages
import numpy as np
import scipy.io
import os
from sklearn.preprocessing import StandardScaler
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset

!unzip features_640.zip
!unzip features_val_640.zip



# Define the mapping of file names to class labels
class_mapping = {
    "n01614925": 0,
    "n01689811": 1,
    "n01753488": 2,
    "n02017213": 3,
    "n02097130": 4,
    "n02165456": 5,
    "n02804414": 6,
    "n02909870": 7,
    "n03706229": 8,
    "n04258138": 9,
}

train_path = 'features_640'
test_path = 'features_val_640'

# Extract unique class names from training and test folders
all_files = os.listdir(train_path) + os.listdir(test_path)
class_names = sorted(set([filename.split(".")[0] for filename in all_files]))

# Create a mapping from class names to numeric labels
class_mapping = {name: idx for idx, name in enumerate(class_names)}

# Initialize empty lists for training data and labels
X_train_list = []
y_train_list = []

# Load training data
for filename in os.listdir(train_path):
    if filename.endswith(".mat"):
        class_name = filename.split(".")[0]
        class_label = class_mapping[class_name]
        file_path = os.path.join(train_path, filename)
        data = scipy.io.loadmat(file_path)
        X_train_list.append(data['feature'])
        y_train_list.append(np.full(data['feature'].shape[0], class_label))

# Concatenate training data
X_train = np.vstack(X_train_list)
y_train = np.concatenate(y_train_list)

# Load test data
X_test_list = []
y_test_list = []

for filename in os.listdir(test_path):
    if filename.endswith(".mat"):
        class_name = filename.split(".")[0]
        class_label = class_mapping[class_name]
        file_path = os.path.join(test_path, filename)
        data = scipy.io.loadmat(file_path)
        X_test_list.append(data['feature'])
        y_test_list.append(np.full(data['feature'].shape[0], class_label))

# Concatenate test data
X_test = np.vstack(X_test_list)
y_test = np.concatenate(y_test_list)

# Standardize training data
scaler = StandardScaler()
X_train_standardized = scaler.fit_transform(X_train)

# Standardize test data using training mean and std
X_test_standardized = scaler.transform(X_test)

# Outputs: X_train_standardized, y_train, X_test_standardized, y_test
print(f"Training set shape: {X_train_standardized.shape}")
print(f"Training labels shape: {y_train.shape}")
print(f"Test set shape: {X_test_standardized.shape}")
print(f"Test labels shape: {y_test.shape}")

#Part A

# Convert data to PyTorch tensors
X_train_tensor = torch.tensor(X_train_standardized, dtype=torch.float32)
y_train_tensor = torch.tensor(y_train, dtype=torch.long)
X_test_tensor = torch.tensor(X_test_standardized, dtype=torch.float32)
y_test_tensor = torch.tensor(y_test, dtype=torch.long)

# Create datasets and dataloaders
train_dataset = TensorDataset(X_train_tensor, y_train_tensor)
train_loader = DataLoader(train_dataset, batch_size=64, shuffle=True)

# Define the linear classifier
class LinearClassifier(nn.Module):
    def __init__(self, input_dim, output_dim):
        super(LinearClassifier, self).__init__()
        self.linear = nn.Linear(input_dim, output_dim)

    def forward(self, x):
        return self.linear(x)

# Initialize the model, loss function, optimizer, and scheduler
input_dim = X_train_tensor.shape[1]
output_dim = len(np.unique(y_train))
model = LinearClassifier(input_dim, output_dim)
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(model.parameters(), lr=0.001)
scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=10, gamma=0.1)

# Training loop
num_epochs = 20
for epoch in range(num_epochs):
    for i, (inputs, labels) in enumerate(train_loader):
        optimizer.zero_grad()
        outputs = model(inputs)
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()
    scheduler.step()
    print(f'Epoch [{epoch+1}/{num_epochs}], Loss: {loss.item():.4f}')

# Evaluate the model
with torch.no_grad():
    train_outputs = model(X_train_tensor)
    _, train_predicted = torch.max(train_outputs, 1)
    train_misclass_error = 1 - (train_predicted == y_train_tensor).sum().item() / len(y_train_tensor)

    test_outputs = model(X_test_tensor)
    _, test_predicted = torch.max(test_outputs, 1)
    test_misclass_error = 1 - (test_predicted == y_test_tensor).sum().item() / len(y_test_tensor)

print(f"Training Misclassification Error: {train_misclass_error:.4f}")
print(f"Test Misclassification Error: {test_misclass_error:.4f}")

#Part B
def probabilistic_pca(X, q):
    n_samples = X.shape[0]
    n_features = X.shape[1]

    # Calculate the mean of the data
    mean = np.mean(X, axis=0)

    # Center the data
    X_centered = X - mean

    # Calculate the covariance matrix
    covariance_matrix = np.cov(X_centered, rowvar=False)

    # Perform eigenvalue decomposition
    eigenvalues, eigenvectors = np.linalg.eigh(covariance_matrix)

    # Sort eigenvalues and eigenvectors in descending order
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # Select the top q principal components
    W = eigenvectors[:, :q]

    # Calculate the variance
    variance = np.mean(eigenvalues[q:])

    return mean, W, variance


# Assuming you have class labels in y_train
unique_classes = np.unique(y_train)
sigma_squared_values = {}
means = {}
for k in unique_classes:
    # Get data for class k
    class_data = X_train_standardized[y_train == k]

    # Perform PPCA
    mean_k, W_k, sigma2_k = probabilistic_pca(class_data, q=20)

    # Store the variance for class k
    sigma_squared_values[k] = sigma2_k
    print(f"Class {k}: sigma2 = {sigma2_k}")

#Part B

def probabilistic_pca(X, q):
    n_samples = X.shape[0]
    n_features = X.shape[1]

    # Calculate the mean of the data
    mean = np.mean(X, axis=0)

    # Center the data
    X_centered = X - mean

    # Calculate the covariance matrix
    covariance_matrix = np.cov(X_centered, rowvar=False)

    # Perform eigenvalue decomposition
    eigenvalues, eigenvectors = np.linalg.eigh(covariance_matrix)

    # Sort eigenvalues and eigenvectors in descending order
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # Select the top q principal components
    W = eigenvectors[:, :q]

    # Calculate the variance
    variance = np.mean(eigenvalues[q:])

    return mean, W, variance


# Assuming you have class labels in y_train
unique_classes = np.unique(y_train)
sigma_squared_values = {}
means = {}  # Initialize means dictionary
Ws = {} # Intialize Ws dictionary to store W_k for each class
for k in unique_classes:
    # Get data for class k
    class_data = X_train_standardized[y_train == k]

    # Perform PPCA
    mean_k, W_k, sigma2_k = probabilistic_pca(class_data, q=20)

    # Store the variance, mean and W for class k
    sigma_squared_values[k] = sigma2_k
    means[k] = mean_k  # Store the mean
    Ws[k] = W_k # Store W
    print(f"Class {k}: sigma2 = {sigma2_k}")

def mahalanobis_distance(x, mean, W, variance):
    x_minus_mu = x - mean
    sigma = np.dot(W, W.T) + variance * np.identity(W.shape[0])
    inv_sigma = np.linalg.inv(sigma)
    distance = np.dot(x_minus_mu.T, np.dot(inv_sigma, x_minus_mu))
    return distance

# Classify test data using Mahalanobis distances
predicted_labels = []
for x in X_test_standardized:
    distances = []
    for k in unique_classes:
        distance = mahalanobis_distance(x, means[k], Ws[k], sigma_squared_values[k])
        distances.append(distance)
    predicted_labels.append(unique_classes[np.argmin(distances)])

# Calculate the misclassification error
ppca_test_misclass_error = 1 - np.sum(np.array(predicted_labels) == y_test) / len(y_test)
print(f"PPCA Test Misclassification Error: {ppca_test_misclass_error:.4f}")

# Classify training data using Mahalanobis distances
predicted_labels_train = []
for x in X_train_standardized:
    distances = []
    for k in unique_classes:
        distance = mahalanobis_distance(x, means[k], Ws[k], sigma_squared_values[k])
        distances.append(distance)
    predicted_labels_train.append(unique_classes[np.argmin(distances)])

# Calculate the misclassification error for training data
ppca_train_misclass_error = 1 - np.sum(np.array(predicted_labels_train) == y_train) / len(y_train)
print(f"PPCA Train Misclassification Error: {ppca_train_misclass_error:.4f}")

def mahalanobis_distance(x, mean, W, variance):
    """Calculates the Mahalanobis distance."""
    x_minus_mu = x - mean
    sigma = np.dot(W, W.T) + variance * np.identity(W.shape[0])
    inv_sigma = np.linalg.inv(sigma)
    distance = np.dot(x_minus_mu.T, np.dot(inv_sigma, x_minus_mu))
    return distance

# Classify training data using Mahalanobis distances
predicted_labels_train = []
for x in X_train_standardized:
    distances = []
    for k in unique_classes:
        distance = mahalanobis_distance(x, means[k], Ws[k], sigma_squared_values[k])
        distances.append(distance)
    predicted_labels_train.append(unique_classes[np.argmin(distances)])

# Calculate the misclassification error for training data
ppca_train_misclass_error = 1 - np.sum(np.array(predicted_labels_train) == y_train) / len(y_train)
print(f"PPCA Train Misclassification Error: {ppca_train_misclass_error:.4f}")