Comparing the HT estimator and the ratio estimator

Here is some R code that should help you do HW assignment 6. You should read the details of assignment 6 in the syllabus carefully before looking at this code.

The only code you will have to change is the last 12 lines in the window which are those after the line set.seed(2004). Of course you will have to consider various designs. It will probably be helpful to look at other functions of popx as possible designs.

getans<-function(est,stnderr,truetot) { err<-abs(est-truetot) lwbd<-est -1.96*stnderr upbd<-est + 1.96*stnderr if(lwbd <= truetot & truetot <= upbd) { cov<-1} else {cov<-0} ans<-c(est,err,lwbd,upbd - lwbd,cov) return(ans) } srstot<-function(smp,popy) { n<-length(smp) N<-length(popy) fpc<-1-n/N ysmp<-popy[smp] est<-N*mean(ysmp) stnderr<-N*sqrt((fpc/n)*var(ysmp)) ans<-getans(est,stnderr,sum(popy)) return(ans) } ratiotot<-function(smp,popy,popx) { n <- length(smp) N<-length(popx) ff<-n/N ysamp<-popy[smp] xsamp<-popx[smp] tx<-sum(popx) trtot<-sum(popy) rhat <- sum(ysamp)/sum(xsamp) est <- rhat * tx dum1<-(N*N*(1-ff))/(n*(n-1)) vartot <- dum1*sum((ysamp-rhat*xsamp)^2) stnderr<-sqrt(vartot) ans<-getans(est,stnderr,sum(popy)) return(ans) } httot<-function(smp,popy,designwts) { wts<-designwts[smp] est<-sum(wts*popy[smp]) n<-length(smp) dum<-sum((n*wts*popy[smp] - est)^2) stnderr<-sqrt((1/(n*(n-1))*dum)) ans<-getans(est,stnderr,sum(popy)) return(ans) } compar3est<-function(popy,popx,design,n) { designwts<-sum(design)/(n*design) smp<-sample(1:length(popy),n,replace=FALSE,prob=design) anssrs<-srstot(smp,popy) ansHT<-httot(smp,popy,designwts) ansratio<-ratiotot(smp,popy,popx) ans<-rbind(anssrs,ansHT,ansratio) return(ans) } compar3estlp<-function(popy,popx,design,n,R) { ans<-matrix(0,3,5) for(i in 1:R){ ans<-ans + compar3est(popy,popx,design,n) } ans<-round(ans/R,digits=3) return(ans) } set.seed(2004) popx<-rlnorm(500,10,.4) hist(popx) popy<-rnorm(500,0.1*popx,3*sqrt(popx)) plot(popx,popy) sum(popy) cor(popx,popy) n<-50 R<-500 design<-popx compar3estlp(popy,popx,design,n,R) design<-rep(1,length(popy)) compar3estlp(popy,popx,design,n,R)

External Data Entry

Enter a Dataset URL :