permb <- function(xu,yu,xp, yp,r, mu=NULL, method="EH",alternative ="two.sided")

{

np=5000

method <- pmatch(method, c("EH", "Maritz"))
    if (is.na(method)) 
        stop("\nbw.method MUST be one of: 'EH', 'Maritz'\n")

 if (missing(r)) {
        r <- cor(xp,yp)
    }

else {
        if ((length(xp) < 1) |(length(yp) < 1)) {
            stop("\n there no paired observations\n")
        }
    }

alt <- pmatch(alternative, c("two.sided","less", "greater"))



if ( method==1)
{

# this function computes the mean of the differences  for the
# unpaired sample


# This part compute permutation test for paired data
paired <- function(x1,x2)

{


n=length(x1)
y1=numeric(n)
y2=numeric(n)



for(i in 1:n)
 
 { u=runif(1)

  if ((u-0.5)<=0) {
 
     
   y1[i]=x1[i]
   y2[i]=x2[i]
  
   } else {
 
 if((u-0.5)>0) { 
 
   y1[i]=x2[i]
   y2[i]=x1[i] 
  }
}
}

axp=mean(y1-y2)



return(axp)

}





unpaired <- function(x1,x2)

{

   m1=length(x1)
   m2=length(x2)


   tc=c(x1,x2)
   n=length(tc)

 
   y=sample(tc,replace=F)
   ax=mean(y[1:m1])
   ay=mean(y[(m1+1):n]) 
  
# the difference between the sample means
  cmean=ax-ay

  return(cmean) 

}

# weight for equal variance


permt <- function(xp,yp,xu,yu,r)

{
n=length(xp)
n1=length(xu)
n2=length(yu)
w=(1/n2+1/n1)/((2-2*r)/n +(1/n2+1/n1))
if(n==0) w=0
if((n1==0)&&(n2==0)) w=1


Tp=paired(xp,yp)

Tu=unpaired(xu,yu)

t10=w*(Tp) +(1-w)*(Tu)

return(t10)
}


t10=replicate(np,permt(xp,yp,xu,yu,r)) 

w=(1/n2+1/n1)/((2-2*r)/n +(1/n2+1/n1))
if(n==0) w=0

if((n1==0)&&(n2==0)) w=1


# the value of the observed test statistic

Tobs=w*(mean(xp)-mean(yp)) +(1-w)*(mean(xu)-mean(yu))

result=c(t10,Tobs)



#The p-value is the chance of seeing a test statistic at least as 
#large as what was observed (.667). This can be #estimated by the
# proportion of permutation statistics that were as 
#large or larger than the observed one

if(alt==1)  { pvalue=sum(abs(t10)>=abs(Tobs))/length(result)

 test="alternative hypothesis: true difference in means is not equal to 0"         

}   

if(alt==2)  { pvalue=sum(t10<=Tobs)/length(result)

 test="alternative hypothesis: true difference in means is less than 0"         

}   

if(alt==3)  { pvalue=sum(t10>=Tobs)/length(result)

 test="alternative hypothesis: true difference in means is greater than 0"         

}   


approach ="Permutation based test on Einsporn and Habtzghi's approach  "


ans=list(approach, test,"pvalue",pvalue)


}


 if(method==2)
{


#  Maritz's Permutation test


permm <- function(xp,yp,xu,yu)

{
n=length(xp)
n1=length(xu)
n2=length(yu)


# 
x=c(xu, rep(0,n2),xp)
y=c(rep(0,n1),yu,yp)




m=n1+n2+n

y1=numeric(m)
y2=numeric(m)

for(i in 1:m)

{
T2=sample(c(x[i],y[i]),replace=F)
y1[i]=T2[1]
y2[i]=T2[2]

}


b=length(y1[y1!=0])
b2=length(y2[y2!=0])

xbar=sum(y1)/b
ybar=sum(y2)/b2



diff=xbar-ybar


return(diff)

}





#Maritz


tm1=replicate(np,permm(xp,yp,xu,yu))


# this part computes the observed value of Marith

n=length(xp)
n1=length(xu)
n2=length(yu)

xl=c(xu, rep(0,n2),xp)
yl=c(rep(0,n1),yu,yp)
bl=c(rep(1,n1),rep(0,n2),rep(1,n))


#ayp <- numeric(np)

xbar1=sum(xl*bl)/sum(bl)
ybar1=(sum(xl)+sum(yl)-sum(xl*bl))/(n1+n2+2*n-sum(bl))


Dm0=xbar1-ybar1



resultm=c(tm1,Dm0)


if(alt==1)  { pvalue=sum(abs(tm1)>=abs(Dm0))/length(resultm)

 test="alternative hypothesis: true difference in means is not equal to 0"         

}   

if(alt==2)  {pvalue=sum(tm1<=Dm0)/length(resultm)

 test="alternative hypothesis: true difference in means is less than 0"         

}   

if(alt==3)  { pvalue=sum(tm1>=Dm0)/length(resultm)

 test="alternative hypothesis: true difference in means is greater than 0"         

}   


approach="Permutation based test on Maritz's approach "
ans=list(approach, test,"pvalue",pvalue)


}


# the value of the test stat using the permutation test

ans


}



# Example
n=20
n1=15
n2=10
r=0.8

xp=rnorm(n)
yp=r*xp+(1-r)*(rnorm(n))

xu=rnorm(n1)
yu=rnorm(n2)
r=0.2
mu=0

permb(xu,yu,xp, yp,r, mu, method="EH", alternative="two.sided")