N = 20
stratumi = list()
stratumi[[1]] = rbind(c(2, 2, 3, 4, 6, 7, 7, 3, 5, 4, 3, 8), 
                      c(17, 15, 14, 15, 18, 11, 10, 16, 18, 12, 11, 15))
stratumi[[2]] = rbind(c(6, 6, 6, 5, 6, 3, 7, 7),
                      c(10, 11, 16, 17, 18, 18, 17, 16))
(Ni = sapply(stratumi, ncol))

uzorak = list()
uzorak[[1]] = rbind(c(2, 6, 7, 4),
                    c(17, 18, 11, 12))
uzorak[[2]] = rbind(c(6, 6, 7),
                    c(11, 18, 17))
(ni = sapply(uzorak, ncol))

xni = sapply(uzorak, function(t) {mean(t[1, ])})
yni = sapply(uzorak, function(t) {mean(t[2, ])})

xstr = sum(Ni * xni)
ystr = sum(Ni * yni)

# kombinovana ocena
(Rc = xstr / ystr)
y = c(stratumi[[1]][2, ], stratumi[[2]][2, ])
(xrc = Rc * mean(y))

# posebna ocena
Ri = xni / yni
ty = sapply(stratumi, function(t) {sum(t[2, ])})
(xrs = sum(Ri * ty) / N)

# prava vrednost
x = c(stratumi[[1]][1, ], stratumi[[2]][1, ])
mean(x)

# disperzije
si2x = sapply(stratumi, function(t) {var(t[1, ])})
si2y = sapply(stratumi, function(t) {var(t[2, ])})
rho = sapply(stratumi, function(t) {cor(t[1, ], t[2, ])})

R = mean(x) / mean(y)
(Dxrc = sum(Ni^2 * (1 - ni/Ni) * (si2x + R^2*si2y - 2*R*rho*sqrt(si2x)*sqrt(si2y)) / ni) / (N^2))
(Dxrs = sum(Ni^2 * (1 - ni/Ni) * (si2x + Ri^2*si2y - 2*Ri*rho*sqrt(si2x)*sqrt(si2y)) / ni) / (N^2))