ustat <- function(x,y){
x <- x[!is.na(x)]
y <- y[!is.na(y)]
m <- length(x)
n <- length(y)
P <- rep(NA,m)
Q <- rep(NA,n)
for (i in 1:m){
P[i] <- sum(y < x[i])
}
for (j in 1:n){
Q[j] <- sum(x < y[j])
}
Pbar <- mean(P)
Qbar <- mean(Q)
V1 <- sum((P-Pbar)^2)
V2 <- sum((Q-Qbar)^2)
Uhat <- (n*Qbar - m*Pbar)/(2*sqrt(V1 + V2 + Pbar*Qbar))
Uhat
}