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




Createpihat<-function(r,w)
{
	pihat0=rep(0,n)
	################# w1 and w2 variables ##########################
		mydata=data.frame(r,w) 
		pihatlogitmodel <- glm(r~w, family=binomial(link="logit"), data=mydata)	
	 #print(pihatlogitmodel )
		newdata = data.frame(w)
		#pihat=predict(pihatlogitmodel, newdata, type="response") #Estimated 
		pihat = data.frame(probs = predict(pihatlogitmodel, 
                                       newdata = newdata, 
                                       type="response"))
		pihat=pihat$probs
		return(pihat)

	
}


Compute_Estimators_eachrep<-function(Datasample,n,N,Xbar, Sumk,Sumu) 
{
	y=Datasample$y
	x=Datasample$x
	u=Datasample$u
	k=Datasample$k
	Pi=k/Sumk
	Phii=Pi*(N-n)/(N-1)+(n-1)/(N-1)
	r=Datasample$r
	w=Datasample$w
	pihat=Createpihat(r,w)
	K=Sumk
	#--------Define for Construct GREG estimator------
		rn=diag(r)
		xn=matrix(x, ncol=1)
		Phiininves=diag(1/Phii)
		#LN=rep(1,N)
		LN=1
		Ln=rep(1,n)
		#LN=matrix(LN,ncol=1)
		Ln=matrix(Ln,ncol=1)
		UN<-c(N,Sumu)
		Un<-c(Ln,u)
		MUN=matrix(UN, ncol=2)
		
		MUn=matrix(Un, ncol=2)
		MSumUN=t(MUN)

	#-------------------------------------------------
	
	
	
	#--------------------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)))
				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%*%xn
				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%*%xn/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--------------------------
			
				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--------------------------
			
				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)
}	

ShowOutput<-function(Estimaors) 
{
	cat("\n\n\t\t\t== Application to real data ==\n")
	cat("----------------------The propose estimators-------------------------- \n")
	cat("\t - MAR mechanism wiht Xbar is known =" ,Estimaors[1] , "\n" )
	cat("\t - MAR mechanism wiht Xbar is unknown =" ,Estimaors[2] , "\n" )
	cat("\t - MCAR mechanism wiht Xbar is known =" ,Estimaors[3] , "\n" )
	cat("\t - MCAR mechanism wiht Xbar is unknown =" ,Estimaors[4] , "\n\n" )
	cat("\n ---------------The variance estimators of the propose estimators--------------------- \n")
	cat("\t - MAR mechanism wiht Xbar is known =" ,Estimaors[5] , "\n" )
	cat("\t - MAR mechanism wiht Xbar is unknown =" ,Estimaors[6] , "\n" )
	cat("\t - MCAR mechanism wiht Xbar is known =" ,Estimaors[7] , "\n" )
	cat("\t - MCAR mechanism wiht Xbar is unknown =" ,Estimaors[8] , "\n\n" )	
}	
#--------------Main function-------------------------
	
#---------------Data Water--------------------------------


r<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,0)

y<-c(1253533,1124353,715855,423005,151111,1389494,10483537,1613897,969195,43266,2890963,719526,1183116,889614,384293,1492636,1198977,0,273722,487488,274694,878693,431775,3815825,1332621,939597,3979414,509964,856492,5166899,778063,1875844,218754,721606,1288233,2812220,722004,0,327288,0)

x<-c(1844988,1548186,1519620,639943,182480,892125,17115581,2411020,868976,76703,4023855,921145,2531038,1192909,598131,2176203,1568352,0,347835,635039,389618,1194184,629823,5351285,2371879,1192060,5283104,953367,1670178,7363236,950121,2836869,302118,931026,1879134,4486879,757496,0,543722,0)

k<-c(61120,67738,35657,22946,7447,74873,379521,77573,49308,3037,99135,39983,57858,54457,24586,67338,49835,19961,14265,25373,14512,46666,21048,185859,60020,45694,183642,28471,37161,94204,35475,74402,13685,48348,70786,143478,39331,14240,18203,12173)

u<-c(1850026,1544901,1520054,661296,179513,1866284,17412038,2401814,1468833,77592,4096479,941949,2568120,1200445,600461,2157585,1652496,692822,338520,642301,405153,1238932,690925,5466569,2272187,1202315,5243405,967033,1694191,7452961,974481,2850660,302886,953179,1888601,4471663,934060,345706,551663,307016)
w<-c(1269570,1123472,713610,428064,147397,1433990,10270556,1610301,1019684,43027,2847088,741428,1196080,891924,384115,1468991,1258054,435326,263630,489046,259211,906650,431552,4051922,1315910,956980,3860605,502254,816439,5147622,791660,1887604,219792,733502,1304705,2757629,724985,273170,326976,224608)

#--------------------------Main function--------------------------------
DataSample=data.frame(r,y,x,k,u,w)

N=74
n=40
Sumx=180040825
Xbar=Sumx/N
Sumk=5112745
Sumu=169325493
Estimaors=Compute_Estimators_eachrep(DataSample,n,N,Xbar, Sumk,Sumu)
ShowOutput(Estimaors) 

#---------------------------------------------------
#source("D:/04_CodeProgram/CodePaper65/CodeRealDataWaterDeand.R")