Programy R

### ESTYMATOR HT

HT_sr<-function(ys,N,pi_m) {

n<-length(ys)
est<-sum(ys/diag(pi_m))/N

eD2_1a<-0
for (i in 1:n) {
eD2_1a<-eD2_1a+((ys[i]/pi_m[i,i])^2)*(1-pi_m[i,i])
}
eD2_1b<-0
for (i in 1:n) {
for (k in 1:n) {
if (i!=k) {
eD2_1b<-eD2_1b+ ys[k]*ys[i]*(pi_m[k,i]-pi_m[k,k]*pi_m[i,i])/(pi_m[k,k]*pi_m[i,i]*pi_m[k,i])
}
}
}
eD2_1<-(eD2_1a+eD2_1b)/(N^2)

eD2_2<-0
for (i in 1:n) {
for (k in 1:n) {
if (k>i) {
eD2_2<-eD2_2+ (pi_m[k,k]*pi_m[i,i]-pi_m[k,i]) *( ( (ys[k]/pi_m[k,k])-(ys[i]/pi_m[i,i]) )^2)/pi_m[k,i]
}
}
}
eD2_2<-eD2_2/(N^2)

print(est)
print(eD2_1)
print(eD2_2)
}


### PRZYKŁAD

pi_m<-matrix(c(0.4,0.01,0.015,0.01,0.3,0.02,0.015,0.02,0.3),ncol=3,nrow=3)

n<-10
N<-100

pi_m2<-matrix( (n*(n-1)/(N*(N-1))),ncol=10,nrow=10)
diag(pi_m2)<-n/N

rm(n)
rm(N)


HT_sr(c(1:3),100,pi_m)
HT_sr(c(1:10),100,pi_m2)


###> HT_sr(c(1:3),100,pi_m)
###[1] 0.1916667
###[1] -0.1078472
###[1] 0.06236111
###> HT_sr(c(1:10),100,pi_m2)
###[1] 5.5
###[1] 0.825
###[1] 0.825