Para problemas categóricos, você usaria a classe CategoricalNB (documentação). Depois de treinar o modelo, as probabilidades estão no atributo feature_log_prob_, só que aqui você vai encontrar o log das probabilidades. Para converter para escala decimal, use a função np.exp.
Exemplo com dados inventados:
rng = np.random.RandomState(1)
X = rng.randint(3, size=(10, 4))
y = np.array([1, 2, 2, 2, 1, 2, 1, 1, 1, 2])
from sklearn.naive_bayes import CategoricalNB
clf = CategoricalNB(force_alpha=True)
clf.fit(X, y)
print(clf.feature_log_prob_)
# [array([[-2.07944154, -0.69314718, -0.98082925],
[-1.38629436, -0.69314718, -1.38629436]]),
array([[-0.69314718, -1.38629436, -1.38629436],
[-0.98082925, -0.98082925, -1.38629436]]),
array([[-0.98082925, -0.98082925, -1.38629436],
[-0.47000363, -1.38629436, -2.07944154]]),
array([[-1.38629436, -0.69314718, -1.38629436],
[-1.38629436, -1.38629436, -0.69314718]])]
print(np.exp(np.array(clf.feature_log_prob_))
# array([[[0.125, 0.5 , 0.375],
[0.25 , 0.5 , 0.25 ]],
[[0.5 , 0.25 , 0.25 ],
[0.375, 0.375, 0.25 ]],
[[0.375, 0.375, 0.25 ],
[0.625, 0.25 , 0.125]],
[[0.25 , 0.5 , 0.25 ],
[0.25 , 0.25 , 0.5 ]]])
O array resultante tem shape (4, 2, 3) porque o exemplo tem 4 variáveis, 2 classes target e 3 categorias por variável.