## Programed by Shang-Yi Lin &  Tsung-Chen Hsieh
## Phd student, Institute of Statistics, National Tsing Hua University
## Date: 2011, August
## A Subroutine for analyzing data


#Function of plot confidence band.
conf.reg=function(x,LCL,UCL,...){
	polygon(c(x,rev(x)),c(LCL,rev(UCL)),...)
}

#Function of calculating the Chao2 estimator.
EstS=function(x,T){
	D=length(which(x!=0))
	Q1=sum(x==1)
	Q2=sum(x==2)
	if(Q2==0){out=D+((T-1)/T)*Q1*(Q1-1)/(2)}
	else {out=D+((T-1)/T)*(Q1^2)/(2*Q2)}
return(out)
}

#Function of summary input data.
summary.dat=function(dat,T){
	D=apply(as.matrix(dat),2,function(x)sum(x>0))
	Q1=apply(as.matrix(dat),2,function(x)sum(x==1))
	Q2=apply(as.matrix(dat),2,function(x)sum(x==2))
	Sest=rep(0,ncol(dat))
	for(i in 1:ncol(dat)){
		Sest[i]=EstS(dat[,i],T[i])
	}
	out=cbind(D,Q1,Q2,Sest)
	colnames(out)=c("D","Q1","Q2","Chao2.est")
	as.data.frame(out)
}

pro<-function(A=A,T=T,k=k)
{
	prod1=1
	temp1=0
	temp2=0
	i=c(1:k)
	temp1=T-A-i+1
	temp2=T-i+1
	prod1=prod(temp1/temp2)
return(prod1)
}

#Function for rarefaction and extropolation 
#Y means data, which is a vector.
#T means the quadratic numbers.
#t means rarified quadratic sequence.
#tt means the extrapolation quadratic samples. 
RareExtro=function(Y,T,t,tt=100){
	D=sum(Y>0)
	Q=rep(0,T)
	Q=sapply(c(1:T),function(j){sum(Y==j)})

	Q0.hat=0
	if(Q[2]>0)  { Q0.hat=((T-1)/T)*(Q[1]^2)/(2*Q[2]) }
	if(Q[2]==0) { Q0.hat=((T-1)/T)*(Q[1]*(Q[1]-1))/2 }

	####Rarefaction part####
	Fun=function(ss){
		Sub=function(x){
			if((T-x)>=ss) {aa=pro(x,T,ss)}
			if((T-x)<ss) {aa=0}
		return(aa)
		}
		a=sapply(Y,Sub)
		S.rare=sum(1-a)
		V.rare=sum((1-a)^2)-S.rare^2/(D+Q0.hat)
		c(ss,S.rare,sqrt(V.rare),S.rare-1.96*sqrt(V.rare),S.rare+1.96*sqrt(V.rare))
	}
	out.rare=t(sapply(t,Fun))


	####Extrapolation part####
	S.exp=rep(0,tt)
	S.exp=Q0.hat*(1-(1-Q[1]/(Q[1]+T*Q0.hat))^(c(1:tt)))
	g.S.exp=0
	g.S.exp=D+S.exp

	##Variance estimator of Extrapolation part
	V=matrix(0,T,T)
	for(i in 1:(T-1)){
		for(j in (i+1):T){
			V[i,j]=-Q[i]*Q[j]/(D+Q0.hat)
		}
	}
	V=V+t(V)
	temp=0
	temp=Q*(1-Q/(D+Q0.hat))
	diag(V)=c(temp)
	Varex<-function(zz)
	{
		dif=rep(0,T)
		A=0
		A=1-( 1 - Q[1]/ (Q[1]+T*Q0.hat) )^zz
		B=0
		B=( 1-Q[1]/(Q[1]+T*Q0.hat) )^(zz-1)
		C=0
		C=(T-1)/(  2*Q[2]* (( 1+(T-1)*Q[1]/(2*Q[2]) )^2)   )
		D=0
		D=(T-1)*(-1*Q[1])/(  2 * (Q[2]^2) * ( ( 1+(T-1)*Q[1]/2/Q[2] )^2 )  )
		dif[1]=1+((T-1)*Q[1]/T/Q[2])*A+Q0.hat*(-zz)*B*C
		dif[2]=1+((T-1)*(-1*Q[1]^2)/T/2/(Q[2]^2))*A+Q0.hat*(-zz)*B*D
		dif[3:T]=1
		var=0
		var=dif%*%V%*%dif
	return(var)
	}
	var.exp=rep(0,tt)
	var.exp=sapply(c(1:tt),Varex)
	out.exp=cbind(c(1:tt)+T,g.S.exp,sqrt(var.exp),g.S.exp-1.96*sqrt(var.exp),g.S.exp+1.96*sqrt(var.exp))

	##Combine rarefaction and extrapolation
	out=rbind(out.rare,out.exp)
	colnames(out)=c("Quadrat","S.hat","sd","95%LCL","95%UCL")
	return(as.data.frame(out))
}