How can I validate that aleatoric uncertainty is well-trained?

I'm trying to implement a neural network that learns aleatoric uncertainty for regression task using pytorch according to Kendall et al.

It seems that the training have done well, which means learning curve converged and predicted mean values fitted to desired ground truth values. In this situation, while predicted mean values can be validated by comparing them to ground truth values using error metric such as MSE or something, I have no idea how to validate uncertainties.

The model outputs some variances that represents uncertainties, but how can I know this values are well-trained?

Here is what I've tried:

Aleatoric uncertainty represents the variance of y given the same Xs. But my dataset rarely have the same Xs. So, I did k-means clustering to collect similar Xs as much as possible, and then calculated standard deviations for each cluster to compare with aleatoric uncertainties. As a result, aleatoric uncertainties are much larger than the standard deviations for every case.

Judging from this result, should it be considered that aleatoric uncertainty was learned incorrectly?

Or is there any better way to evaluate aleatoric uncertainty?

Appendix:

I have 3 features, and I'm using MinMaxScaler because they have different scales.

sc_X = MinMaxScaler()
sc_Y = MinMaxScaler()
dataset_X = op_7GT_max[gt_input_cols]
dataset_Y = op_7GT_max[gt_output_cols]
X_train, X_test, y_train, y_test = train_test_split(dataset_X, dataset_Y, test_size=0.2)
X_train_sc = sc_X.fit_transform(X_train.values)
y_train_sc = sc_Y.fit_transform(y_train.values)
X_test_sc = sc_X.transform(X_test.values)
y_test_sc = sc_Y.transform(y_test.values)

Here is my model:

class Model(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(nn.Linear(3, 16), nn.ReLU(), nn.Dropout(0.2),
                                 nn.Linear(16, 32), nn.ReLU(), nn.Dropout(0.2),
                                 nn.Linear(32, 16), nn.ReLU())
    
        self.out_mean = nn.Sequential(nn.Linear(16, 1))
        self.out_logstd = nn.Sequential(nn.Linear(16, 1))
    
    def get_loss(self, y, mean, logstd):

        negative_log_likelihood = logstd + 0.5 * np.log(2 * np.pi) + 0.5 * ((y - mean) ** 2) / 
                                                                  (torch.exp(logstd) ** 2)
        return torch.mean(negative_log_likelihood)
    
    def forward(self, x):
        features = self.net(x)
        mean = self.out_mean(features)
        logstd = self.out_logstd(features)
    
        return mean, logstd

    def inference(self, x, test_num=1000):
    
        mus = []
        sigmas = []
        for _ in range(test_num):
            mu, sigma = self.forward(x)
            mus.append(mu)
            sigmas.append(sigma)
    
        mus = torch.cat(mus, dim=1)
        sigmas = torch.cat(sigmas, dim=1)
    
        mean_mu = mus.mean(dim=1)
        std_mu = mus.std(dim=1)
    
        mean_sigma = sigmas.mean(dim=1)
        std_sigma = sigmas.std(dim=1)
    
        return mean_mu, std_mu, mean_sigma, std_sigma

Here is how I get predicted value and uncertainties, and how they are inverse transformed. _mean_mu_inv is predicted value, _std_mu_inv is epistemic uncertainty, and _mean_sigma_inv is aleatoric uneratainty which I interested in.

_X_test = torch.from_numpy(X_test_sc).float().to(device)

# Input X_test to model and get mean and logstd of Y
mean_mu, std_mu, mean_sigma, std_sigma = model.inference(_X_test)
_mean_mu = mean_mu.data.cpu().numpy()
_std_mu = std_mu.data.cpu().numpy()
_mean_sigma = np.exp(mean_sigma.data.cpu().numpy())
print(_mean_mu.shape, _std_mu.shape, _mean_sigma.shape)

# Inverse scaling
_mean_mu_inv = sc_Y.inverse_transform(_mean_mu.reshape(-1, 1)).flatten()
_std_mu_inv = abs(sc_Y.data_max_ - sc_Y.data_min_) * _std_mu
_mean_sigma_inv = abs(sc_Y.data_max_ - sc_Y.data_min_) * _mean_sigma