Solve integral of $\Phi(a+bx)^k(1-\Phi(a+bx))^{n-k}\phi(x)$

Question

I want to know if there is a way of solving the integral $$ \int_{\mathbb{R}} \Phi(a+bx)^k(1-\Phi(a+bx))^{n-k}\phi(x)dx $$ or approximate numerically it with a given error. The closest question (answered) I found is gaussian integral of power of cdf : $\int_{-\infty}^{+\infty} \Phi(x)^n \cdot \phi(a+bx) \cdot dx$ which gives a very nice and complete result, and also some approximations. The logic behind this integral is that I have a r.v. $K|Z\sim Binomial(n,p(Z))$ where $Z\sim N(0,1)$, with $$ p(Z) = \Phi(a+bZ). $$ I want to know the unconditional distribution of $K$ through $$ \begin{aligned} f(k) &= \int_{\mathbb{R}} f(k|z)f(z)dz\\ & = \int_{\mathbb{R}} \binom{n}{k}p(z)^k(1-p(z))^{n-k}\phi(z)dz \end{aligned} $$ with $p(z) = \Phi(a+bz)$. I have an idea of how to solve it for $n\to \infty$, but I would like to know if there are other results. Thank you!

Answer 1

Using R and the default precision (which can be tweaked) here is a numerical approach to find the probability mass function given $a$, $b$, and $n$.

distribution <- function(a, b, n) {
Define integrand function
integrand <- function(x, a, b, n, k) {
    p <- pnorm(a+bx)
    p^k  (1-p)^(n-k) * exp(-x^2/2)/sqrt(2*pi)
  }
Find probabilities for k=0,1,...,n
prob <- rep(NA, n + 1)
  abs.error <- rep(NA, n + 1)
  for (i in 0:n) {
      temp <- integrate(integrand, -Inf, Inf, a=a, b=b, n=n, k=i)
      prob[i+1] <- choose(n, i)*temp$value
      abs.error[i+1] <- temp$abs.error 
  }
  data.frame(k=0:n, prob=prob, abs.error=abs.error)
}
distribution(a, b, n)
distribution(1, 3, 5)
  k       prob    abs.error
1 0 0.24088017 5.692414e-05
2 1 0.07123875 6.393196e-06
3 2 0.05928153 5.117496e-05
4 3 0.06283731 3.721303e-05
5 4 0.08669902 7.994856e-06
6 5 0.47906321 4.687144e-05

Below are two symbolic results. The first is a closed-form for the mean of the unconditional distribution of $K$ (which involves switching the order of summation and integration):

$$\mu=\sum_{k=0}^n k \int_{-\infty}^\infty \binom{n}{k}\Phi(a+b z)^k (1-\Phi(a+b z))^{n-k}\frac{e^{-z^2/2}}{\sqrt{2\pi}}dz$$

$$=\int_{-\infty}^\infty \sum_{k=0}^n k \binom{n}{k}\Phi(a+b z)^k (1-\Phi(a+b z))^{n-k}\frac{e^{-z^2/2}}{\sqrt{2\pi}}dz$$

$$=\int_{-\infty}^\infty n \Phi(a+b z)\frac{e^{-z^2/2}}{\sqrt{2\pi}}dz$$ $$=n\Phi\left(\frac{a}{\sqrt{1+b^2}}\right)$$

The variance of $K$ is found (also using the list of gaussian integrals)

$$\sigma^2=n \left(n \left(1-\Phi \left(\frac{a}{\sqrt{b^2+1}}\right)\right) \Phi \left(\frac{a}{\sqrt{b^2+1}}\right)-2 (n-1) T\left(\frac{a}{\sqrt{b^2+1}},\frac{1}{\sqrt{2 b^2+1}}\right)\right)$$

where $\Phi(.)$ is the standard normal cumulative distribution function and $T(.,.)$ is Owen's T function.