function (x = c(5, 7, 16, 12), y = c(7, 9, 15, 9)) 
{
    d = x + y
    n = sum(x)
    m = sum(y)
    N = m + n
    midrank = cumsum(d) - (d - 1)/2
    Ws = sum(x * midrank)
    meanWs = n * (N + 1)/2
    varWs = m * n * (N + 1)/12 - m * n * sum(d^3 - d)/(12 * N * 
        (N - 1))
    pval = 1 - pnorm((Ws - meanWs)/sqrt(varWs))
    out = c(Ws, meanWs, sqrt(varWs), pval)
    names(out) = c("Ws.star", "mean.Ws.star", "sdev.Ws.star", 
        "pval.Ws.star")
    out
}