function (x = c(5, 7, 16, 12), 
          y = c(7, 9, 15, 9),Nsim=1000){
    d = x + y; n = sum(x); m = sum(y); ell = length(x)
    N = m + n
    midrank = cumsum(d) - (d - 1)/2
    midrank.vec=NULL
    for(j in 1:ell){
      midrank.vec=c(midrank.vec,rep(midrank[j],d[j]))}
    Ws = sum(x * midrank)
    Ws.vec=NULL
    for(i in 1:Nsim){
        Ws.vec[i]=sum(sample(midrank.vec,n,replace=F))
    }
    pval=mean(Ws.vec>=Ws)
    out = c(Ws,  pval)
    names(out) = c("Ws.star", "pval.Ws.star")
    out
}