Can I compare the output probabilities of two machine learning models? I'm sorry if this is a silly question.
Suppose there are two logistic regression models $M_1$ and $M_2$ trained on the same (or similar) dataset, and their outputs of given input $x$ are $P_{M_1}(y \mid x)$ and $P_{M_2}(y \mid x)$, respectively.
What I feel confused about is, for any given input $x$, could I determine which model's prediction is more reliable by comparing this output probability? For example, if $P_{M_1}(y=k \mid x) > P_{M_2}(y=k \mid x)$, then the output of $M_1$ is more reliable? If cannot, how to compare these two models and make a decision that which one should be chosen?

This sounds appealing at first. After all, if one model predicts $0.6$ and another $0.9$, you’d trust the confident mode making the confident $0.9$ prediction and not the wishy-washy $0.6$ prediction, right? The trouble is that the right answer might be that your case has ambiguity, and about $60\%$ of such cases will go for one class and $40\%$ for the other. You might be interested in evaluating model calibration and performance on strictly proper scoring-rules, such as log loss and Brier score, that seek out the correct probability predictions.
Following an answer by Stéphane Laurent, let's do a simulation in R to see what it means to predict the correct (unobservable) probability.
set.seed(2021)
N <- 1000
x <- runif(N, -2, 2)
z <- x
pr <- 1/(1 + exp(-z)) # This is the true probability!
y <- rbinom(N, 1, pr)
L <- glm(y ~ x, family = binomial)

In this simulation, the y variable mimics the discrete categories that we would observe. However, there are y values of $0$ that correspond to $P(0)<0.5$, just by the luck of the draw. That is, they turned out to be $0$, even though, given x, the result is more likely to be $1$ than $0$. Therefore, we want our model to predict probability values that come close to pr. In expected value, strictly proper scoring rules like log loss and Brier score are minimized by predicting pr.
