GPU 上具有 BIC 或 AIC 的高斯混合模型

Question

scikit-learn 中有 GMM 的 bic/aic 标准，但我想将我的数据拟合到 GPU 上。 我发现 GMM 在 CuPy（cuda numpy 包装器）中实现，但它没有 bic/aic 标准。（https://github.com/cupy/cupy/blob/master/examples/gmm/gmm.py）

如何对该代码实施 bic/aic 标准？或者有好的图书馆吗？请帮忙。

Answer 1

你在 pytorch 中的代码的想法，也许并不完美，但一个起点：

import torch
import math

def bic_aic(log_likelihood, n_params, n_samples):
    bic = -2 * log_likelihood + n_params * math.log(n_samples)
    aic = 2 * n_params - 2 * log_likelihood
    return bic, aic

def estimate_log_prob(X, inv_cov, means, cov_type='full'):
    n_features = X.shape[1]
    if cov_type == 'full':
    log_det = torch.sum(torch.log(inv_cov))
    precisions = inv_cov ** 2
    elif cov_type == 'diag':
    log_det = torch.sum(torch.log(inv_cov))
    precisions = inv_cov
    elif cov_type == 'spherical':
    log_det = n_features * torch.log(inv_cov)
    precisions = inv_cov

    log_prob = torch.sum((means ** 2 * precisions), 1) - 2 * torch.matmul(X, (means * precisions).T) + torch.matmul(X ** 2, precisions.T)
    return -0.5 * (n_features * math.log(2 * math.pi) + log_prob) + log_det

def e_step(X, inv_cov, means, weights, cov_type='full'):
    weighted_log_prob = estimate_log_prob(X, inv_cov, means, cov_type=cov_type) + torch.log(weights)
    log_prob_norm = torch.logsumexp(weighted_log_prob, dim=1)
    log_resp = weighted_log_prob - log_prob_norm[:, None]
    return torch.mean(log_prob_norm), log_resp

def m_step(X, resp, cov_type='full'):
    nk = torch.sum(resp, axis=0)
    means = torch.matmul(resp.T, X) / nk[:, None]
    
    if cov_type == 'full':
    covariances = torch.matmul((resp * (X - means[0]).T), (X - means[0])) / nk[0]
    elif cov_type == 'diag':
    covariances = torch.sum(resp * (X - means) ** 2, axis=0) / nk
    elif cov_type == 'spherical':
    covariances = torch.sum(resp * (X - means) ** 2) / (nk * X.shape[1])
    
    return nk / len(X), means, covariances

def train_gmm(X, max_iter, tol, means, covariances, cov_type='full'):
    lower_bound = -float('inf')
    converged = False

    n_samples, n_features = X.shape
    n_params = 3 * n_features  # Modify based on your cov_type
    
    weights = torch.tensor([0.5, 0.5], dtype=torch.float32)
    inv_cov = 1 / torch.sqrt(covariances)
    
    for n_iter in range(max_iter):
    prev_lower_bound = lower_bound
    log_prob_norm, log_resp = e_step(X, inv_cov, means, weights, cov_type=cov_type)
    weights, means, covariances = m_step(X, torch.exp(log_resp), cov_type=cov_type)
    
    inv_cov = 1 / torch.sqrt(covariances)
    lower_bound = log_prob_norm

    bic_score, aic_score = bic_aic(lower_bound.item(), n_params, n_samples)
    print(f"Iteration {n_iter}, BIC: {bic_score}, AIC: {aic_score}")

    change = lower_bound - prev_lower_bound
    if abs(change) < tol:
        converged = True
        break

    if not converged:
    print('Failed to converge. Increase max-iter or tol.')
    
    return inv_cov, means, weights, covariances

# Synthetic data and initial parameters (Replace with your own)
X = torch.rand(500, 2)
initial_means = torch.rand(2, 2)
initial_covariances = torch.rand(2, 2)

# Train GMM
max_iter = 100
tol = 1e-3
inv_cov, means, weights, covariances = train_gmm(X, max_iter, tol, initial_means, initial_covariances, cov_type='full')