function(Upper=T,sim=T,Nsim=10000){
xPsych=c(-14, 3, 1, -16, -21, 7, -7, -13, -22, -17, -14, -8, 7, -18, 
-13, -18, -9, -22, -25, -24, -18, -13, -13, -18, -5)
yPsych=c(-18, -16, -9, -14, 3, -9, -16, 10, -11, -3, -13, -21, -2, -11, 
-16, -12, -13, -6, -9, -7, -11, -9)
x=xPsych
y=yPsych
z=c(x,y)
rz=rank(z)
m=length(x)
n=length(y)
N=m+n
R=sort(rz[1:m])
S=sort(rz[(m+1):(m+n)])
Ws.star=sum(S)
EWs.star=n*(N+1)/2
varWs.star=n*(N-n)*var(rz)/N
p.val.sim=NULL
if(sim==T){
temp=NULL
   for(i in 1:Nsim){
      temp[i]=sum(sample(rz,n,replace=F))
}
if(Upper==T){
p.val.sim=mean(temp>=Ws.star)}else{
p.val.sim=mean(temp<=Ws.star)
}
}
if(Upper==T){
p.val=1-pnorm((Ws.star-EWs.star)/sqrt(varWs.star))}else{
p.val=pnorm((Ws.star-EWs.star)/sqrt(varWs.star))
}
out=signif(c(Ws.star,EWs.star,varWs.star,p.val,p.val.sim),4)
names(out)=c("Ws.star","EWs.star","varWs.star","p.val","p.val.sim")
out
}