# Function "Smirnov2x2.test"  
# plots extreme and rejection regions (one-sided, nominal alpha=0.05)
# for the test with the Smirnov type test statistic
# calculating p-value and the actual size
# for a pair of 2x2 contingency tables 
#  x1,x2	y1,y2
#  x3,x4	y3,y4
# Author: Anastasia Ivanova <aivanova@bios.unc.edu>

Smirnov2x2.test_function(x1,x2,x3,x4,y1,y2,y3,y4)
{
xx1_min(x1+x2,x1+x3)
yy1_min(y1+y2,y1+y3)
v1_rep(c(0:xx1),yy1+1)
v2_as.vector(t(matrix(rep(c(0:yy1),xx1+1),ncol=(xx1+1))))
nprob_dhyper(v1,x1+x2,x3+x4,x1+x3)*dhyper(v2,y1+y2,y3+y4,y1+y3)
Gamma1_cbind(v1,v2,nprob)
observed_as.numeric(Gamma1[,1] == x1 & Gamma1[,2] == y1)
Gamma1_cbind(Gamma1,observed)
x_v1
y_v2

par(oma=c(0, 0, 3, 0),mfcol=c(1,2))
minx_min(Gamma1[,1])
maxx_max(Gamma1[,1])
miny_min(Gamma1[,2])
maxy_max(Gamma1[,2])

Smirnov_pmax(v1/(x1+x2)-(x1+x3-v1)/(x3+x4),v2/(y1+y2)-(y1+y3-v2)/(y3+y4),0)
rankSmirnov_rank(Smirnov)
Gamma1_cbind(Gamma1,Smirnov,rankSmirnov)
Smirnovobs_pmax(x1/(x1+x2)-(x1+x3-x1)/(x3+x4),y1/(y1+y2)-(y1+y3-y1)/(y3+y4),0)
Smirnov.region_Gamma1[Smirnov >= Smirnovobs ,]
noSmirnov.region_Gamma1[Smirnov < Smirnovobs ,]
Smirnov.p_sum(Smirnov.region[,3])
plot(Gamma1[,1],Gamma1[,2],xlab="X11",ylab="Y11",type="n",lab=c(maxx-minx+1,maxy-miny+1,2))
title(main=paste("Smirnov Test","(p=",round(Smirnov.p,digits=4),")"))
text(Smirnov.region[,1],Smirnov.region[,2],round(Smirnov.region[,6],digits=0))
if (is.array(noSmirnov.region)) points(noSmirnov.region[ ,1],noSmirnov.region[ ,2], type="p",pch=1) else
    if (length(noSmirnov.region)>0) points(noSmirnov.region[1],noSmirnov.region[2], type="p",pch=1) 
#points(noSmirnov.region[,1],noSmirnov.region[,2],pch=1)
points(x1,y1, type="p", pch=1, cex=2)

size_0
porog_max(Smirnov)
Smirnov1_Smirnov
while (size <= 0.05) 
{ porog_max(Smirnov1) 
  size_size+sum(Gamma1[Smirnov == porog,3])
  Smirnov1_Smirnov1[Smirnov1 != porog]
}
Smirnov.rejreg_Gamma1[Smirnov > porog,]
noSmirnov.rejreg_Gamma1[Smirnov <= porog,]
if (is.array(Smirnov.rejreg)) Smirnov.size_sum(Smirnov.rejreg[ ,3]) else if (length(Smirnov.rejreg)>0) Smirnov.size_sum(Smirnov.rejreg[3]) else Smirnov.size_0

plot(Gamma1[,1],Gamma1[,2],xlab="X11",ylab="Y11",type="n",lab=c(maxx-minx+1,maxy-miny+1,2))
title(main=paste("Smirnov Test","(size=",round(Smirnov.size,digits=4),")"))
if (is.array(Smirnov.rejreg)) text(Smirnov.rejreg[,1],Smirnov.rejreg[,2],round(Smirnov.rejreg[,6],digits=0)) else
    if (length(Smirnov.rejreg)>0) text(Smirnov.rejreg[1],Smirnov.rejreg[2],round(Smirnov.rejreg[6],digits=0))
if (is.array(noSmirnov.rejreg)) points(noSmirnov.rejreg[ ,1],noSmirnov.rejreg[ ,2], type="p",pch=1) else
    if (length(noSmirnov.rejreg)>0) points(noSmirnov.rejreg[1],noSmirnov.rejreg[2], type="p",pch=1) 

points(x1,y1, type="p", pch=1, cex=2)

mtext(outer=T,paste("Smirnov Test for (",x1,",",x2,",",x3,",",x4,");(",y1,",",y2,",",y3,",",y4,")"),side=3,cex=1.3)
c(Smirnov.p,Smirnov.size)
}

Smirnov2x2.test(4,14,0,19,11,15,1,25)