function (n = 6, m = 4,a=6,PDF=F) 
{
if(PDF==T) pdf(file=paste("contcorrn",n,"m",m,"a",a,".pdf",sep=""),width=11)
    N = m + n
    mean.Ws = n * (N + 1)/2 - n * (n + 1)/2
    var.Ws = m * n * (N + 1)/12
    x = 0:(m * n)
    x1=0:a
    probs = dwilcox(x, n, m)
    mp=max(probs)
    plot(x, probs, type = "h", lty = 2,xlim=c(-5,m*n+5),ylim=c(0,mp),
      xlab="a",ylab=expression(P(W[XY]==a)),axes=F)
axis(1)
axis(2)
    abline(h=0)
    points(x, probs, pch = 16, cex = 0.7)
    rect(x - 0.5, rep(0, length(x)), x + 0.5, probs)
    rect(x1-.5,rep(0,length(x1)),x1+.5,probs[1:length(x1)],density=10)
    rect(6,0,6.5,probs[length(x1)],density=10,angle=-45)
    z = seq(-5, m * n + 5, length.out = 200)
    phiz = dnorm(z, mean.Ws, sqrt(var.Ws))
    lines(z, phiz, col = "blue")
    text(a-.2,.995*mp,substitute("true  "~P(W[XY]<= ax)==px,list(ax=a,
             px=signif(sum(probs[1:length(x1)]),4))),adj=1)
    abline(v=a,lty=2)
    abline(v=a+.5)
    p1=pnorm((a-mean.Ws)/sqrt(var.Ws))
    p2=pnorm((a+.5-mean.Ws)/sqrt(var.Ws))
    text(a-.2,.78*mp,substitute(P(W[XY]<= ax)%~~%px,list(ax=a,
             px=signif(p1,4))),adj=1)
text(a-.2,.6*mp,substitute(P(W[XY]<= ax)%~~%px,list(ax=a,
             px=signif(p2,4))),adj=1)
text(a-.2,.65*mp,"with continuity correction",adj=1,cex=.9)
text(a-.2,.83*mp,"without continuity correction",adj=1,cex=.9)
text(a-.2,.9*mp,"normal approximation",adj=1)
arrows(a-5,.57*mp,a+.5-.02,.57*mp,angle=15,length=.1)
    arrows(a-6,.75*mp,a-.02,.75*mp,angle=15,length=.1)
text(m*n,mp*.95,paste("n =",n," ,  m =",m),adj=1)
if(PDF==T) dev.off()
}