# Code Program 0
	#rm(list=ls(all=TRUE))
	#options(digits =22)
	library("sampling")
	library("LaplacesDemon")
#-----------------Function---------------------------------------



DefineddataSample<-function(DataP,n,N,p) 
{		
		y=DataP$y		#Study variable
		x=DataP$x		#Auxiliary variable
		k=DataP$k		#Size variable	
		u=DataP$u		#Caribration variable
		w1=DataP$w		# Auxiliary variable use estimated response probability
		
		r=DataP$r
		Sumk=sum(k)

		Pi=k/Sumk
		Phiipop=Pi*(N-n)/(N-1)+(n-1)/(N-1)

		ysample=rep(0,n)
		xsample=rep(0,n)
		ksample=rep(0,n)
		usample=rep(0,n)		
		wsample=rep(0,n)

		Pisample=rep(0,n)
		Phiisample=rep(0,n)
		rsample=rep(0,n)
		
		
		s=UPmidzuno(Phiipop)
		nsample=(1:N)[s==1]	
		ysample=y[nsample]
		xsample=x[nsample]
		ksample=k[nsample]
		usample=u[nsample]
		rsample=r[nsample]		
		wsample=w1[nsample]
		Pisample=Pi[nsample]
		Phiisample=Phiipop[nsample]
	 	datasample1=data.frame(ysample ,xsample,usample, ksample ,wsample ,Pisample, Phiisample ,rsample, nsample)
		return(datasample1)
}

Createpihat<-function(Data)
{
	n=length(Data$rsample)
	r=Data$rsample
	w=Data$wsample


	pihat0=rep(0,n)
		mydata=data.frame(r,w) 
		pihatlogitmodel <- glm(r~w, family=binomial(link="logit"), data=mydata)	
		newdata = data.frame(w)
		pihat = data.frame(probs = predict(pihatlogitmodel, 
                                       newdata = newdata, 
                                       type="response"))
		pihat=pihat$probs
		return(pihat)

	
}


Compute_Estimators<-function(Dpop,n,N,p,M) 
{
	YMAR2M=rep(0,M)
	YMAR3M=rep(0,M)
	YMCAR2M=rep(0,M)
	YMCAR3M=rep(0,M)
	VMAR2M=rep(0,M)
	VMAR3M=rep(0,M)
	VMCAR2M=rep(0,M)
	VMCAR3M=rep(0,M)	
	
	for(i in 1:M)
	{
		Ds=DefineddataSample(Dpop,n,N,p) 
		DataEstimators=Compute_Estimators_eachrep(Ds,n,N,p,Dpop$x,Dpop$u, Dpop$k)
		YMAR2M[i]=DataEstimators[1]
		YMAR3M[i]=DataEstimators[2]
		YMCAR2M[i]=DataEstimators[3]
		YMCAR3M[i]=DataEstimators[4]
		VMAR2M[i]=DataEstimators[5]
		VMAR3M[i]=DataEstimators[6]
		VMCAR2M[i]=DataEstimators[7]
		VMCAR3M[i]=DataEstimators[8]
	
	}
	DEstimators=data.frame(YMAR2M,YMAR3M,YMCAR2M,YMCAR3M,VMAR2M,VMAR3M,VMCAR2M,VMCAR3M  )
	return(DEstimators)
}

Compute_Estimators_eachrep<-function(Datasample,n,N,p,xpop,upop, kpop) 
{
	y=Datasample$ysample
	x=Datasample$xsample
	u=Datasample$usample
	k=Datasample$ksample
	Phii=Datasample$Phiisample
	r=Datasample$rsample
	w1=Datasample$w1sample
	w2=Datasample$w2sample
	pihat=Createpihat(Datasample)
	K=sum(kpop)
	#--------Define for Construct GREG estimator------
		rn=diag(r)
		xn=matrix(x, ncol=1)
		Phiininves=diag(1/Phii)
		LN=rep(1,N)
		Ln=rep(1,n)
		LN=matrix(LN,ncol=1)
		Ln=matrix(Ln,ncol=1)
		UN<-c(LN,upop)
		Un<-c(Ln,u)
		MUN=matrix(UN, ncol=2)
		MUn=matrix(Un, ncol=2)
		MSumUN=t(MUN)%*%LN	
	#-------------------------------------------------
	
	
	
	#--------------------Compute value of new ratio estimator-------------------------
		#--------(1) MAR mechanism--------------------------
			#----------(1.1) Case known Xbar--------------------------------
				YbarRhat=1/N*(sum ( (r*y)/ (Phii*pihat)))
				XbarRhat=1/N*(sum ( (r*x)/ (Phii*pihat)))
				Xbar=mean(xpop)
				Rhat2= YbarRhat/XbarRhat
				YMAR2=  N*(YbarRhat/XbarRhat)*Xbar			# New ratio estimator
			#----------(1.2) Case unknown Xbar--------------------------------
				pininverse=diag(1/pihat)
				MSumUn=t(MUn)%*%pininverse%*%Phiininves%*%rn%*%Ln
				beta0=t(MUn)%*%rn%*%pininverse%*%Phiininves%*%MUn
				beta1=t(MUn)%*%rn%*%pininverse%*%Phiininves%*%x
				betaphat=solve(beta0)%*%beta1
				betahat1= betaphat
				ComputeXbarGREGhat=(1/N)*(sum((r*x)/(Phii*pihat))+t( MSumUN-MSumUn)%*%betaphat)
				XbarGREGhat=ComputeXbarGREGhat[1][1]
				YMAR3=N*(YbarRhat/XbarRhat)*XbarGREGhat	# New ratio estimator

		#--------(2) MCAR mechanism--------------------------
			#----------(2.1) Case known Xbar--------------------------------
				#phat=sum(r/Phii)/sum(1/Phii)
				phat=sum(r)/n
				

				YbarRhat=1/N*(sum ( (r*y)/ (Phii*phat)))
				XbarRhat=1/N*(sum ( (r*x)/ (Phii*phat)))
				Rhat2prime = YbarRhat/XbarRhat
				YMCAR2=N*(YbarRhat/XbarRhat)*Xbar			# New ratio estimator
			#----------(2.2) Case unknown Xbar--------------------------------				
			
				MSumUn=t(MUn)%*%Phiininves%*%rn%*%Ln/phat
				beta0=t(MUn)%*%rn%*%Phiininves%*%MUn/phat
				beta1=t(MUn)%*%rn%*%Phiininves%*%x/phat
				betaphat=solve(beta0)%*%beta1
				betahat2= betaphat
				ComputeXbarGREGhat=(1/N)*(sum((r*x)/(Phii*phat))+t( MSumUN-MSumUn)%*%betaphat)
				XbarGREGhat=ComputeXbarGREGhat[1][1]
				YMCAR3=N*(YbarRhat/XbarRhat)*XbarGREGhat	# New ratio estimator				
		Yhat<-c(YMAR2,YMAR3,YMCAR2,YMCAR3)
		
			
	#--------------------Compute variance estimators-------------------------			
		#--------(1) MAR mechanism--------------------------
			#----------(1.1) Case known Xbar--------------------------------
				ubetap=MUn%*%betahat1
				A2hat=y-Rhat2*x
				Bhat=x-Rhat2*ubetap
				D=(1-Phii)/Phii
				Dhat=r*D/(Phii*pihat)
				E=(1-pihat)/pihat
				Ehat=r*E/(Phii*pihat)
				VMART21= sum( Dhat*A2hat^2+Ehat*y^2)
				
				VMART31= sum( (Dhat+Ehat)*Bhat^2)
		#--------(2) MCAR mechanism--------------------------
			#----------(2.1) Case known Xbar--------------------------------
				ubetap=MUn%*%betahat2
				A2hatprime=y-Rhat2prime*x
				Bhatprime=x-Rhat2prime*ubetap
				Dhatprime=r*D/(Phii*phat)
				Eprime= (1-phat)/phat
				Eprimehat=r*E/(Phii*phat)
				VMCART21= sum( Dhatprime*A2hatprime^2+Eprimehat*y^2)
				VMCART31= sum( (Dhatprime+Eprimehat)*Bhatprime^2)
	#-------------------Compute Seconde Term of Variance-------------------
	VMART22=0
	VMCART22=0
	VMART32=0
	VMCART32=0	
	
	for(i in 1:n)
	{
		for(j in 1:n)
		{
			if(j!=i)
			{
				Phi_ij=((k[i]+k[j])/K)*((N-n)*(n-1))/((N-1)*(N-2))+((n-1)*(n-2))/((N-1)*(N-2))
				Dij=(Phi_ij-Phii[i]*Phii[j])/(Phii[i]*Phii[j])
				Dijhat=r[i]*r[j]*Dij/(Phi_ij*pihat[i]*pihat[j])
				
				Dijhatprime= r[i]*r[j]*Dij/(Phi_ij*(phat^2))
				VMART22=VMART22+Dijhat*A2hat[i]*A2hat[j]
				VMCART22=VMCART22+Dijhatprime*A2hatprime[i]*A2hatprime[j]				
				VMART32=VMART32+Dijhat*Bhat[i]*Bhat[j]
				VMCART32=VMCART32+Dijhatprime*Bhatprime[i]*Bhatprime[j]
			}
		}
	}
	VMAR2= VMART21 + VMART22
	VMCAR2= VMCART21 + VMCART22
	VMAR3= VMART31 + VMART32
	VMCAR3= VMCART31 + VMCART32
	Vhat <-c(VMAR2,VMAR3,VMCAR2,VMCAR3) 
	Estimator <-c(Yhat, Vhat)
}	

Create_Pi<-function(wpop) 
{
	w=wpop
	
	a=-0.1+ 0.04*w
	b=exp(-a)
	Pi= 1/(1+b)
	#print(Pi)
	return(Pi)
}


ComputeTrueVariance<-function(DataP,p,n,N)
{
	y=DataP$y		#Study variable
	x=DataP$x		#Auxiliary variable
	k=DataP$k		#Size variable	
	u=DataP$u		#Caribration variable
	w=DataP$w		# Auxiliary variable use estimated response probability
					# Auxiliary variable use estimated response probability
	r=DataP$r
	Truepi=DataP$Pi
	Pi=k/sum(k)
	K=sum(k)
	Phii=Pi*(N-n)/(N-1)+(n-1)/(N-1)
	D=(1-Phii)/Phii
#-------Ratio Estimator-----------------------	

	R=mean(y)/mean(x)
	A=y-R*x



	
#--------------Known X bar---------------------------------
	VMAR21= sum(D*A^2)
		
#-------------Unknown X bar----------------------------------
	LN=rep(1,N)
	UN<-c(LN,u)
	MUN=matrix(UN, ncol=2)
	MSumUN=t(MUN)%*%LN
	Beta=solve(t(MUN)%*%MUN)%*%(t(MUN)%*%x)
	#print(Beta)
	UBeta=MUN%*%Beta
	B=x-R*UBeta
	VMAR31= sum(D*B^2)
	
	
#---------------------Second Term----------------------------------
	VMAR22=0
	VMAR32=0
	for(i in 1:N)
	{
		for(j in 1:N)
		{
			if(j!=i)
			{
				Phi_ij=((k[i]+k[j])/K)*((N-n)*(n-1))/((N-1)*(N-2))+((n-1)*(n-2))/((N-1)*(N-2))
				Dij=(Phi_ij-Phii[i]*Phii[j])/(Phii[i]*Phii[j])
				VMAR22=	VMAR22+  Dij*A[i]*A[j]
				VMAR32=	VMAR32+  Dij*B[i]*B[j]

			}
		}
	}
	
#---------------------Third Term----------------------------------	
	E=(1-Truepi)/Truepi
	Eprime=(1-p)/p
	VMAR23=sum(E*y^2)
	VMAR33=sum(E*B^2)
	VMAR23prime=sum(Eprime*y^2)
	VMAR33prime=sum(Eprime*B^2)

#------------------------------------------------------------------
	VMAR2=VMAR21+VMAR22+VMAR23
	VMAR3=VMAR31+VMAR32+VMAR33
	VMCAR2=VMAR21+VMAR22+VMAR23prime
	VMCAR3=VMAR31+VMAR32+VMAR33prime	
	TrueVariance<-c(VMAR2,VMAR3,VMCAR2,VMCAR3)
	return(TrueVariance)
}

ComputeEff<-function(Truevalue,Estimator,M) 
{

	EsimYhat=sum(Estimator)/M
	RBias=EsimYhat-Truevalue
	RBsimYhat= abs(RBias)/Truevalue
	#MSE=1/(M-1)*sum((Yhat-Y)^2)
	RRMSE=sqrt(1/(M-1)*sum((Estimator-Truevalue)^2))/Truevalue
	SimRelative<-c(RBsimYhat,  RRMSE)
    return(SimRelative)
}	

CompareEfficiency<-function(TrueV,DEstimaors,SumY,M) 
{
	VMAR2=TrueV[1]
	VMAR3=TrueV[2]
	VMCAR2=TrueV[3]
	VMCAR3=TrueV[4]	
	
	YMAR2M=DEstimaors[1]
	YMAR3M=DEstimaors[2]	
	YMCAR2M=DEstimaors[3]
	YMCAR3M=DEstimaors[4]
	VMAR2M=DEstimaors[5]
	VMAR3M=DEstimaors[6]	
	VMCAR2M=DEstimaors[7]
	VMCAR3M=DEstimaors[8]
	EffYMAR2=ComputeEff(SumY,YMAR2M,M) 
	EffYMCAR2=ComputeEff(SumY,YMCAR2M,M) 
	EffYMAR3=ComputeEff(SumY,YMAR3M,M) 
	EffYMCAR3=ComputeEff(SumY,YMCAR3M,M) 

	
	EffVMAR2=ComputeEff(VMAR2,VMAR2M,M) 
	EffVMCAR2=ComputeEff(VMCAR2,VMCAR2M,M)
	EffVMAR3=ComputeEff(VMAR3,VMAR3M,M)
	EffVMCAR3=ComputeEff(VMCAR3,VMCAR3M,M)
	
	
	print(EffYMAR2[2])
	print(EffYMCAR2[2])
	print(EffYMAR3[2])
	print(EffYMCAR3[2])
	cat("---------------------------------------","\n")
	print(EffVMAR2[2])
	print(EffVMCAR2[2])
	print(EffVMAR3[2])
	print(EffVMCAR3[2]) 
	
}
GenPopulation<-function(N) 
{
	rho <- 0.70
	mu1 <-15; s1 <- 1
	mu2 <-5; s2 <- 1
# Parameters for bivariate normal distribution
	mu <- c(mu1,mu2) # Mean
	sigma <- matrix(c(s1^2, rho, rho, s2^2),2) # Covariance matrix
	xu <- mvrnorm(N, mu = mu, Sigma = sigma ) # from MASS package
	colnames(xu) <- c("x","u")
	xu=data.frame(xu)	   
	k<-rgamma(N, scale=10, shape=5)
	w<-rgamma(N, scale=5, shape=10)
	x=xu$x
	u=xu$u
	y <- abs(0.2*x+0.1*w+2*k + 3.7*sqrt(k) +2*u+ rnorm(N) )
	data=data.frame(y,x,w,k,u)
	return(data)
}	
#----------------------------------------------------------------------------

#--------------Main function-------------------------
library(MASS)
set.seed(123)
N=3000
n=1200
M=1000
Dpop=GenPopulation(N)
Pi=Create_Pi(Dpop$w)
#print(Pi)
r=rbern(N, prob = Pi)
#print(r)
print(n)

Dpop=data.frame(Dpop,r,Pi)
p=mean(Dpop$r)
print(p)
TrueV=ComputeTrueVariance(Dpop,p,n,N)
	
DEstimaors=Compute_Estimators(Dpop,n,N,p,M)

ComEff=CompareEfficiency(TrueV,DEstimaors,sum(Dpop$y),M)