import numpy as np
from sklearn.linear_model import BayesianRidge
from sklearn import datasets

# on main
X, y = datasets.make_regression(n_samples=10, n_features=20)
n_features = X.shape[1]
reg = BayesianRidge(fit_intercept=False).fit(X, y)
covariance_matrix = np.linalg.inv(
    reg.lambda_ * np.identity(n_features) + reg.alpha_ * np.dot(X.T, X)
)
np.allclose(reg.sigma_, covariance_matrix)