N = 4
n = 2

br_tigrova = c(11, 20, 23, 54)
Mi = c(100, 200, 300, 500)
M = sum(Mi)

# u uzorak odabrane 1. i 2. jedinka
indeks = c(1,2)
uzorak = br_tigrova[indeks]
uzorak

pi = Mi/M
pi
sum(pi)

# ako sa ponavljanjem
t_hh = sum(uzorak/pi[indeks])/n
t_hh
t = sum(br_tigrova)
t

D_t_hh = sum(pi*(br_tigrova/pi-t)^2)/n
D_t_hh

D_t_hh_ocena = sum((uzorak/pi[indeks]-t_hh)^2)/(n*(n-1))
D_t_hh_ocena

# ako bez ponavljanja
# za HT ocenu potrebne su nam v-ce ukljucenja
# uzorak bp -> ne vazi pii = (1-pi)^n u opstem slucaju!
# pii = suma v-ca svih uzoraka koji sadrze i-tu jedinku
uzorci = list()
obelezje_na_uzorku = list()
i = 1
for(j in 1:N) {
  for(k in 1:N) {
    if(j<k) {
      uzorci[[i]] = c(j,k)
      obelezje_na_uzorku[[i]] = c(br_tigrova[j],br_tigrova[k])
      i = i+1
    }
  }
}
uzorci
obelezje_na_uzorku

# v-ca uzorka (i,j) -> 
# P(odaberemo i)P(odaberemo j|odabrali smo i) + P(odaberemo j)P(odaberemo i|odabrali smo j)
# = (Mi/M)*(Mj/(M-Mi)) + (Mj/M)*(Mi/(M-Mj)) 
vca_uzorka = function(i,j) {
  (Mi[i]/M)*(Mi[j]/(M-Mi[i])) + (Mi[j]/M)*(Mi[i]/(M-Mi[j]))
}

p_uzorka = c()
for(i in 1:length(uzorci)) {
  p_uzorka[i] = vca_uzorka(uzorci[[i]][1], uzorci[[i]][2])
}
p_uzorka
sum(p_uzorka)

pii = rep(0,N)
for(i in 1:N) {
  for(j in 1:length(uzorci)) {
    if(i %in% uzorci[[j]]) {
      pii[i] = pii[i]+p_uzorka[j]
    }
  }
}
pii
# provera - sum(pii)=n, za uzorak bez ponavljanja
# (nekad moze da vrati FALSE iako zapravo jeste =, zbog zaokruzivanja)
sum(pii)==n

t_ht = sum(uzorak/pii[indeks])
t_ht
sum(br_tigrova)

D_t_ht = sum((1-pii)*br_tigrova^2/pii)
for(i in 1:N) {
  for(j in 1:N) {
    if(i!=j) {
      pij = vca_uzorka(i,j)
      D_t_ht = D_t_ht + (pij-pii[i]*pii[j])*(br_tigrova[i]*br_tigrova[j])/(pii[i]*pii[j])
    }
  }
}
D_t_ht

D_t_ht_ocena = sum((1-pii[indeks])*uzorak^2/(pii[indeks]^2))
for(i in indeks) {
  for(j in indeks) {
    if(i!=j) {
      pij = vca_uzorka(i,j)
      D_t_ht_ocena = D_t_ht_ocena + (pij-pii[i]*pii[j])*(uzorak[i]*uzorak[j])/(pii[i]*pii[j]*pij)
    }
  }
}
D_t_ht_ocena

D_t_ht_ocena_syg = 0
for (i in indeks) {
  for(j in indeks) {
    if(i<j) {
      D_t_ht_ocena_syg = D_t_ht_ocena_syg + (pii[i]*pii[j]-pij)*((uzorak[i]/pii[i]-uzorak[j]/pii[j])^2)/pij
    }
  }
}
D_t_ht_ocena_syg