R/pblocks.R

In rkt: Mann-Kendall Test, Seasonal and Regional Kendall Tests

pblocks <- function(BL,S){

# dim(BL) is always 6 decreasing order

S<-abs(S)
NL<-BL*(BL-1)/2+1
N<-BL[1]
N<-N*(N-1)/2+1
BS<-rep(0,N*6)
dim(BS)<-c(N,6)
BP<-rep(1,N*6)
dim(BP)<-c(N,6)
for (B in 1:6)
	{
	for (I in 1:(BL[B]*(BL[B]-1)/2+1)) BS[I,B]<-BL[B]*(BL[B]-1)/2-(I-1)*2
	}
BP[1,1]<-tauprob(BL[1],BS[1,1])
for (I in 2:(BL[1]*(BL[1]-1)/2+1)) 
	{
	if (BS[I,1]>=0)
		{
		BP[I,1]<-tauprob(BL[1],BS[I,1])-tauprob(BL[1],BS[I-1,1])
		}	
	else
		{
		BP[I,1]<-BP[BS[,1]==-BS[I,1],1]
		}
	}

for (B in 2:6)
	{
	if (BL[B]>0)
		{
		if (BL[B]==BL[B-1])
			{
			BP[,B]<-BP[,B-1]
			}
			else
			{
			BP[1,B]<-tauprob(BL[B],BS[1,B])
			for (I in 2:(BL[B]*(BL[B]-1)/2+1)) 
				{
				if (BS[I,B]>=0)
					{
					BP[I,B]<-tauprob(BL[B],BS[I,B])-tauprob(BL[B],BS[I-1,B])
					}	
				else
					{
					BP[I,B]<-BP[BS[,B]==-BS[I,B],B]
					}
				}
			}
		}
	}
PTOT<-0
for (I1 in 1:NL[1])
	{
for (I2 in 1:NL[2])
	{
for (I3 in 1:NL[3])
	{
for (I4 in 1:NL[4])
	{
for (I5 in 1:NL[5])
	{
for (I6 in 1:NL[6])
	{
	thisS<-BS[I1,1]+BS[I2,2]+BS[I3,3]+BS[I4,4]+BS[I5,5]+BS[I6,6]
	if (thisS>=S) PTOT<-PTOT+BP[I1,1]*BP[I2,2]*BP[I3,3]*BP[I4,4]*BP[I5,5]*BP[I6,6]
	}
	}
	}
	}
	}
	}
return(PTOT)
}