In Gaussian Process binary classification, why are sigmoid functions preferred over Gaussian functions? I am currently studying "Gaussian Processes for Machine Learning", and in chapter 3 they state that the posterior $p(y_*|X,\mathbf{y},\mathbf{x}_*)$ (eq. 3.10) and the latent variable posterior $p(f_*|X,\mathbf{y},\mathbf{x}_*)$ (eq. 3.9) cannot generally be solved analytically, due to the sigmoid likelihoods in (3.9) and the sigmoid function in (3.10). To save people from having to look up the equations, they are as follows:
$$
\begin{align}
p(y_*=+1|X,\mathbf{y},\mathbf{x}_*) &= \int\sigma(f_*)\,p(f_*|X,\mathbf{y},\mathbf{x}_*)\,df_*\quad\quad&\mbox{(3.10)} \\
p(f_*|X,\mathbf{y},\mathbf{x}_*) &= \int p(f_*|X,\mathbf{x}_*,\mathbf{f})\,p(\mathbf{f}|X,\mathbf{y})\,d\mathbf{f}&\mbox{(3.9)}
\end{align}
$$
My main question is: for binary classification with $f$ modelled as a Gaussian Process, why use sigmoid functions at all (in either equation) instead of the Gaussian function
$$
p(y=+1\,|\,f(\mathbf{x}))=g(f(\mathbf{x}))\triangleq\exp\left\{-\frac{f^2(\mathbf{x})}{2}\right\} \enspace?
$$
This would lead to closed-form solutions to both integrals. The Gaussian function is not monotonic, like sigmoid functions, but GPs can generate functions with multiple turning points, so monotonicity seems unnecessary. To ensure that (3.10) converges to $\frac{1}{2}$ when $\mathbf{x_*}$ is far from the training data, it would presumably suffice to give the prior $p(\mathbf{f}|X)$ a mean:
$$
\begin{align}
\mathbb{E}[\mathbf{f}|X] &= \omega\mathbf{1}_n \\
\omega&=\sqrt{-2\ln\frac{1}{2}} \enspace,
\end{align}
$$
where $\mathbf{1}_n$ is a vector of $n$ $1$'s and $n$ is the number of training samples, since:
$$
g\left(\omega\right)=\frac{1}{2}\enspace.
$$
In contrast to the behaviour of sigmoid likelihoods, Gaussian likelihoods would favour large (positive or negative) entries in $\mathbf{f}$ for negatively labelled input points, and small entries in $\mathbf{f}$ for positively labelled points.
Would Gaussian functions lead to problems that do not occur with sigmoids? Are there any papers in which Gaussian functions have been used in binary GP classification instead of sigmoids?
Update, 25th May 2017
On further reflection, the non-zero prior mean suggested above also helps to resolve the ambiguity about what the sign of $f$ should be ($g$ does not favour either sign; $g(f(\mathbf{x}))=g(-f(\mathbf{x}))$ ). Resolving this ambiguity seems to be important, because if the mean of the prior, $p(\mathbf{f}|X)$, was zero, then the mean of $p(\mathbf{f}|X,\mathbf{y})$ would also be zero under a likelihood defined by $g$, as the prior and likelihood would both be even functions of $\mathbf{f}$. I.e.:
$$
\begin{align}
p(\mathbf{y}|\mathbf{f})&=\prod_{i=1}^n p(\mathbf{y}_i|\mathbf{f}_i) \\
p(\mathbf{y}_i|\mathbf{f}_i) &= 
\begin{cases}
g(\mathbf{f}_i) & ,\;\mathbf{y}_i=+1 \\
1-g(\mathbf{f}_i) & ,\;\mathbf{y}_i=-1
\end{cases} \\
\therefore \mathbb{E}[\mathbf{f}|X]=\mathbf{0} \enspace\rightarrow\enspace p(-\mathbf{f}|X,\mathbf{y})
&=\frac{p(\mathbf{y}|-\mathbf{f})p(-\mathbf{f}|X))}{p(\mathbf{y}|X)}
=\frac{p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|X))}{p(\mathbf{y}|X)}
=p(\mathbf{f}|X,\mathbf{y}) \enspace.
\end{align}
$$
If the mean of $p(\mathbf{f}|X,\mathbf{y})$ was zero, the training set labels $\mathbf{y}$ would not provide any information about the query point label $y_*$, so clearly we must not allow this. So in addition to defining $\mathbb{E}[\mathbf{f}|X]=\omega\mathbf{1}_n$, perhaps we should further bias $p(\mathbf{f}|X,\mathbf{y})$ towards positive $\mathbf{f}$ by giving the prior $p(\mathbf{f}|X)$ relatively small standard deviations, e.g. $\sqrt{k(x,x)}=\frac{\omega}{\beta}$, where $k$ is the covariance function and $\beta\in[2,3]$. If we do this, we should probably also scale up $g$'s argument, so that $\mathbf{f}$ will not have to be improbably far from the prior mean to produce small values of $g$:
$$
g(f(\mathbf{x});s)=\exp\left\{-\frac{f^2(\mathbf{x})}{2s^2}\right\}\enspace,
$$
where $s<1$.
Would this be a reasonable way to fix the $f$ sign ambiguity problem?
 A: I believe they mention this in the footnote to chapter 3 (first page)

One may choose to ignore the discreteness of the target values, and use a regression treatment, where all targets happen to be say ±1 for binary classification. This is known as least-squares classification, see section 6.5.

Looking at 6.5 http://www.gaussianprocess.org/gpml/chapters/RW6.pdf they mention the advantage of using sigmoid functions is that the outputs can be be interpreted probabilistically (ie, the probability that an example has a positive response).
A: The problem with this approach is that the number of terms in $p(\mathbf y|\mathbf f)$ would grow exponentially with the number of negatively-labelled points in the training set, so the closed-form solution to (3.9) would have exponential time complexity. More specifically, if we assume, without loss of generality, that
$$
\mathbf y_1=\ldots=\mathbf y_a=-1 \enspace,\enspace \mathbf y_{a+1}=\ldots=\mathbf y_n=+1 \enspace,
$$
then
$$
p(\mathbf y|\mathbf f) 
  = \left(\prod_{i=1}^a (1-g(\mathbf f_i))\right)
    \prod_{i=a+1}^n g(\mathbf f_i) \enspace.
$$
To obtain a closed-form solution to (3.9), we have to expand the first product into a sum of (unnormalised) Gaussian functions, so that we can integrate each one separately:
$$
\prod_{i=1}^a (1-g(\mathbf f_i)) 
  = \sum_{I\in \mathcal{P}\{1,\ldots,a\}} (-1)^{|I|}\exp\left\{
      -\frac{1}{2}\sum_{i\in I}\mathbf f^2_i
    \right\} \enspace.
$$
There are $2^a$ sets in the power set $\mathcal P\{1,\ldots,a\}$ of the negatively-labelled-point indices $\{1,\ldots,a\}$, so solving (3.9) would involve computing $2^a$ Gaussian integrals.
