uneq

##function that computes HH estimator and its variance estimator

VHH=function(y,x,n){
  p=x/sum(x)
  N=length(y)
  sa=sample(1:N,n,replace = TRUE, prob =p )
  tauh=(1/n)*sum(y[sa]/p[sa])
  vhat=(1/(n*(n-1)))*sum((y[sa]/p[sa]-tauh)^2)
  res=c(est=tauh,varh=vhat)
  res}

# simulation on trees data

VHH(trees$Height,trees$Girth,10)

##       est      varh 
##  2295.966 21428.932

sum(trees$Height)

## [1] 2356

p=trees$Girth/sum(trees$Girth)
(1/10)*sum(p*((trees$Height/p)-sum(trees$Height))^2)

## [1] 22393.9

v=numeric(0)
for(i in 1:10000)
  v=rbind(v,VHH(trees$Height,trees$Girth,10))
res=apply(v,2,mean)

#estimator of Total and its variance estimator

res

##       est      varh 
##  2357.587 22422.108

#Histogram of estimator and its variance estimator

mean1=v[,1]
varhat=v[,2]
hist(mean1)

hist(varhat)