PowerSignRank = function(alpha=.05,N=10,Delta=1,scale=1,
                       dist="norm",df=1,Nsim=10000){
# PowerSignRank computes the power of the signed-rank Wilcoxon test
# for significance level alpha, sample size N, shift parameter Delta,
# using Nsim simulated samples of size N from a distribution with 
# random sample generator function "rdist" or "runif(N,0,2)-1" when "dist"="unif".
# The shift parameter is added to the sample after it has been multiplied by scale.
# It also computes the corresponding two normal approximations.
# In this function it is assumed that large values of the Wilcoxon signed rank test
# statistic are significant.
# Here "dist" can be:
# "norm" : sample from the standard normal distribution,
# "t"    : sample from the Student t-distribution with df degrees of freedom,
# "logis": sample from the standard logistic distribution
# "unif" : sample from a uniform(-1,1) distribution
# "dexp" : sample from a standard double exponential distribution (internally defined).
#  When scale is different from one each such sample is multiplied by scale.
#  Each generated sample (after possible scaling) is shifted by adding Delta.
#
# This function produces a list of 5 vectors. The first contains the value for dist.
# The second contains the value for input.n=c(N,Nsim), the third contains the 
# input.parms=c(df,scale,Delta) where df is only used when dist="t". 
# The fourth vector contains pvec=c(p,p1,p2,meanVs,varVs) which are used in the 
# first normal approximation. The fifth vector contains 
# power.vec=c(alpha,alpha.c,c.alpha,power.sim,power.norm1,power.norm2),
# where alpha=nominal significance level, alpha.c= achieved level using
# c.alpha=critical level,power.sim=simulated power of the Wilcoxon signed rank test,
# and power.norm1 and power.norm2 are the two normal approximations. In power.norm1
# we use either the exact values of p, p1, p2 (when dist="norm") or we use
# simulated estimates for them for the other distributions.
# 
#================================================================================
rdexp = function(n){
   u=runif(n)
   sel1=u<=.5
   sel2=u>.5
   x=rep(0,length(u))
   x[sel1] = log(2*u[sel1])
   x[sel2] = -log(2*(1-u[sel2]))
   x
}
if(dist=="norm"){
   ell0=1/(scale*sqrt(2*pi))
   ell0.star=1/(scale*2*sqrt(pi))
}
if(dist=="unif"){
   ell0=1/(2*scale)
   ell0.star=ell0
}
if(dist=="logis"){
   ell0=1/(scale*4)
   ell0.star=1/(scale*6)
}
if(dist=="dexp"){
   ell0=1/(scale*2)
   ell0.star=1/(scale*4)
} 
if(dist=="t"){
   ell0=gamma((df+1)/2)/(gamma(df/2)*sqrt(pi*df)*scale)
   ell0.star=(gamma((df+1)/2)/gamma(df/2))^2 *
           gamma(df+1/2)/gamma(df+1) / (scale*sqrt(df*pi))  
} 
k=qsignrank(1-alpha,N)
c.alpha=k+1
if(dist=="unif"){Zmat=matrix(runif(Nsim*N,0,2)-1,ncol=Nsim)}
if(dist=="t"){Zmat=matrix(rt(Nsim*N,df,0),ncol=Nsim)}
if(dist!="t" & dist!="unif"){ Zmat=matrix(get(paste("r",dist,sep=""))(Nsim*N),ncol=Nsim)}
if(scale!=1){Zmat=Zmat*scale+Delta}else{Zmat=Zmat+Delta}
if(dist=="norm"){
  if(!exists("pmnorm")) library(mnormt)
  p=pnorm(Delta/scale)
  p1=pnorm(sqrt(2)*Delta/scale)
  p2=pmnorm(c(sqrt(2)*Delta/scale,sqrt(2)*Delta/scale),
                c(0,0),varcov=matrix(c(1,.5,.5,1),ncol=2))
}else{
  p2=mean(Zmat[1,]+Zmat[2,]>0 & Zmat[1,]+Zmat[3,]>0)
  p1=mean(Zmat[1,]+Zmat[2,]>0)
  p=mean(Zmat[1,]>0)
}
if(Delta==0){
  p  =.5
  p1 =.5
  p2 = 1/3
}
meanVs=N*(N-1)*p1/2+N*p
varVs=N*(N-1)*(N-2)*(p2-p1^2)+(N*(N-1)/2)*(2*(p-p1)^2+3*p1*(1-p1))+N*p*(1-p)
meanVs0=N*(N+1)/4
varVs0=N*(N+1)*(2*N+1)/24
Rmat=apply(abs(Zmat),2,rank)
RZmat=Rmat*(Zmat>0)
Vs=apply(RZmat,2,sum)
alpha.c=1-psignrank(k,N)
power.sim=mean(Vs>=c.alpha)
power.norm1=pnorm((meanVs-c.alpha+.5)/sqrt(varVs))
power.norm2=pnorm((N*(N-1)*ell0.star+N*ell0)*Delta/sqrt(N*(N+1)*(2*N+1)/24)-
   qnorm(1-alpha.c))
input.N=c(N,Nsim)
input.parms=c(df,scale,Delta)
names(input.N)=c("N","Nsim")
names(input.parms)=c("df","scale","Delta")
pvec=c(p,p1,p2,meanVs,varVs)
names(pvec)=c("p","p1","p2","meanVs","varVs")
power.vec=c(alpha,alpha.c,c.alpha,power.sim,power.norm1,power.norm2)
names(power.vec)=c("alpha","alpha.c","c.alpha","power.sim","power.norm1","power.norm2")
list(dist=dist,input.N=input.N,input.parms=input.parms,pvec=pvec,power.vec=power.vec)

}