Log of Mu 281 Ratio Estimation Example
=================================== 9/21/04

> cor(mu284[,c(2,3,9)])
            P85       P75     REV84
P85   1.0000000 0.9983968 0.9781463
P75   0.9983968 1.0000000 0.9745424
REV84 0.9781463 0.9745424 1.0000000

## Now compare empirical MSE of 1000 ordinary 
## estimators of mu281 Y = P85 total with ratio 
## estimator based on X=P75 total or X=REV84 total
## from SRS's with n=50 from mu281

## Recall that indSm are the indices in 1:284 of
## the 281 smaller cities

> thatmat <- array(0, dim=c(1000,3), dimnames=
    list(NULL, c("t_P85","tr_P75","tr_REV84")))
> for (i in 1:1000) {
    isam <- sample(indSm, 50)
    ysbar <- mean(mu284$P85[isam])
    thatmat[i,] <- ysbar/c(1, mean(mu284$P75[isam]),
       mean(mu284$REV84[isam])) }
  thatmat <- thatmat %*% diag(c(281, sum(mu284$P75[indSm]),
       sum(mu284$REV84[indSm])))
  round(apply( thatmat - sum(mu284$P85[indSm]), 2, 
       function(col) mean(col^2) ))
[1] 750782   7461  129105

## Empirical ratios of two latter MSE's to first are
[1] 0.00994 0.17196

## Theoretical values:
> round(mean( (data.matrix(mu284[indSm,c(2,3)]) %*% c(1, 
    -mean(mu284$P85[indSm])/mean(mu284$P75[indSm])) )^2 )/
    (280*var(mu284$P85[indSm])/281),5)
[1] 0.0099
> round(mean( (data.matrix(mu284[indSm,c(2,9)]) %*% c(1, 
    -mean(mu284$P85[indSm])/mean(mu284$REV84[indSm])) )^2 )/
    (280*var(mu284$P85[indSm])/281),5)
[1] 0.16605

## So the theoretical ratios of MSE values are  
##   very close to the empirical ones !!!
## and both sets are MUCH less than 1, 
##   so ratio estimation works !!