Why do we model noise in linear regression but not logistic regression? The canonical probabilistic interpretation of linear regression is that $y$ is equal to $\theta^Tx$, plus a Gaussian noise random variable $\epsilon$.
However, in standard logistic regression, we don't consider noise (e.g. random bit flips with probability $p$) of the label $y$. Why is that?
 A: Short answer: we do, just implicitly.

A possibly more enlightening way of looking at things is the following.
In Ordinary Least Squares, we can consider that we do not model the errors or noise as $N(0,\sigma^2)$ distributed, but we model the observations as $N(x\beta,\sigma^2)$ distributed.
(Of course, this is precisely the same thing, just looking at it in two different ways.)
Now the analogous statement for logistic regression becomes clear: here, we model the observations as Bernoulli distributed with parameter $p(x)=\frac{1}{1+e^{-x\beta}}$.
We can flip this last way of thinking around if we want: we can indeed say that we are modeling the errors in logistic regression. Namely, we are modeling them as "the difference between a Bernoulli distributed variable with parameter $p(x)$ and $p(x)$ itself".
This is just very unwieldy, and this distribution does not have a name, plus the error here depends on our independent variables $x$ (in contrast to the homoskedasticity assumption in OLS, where the error is independent of $x$), so this way of looking at things is just not used as often.
A: To supplement Stephan's answer, similar to how in linear regression the target $y$ is generated by a ''systematic'' component involving $x$ and an independent ''noise'' component,  in logistic regression (and softmax regression more generally) you can actually also think of the target $y$ as computed by the following operation involving $x$ and some noise $\epsilon$:
$$ y = \arg \max_{i \in \{0, 1\}} (\alpha_i  + \epsilon_i)$$
where $\alpha_0 = 0, \alpha_1 = \theta^T x$, and $\epsilon_0, \epsilon_1$ are both independent "noise" variables following $\text{Gumbel}(0,1)$ distribution; you can check that this way $y$ follows Bernoulli with $\mathbb{P}(y=1|x)= 1/(1+e^{-\theta^T x})$ as desired.
This way of sampling from a categorical (in this case Bernoulli) distribution is widely known as the Gumbel-max trick in machine learning: https://lips.cs.princeton.edu/the-gumbel-max-trick-for-discrete-distributions/   (The basic idea comes from the reparameterization trick. There's also a closely related Gumbel-softmax trick that essentially turns the above $\arg \max$ operation of Gumbel-max differentiable).
