As far as I know, the sigmoid is often used as the activation function of the output layer mainly because it is a convenient way of producing an output $p \in [0, 1]$, which can be interpreted as a probability, although that can be misleading or even wrong (if you interpret it as an uncertainty too).
You may require the output of the neural network to be a probability, for example, if you use a cross-entropy loss function, although you could in principle produce only $0$s or $1$s. The probability $p$ can then be used to decide the class (or label) of the input. For example, if $p > \alpha$, then you purposedly decide that the input belongs to class $1$, otherwise, it belongs to class $0$. The parameter $\alpha$ is called the classification (or decision) threshold. The choice of this threshold can actually depend on the problem and it is one of the reasons people use the AUC metric, i.e. to avoid choosing this classification threshold.
Can neural networks with a sigmoid as the activation function of the output layer approximate continuous functions? Is there an analogue to the universal approximation theorem for this case?
The most famous universal approximation theorem for neural networks assumes that the activation functions of the units of the only hidden layer are sigmoids, but it does not assume that the output of the network will be squashed to the range $[0, 1]$. To be more precise, the UAT (theorem 2 of Approximation by Superpositions of a Sigmoidal Function, 1989, by G. Cybenko) states
Let $\sigma$ be any continuous sigmoidal function. Then finite sums of the
G(x)=N∑j=1αjσ(yTjx+θj)
are dense in $C(I_n)$.
In other words, given any $f \in C(I_n)$ and $\epsilon > 0$, there is a sum, $G(x)$, of the above form, for which
|G(x)−f(x)|<ϵ
Here, $f$ is the continuous function that you want to approximate, $G(x)$ is a linear combination of the outputs of $N$ (which should be arbitrarily big) units of the only hidden layer, $I_n$ denotes the $n$-dimensional unit cube, $[0, 1]^n$, $C(I_n)$ denotes the space of continuous functions on $I_n$, $x \in I_n$ (so the assumption is that the input to the neural network is an element of $[0, 1]^n$, i.e. a vector $x \in \mathbb{R}^n$, whose entries are between $0$ and $1$) and $y_j$ and $\theta_j$ are respectively the weights and bias of the $j$ unit. The assumption that $f$ is a real-valued function means that $f$ can take any value on $\mathbb{R}$ (i.e. $f: [0, 1]^n \rightarrow \mathbb{R}$). You should note that $G(x)$ is the output of the neural network, which is a combination (where the coefficients are $\alpha_j$) of the outputs of the units in the only hidden layer, so there's no restriction on the output of $G(x)$, unless you restrict $\alpha_i$ (but, in this theorem, there's no restriction on the values $\alpha_j$ can take).
Of course, if you restrict the output of the neural networks to the range $[0, 1]$, you cannot approximate all continuous functions of the form $f: [0, 1]^n \rightarrow \mathbb{R}$ (because not all of these functions will have the codomain $[0, 1]$)! However, the sigmoid has an inverse function, i.e. the logit, so you can reverse the output of such a neural network. So, in this sense (i.e. by reversing the output of the sigmoid), a neural network with a sigmoid as the activation function of the output layer can potentially approximate any continuous function too.
The UAT above only states the existence of $G(x)$ (i.e. it's an existence theorem). It doesn't tell you how you can find $G(x)$. So, if you use a sigmoid as the activation function of the output layer or not is a little bit orthogonal to the universality of neural networks.
