Expected values normal order statistics

Question

Given this equation pertaining to algorithm AS 177 (from Royston, 1982):

In a sample of size $n$ the expected value of the $r$th largest order statistic is given by $$E(r,\,n)=\frac{n!}{(r-1)!(n-r)!}\int_{-\infty}^\infty x\{1-\Phi(x)\}^{r-1}\{\Phi(x)\}^{n-r}\phi(x)dx,$$ where $\phi(x)=\frac{1}{\sqrt{2\pi}}\exp(-\tfrac12x^2)$ and $\Phi(x)=\int_{-\infty}^z\phi(z)dz$.

I have been trying to code this (really the expected normal order statistic of the Terry-Hoeffding Test) into Python with some difficulty. I can code everything up to the integral and I can also get the CDF and PDF functions. Where I am confused is the $x$ (the actual value?) before the curly brackets and the $dx$ that follows the PDF.
Knowing little of mathematical notation, I am guessing I multiply the first fraction by $x$ times everything in the 2 sets of curly brackets and PDF but that doesn't seem correct. Also everything after "where" is Greek to me (literally and figuratively).

Any help in helping me understand this equation so that I can convert it to code would be greatly appreciated.

Answer 1

This should do it, once you install scipy:

from scipy.integrate import quad
from scipy.special import binom
from scipy.stats import norm

inf, phi, Phi = float('inf'), norm.pdf, norm.cdf

def E(r, n):
    def f(x):   
        F = Phi(x)
        return x*(1-F)**(r-1)*F**(n-r)*phi(x)
    return r*binom(n, r)*quad(f, -inf, inf)[0]

See here, here & here.