If p = G(x_i \beta), then Wolfram alpha says
\int_0^1 p^y (1 - p)^{1 - y} dy = \frac{1 - 2 p}{2 \text{tan}^{-1}(1 - 2 p)}
So since we can normalize the distribution over y, then it’s just a regular likelihood and it will hopefully sample in Stan.
I think there should be another exp in the denominator of G.