Dzeknajf
Implementiramo funkciju koja racuna Dzeknajf pseudovrednosti. Prosledjujemo joj vektor podataka i funkciju koja racuna ocenu na osnovu podataka.
pseudo_jack <- function(x, h) {
if (is.vector(x)) {
n <- length(x)
} else {
n <- nrow(x)
}
leave_one_out <- sapply(1:n, function(i) h(subset(x, 1:n != i)))
pseudo_values <- n*h(x) - (n-1)*leave_one_out
pseudo_values
}Kada imamo pseudo vrednosti, lako racunamo sve sto nam je potrebno za Dzeknajf.
Dzeknajf ocenu funkcionala h, ocenu pristrasnosti ocene \(\hat h\) i ocenu disperzije ocene \(\hat h\) racunamo sa:
h_jack <- function(x, h) mean(pseudo_jack(x, h))
b_jack <- function(x, h) h(x) - h_jack(x, h)
V_jack <- function(x, h) {
pj <- pseudo_jack(x, h)
var(pj) / length(pj)
}Generisimo podatke iz rada Efron-a i Stein-a o LSAT rezultatima u americi.
X_efron <- data.frame(LSAT = c(576, 635, 558, 578, 666, 580, 555, 661, 651, 605, 653, 575, 545, 572, 594),
GPA = c(3.39, 3.30, 2.81, 3.03, 3.44, 3.07, 3.00, 3.43, 3.36, 3.13, 3.12, 2.74, 2.76, 2.88, 2.96))\(S\) ce biti uzoracka korelacija izmedju LSAT i GPA i izracunacemo ocene kvaliteta te ocene.
S <- function(X) cor(X$LSAT, X$GPA)
cat("S = ", S(X_efron), '\n',
"h_jack = ", h_jack(X_efron, S), '\n',
"b_jack = ", b_jack(X_efron, S), '\n',
"V_jack = ", V_jack(X_efron, S),
sep = "")S = 0.7763745
h_jack = 0.7828481
b_jack = -0.006473623
V_jack = 0.02031156Efron i Stein u svom radu dokazuju da je ocena V_jack pozitivno pristrasna i daju mogucu korekciju te pristrasnosti, koju cemo implementirati.
V_jack_corr <- function(x, h, vjack = NULL) {
if (is.vector(x)) {
n <- length(x)
} else {
n <- nrow(x)
}
if (is.null(vjack))
vjack <- V_jack(x, h)
Qii <- apply(combn(1:n, 2), 2, function(ij) {
n*h(x) - (n-1)*(h(subset(x, 1:n != ij[1])) + h(subset(x, 1:n != ij[2]))) +
(n-2)*(h(subset(x, !(1:n %in% ij))))
})
vjack - ((n*(n-1)/2 - 1) / (n*(n+1)))*var(Qii)
}Korigovana ocena disperzije ocene S je
cat("V_jack_corr = ", V_jack_corr(X_efron, S), sep = "")V_jack_corr = 0.01818288Testirajmo rezultate za neki vestacki uzorak za koji znamo raspodelu
set.seed(825)
n <- 50
X <- rnorm(n)
cat("Var = ", var(X), '\n',
"h_jack = ", h_jack(X, var), '\n',
"b_jack = ", b_jack(X, var), '\n',
"var(var) = ", 2/(n-1), '\n',
"V_jack = ", V_jack(X, var), '\n',
"V_jack_corr = ", V_jack_corr(X, var),
sep = "")Var = 1.029494
h_jack = 1.029494
b_jack = 1.332268e-15
var(var) = 0.04081633
V_jack = 0.05423948
V_jack_corr = 0.0533997Formiramo listu ocena za koje cemo porediti disperzije i dzeknajf ocene istih
estimates <- list(mean = mean, var = var, inv_mean = function(x) 1/mean(x), cv = function(x) sd(x)/mean(x))Izracunamo ocene, prave disperzije i dzeknajf ocene tih disperzija za normalnu, eksponencijalnu i geometrijsku raspodelu.
n <- 50
N <- 1000
set.seed(406)
normal <- rnorm(n, 1)
norm_results <- sapply(estimates, function(f)
c(est = f(normal),
var = var(replicate(N, f(rnorm(n, 1)))),
V_jack = V_jack(normal, f))
)
set.seed(406)
exponential <- rexp(n)
exp_results <- sapply(estimates, function(f)
c(est = f(exponential),
var = var(replicate(N, f(rexp(n)))),
V_jack = V_jack(exponential, f))
)
set.seed(406)
geometric <- rgeom(n, 0.5)
geom_results <- sapply(estimates, function(f)
c(est = f(geometric),
var = var(replicate(N, f(rgeom(n, 0.5)))),
V_jack = V_jack(geometric, f))
)i prikazemo rezultate:
results <- round(t(cbind(norm_results, exp_results, geom_results)), 3)
rownames(results) <- c("N(1,1) mean",
"N(1,1) var",
"N(1,1) 1/mean",
"N(1,1) sd/mean",
"Exp(1) mean",
"Exp(1) var",
"Exp(1) 1/mean",
"Exp(1) sd/mean",
"Geom(0.5) mean",
"Geom(0.5) var",
"Geom(0.5) 1/mean",
"Geom(0.5) sd/mean")
colnames(results) <- c("estimate", "true variance", "jackknife variance")
results estimate true variance jackknife variance
N(1,1) mean 0.927 0.018 0.030
N(1,1) var 1.504 0.038 0.059
N(1,1) 1/mean 1.078 0.022 0.041
N(1,1) sd/mean 1.322 0.036 0.064
Exp(1) mean 1.117 0.020 0.025
Exp(1) var 1.271 0.160 0.089
Exp(1) 1/mean 0.895 0.020 0.017
Exp(1) sd/mean 1.009 0.018 0.009
Geom(0.5) mean 0.900 0.039 0.040
Geom(0.5) var 2.010 0.699 0.419
Geom(0.5) 1/mean 1.111 0.049 0.069
Geom(0.5) sd/mean 1.575 0.041 0.046Interesantan rezultat!
Generisacemo uzorak od 10 hiljada dzeknajf ocena disperzije za popravljenu uzoracku disperziju nad uzorkom obima 50 iz standardne normalne raspodele i analizirati raspodelu dzeknajf ocene.
set.seed(13)
# generisemo uzorak
jack_smp <- replicate(10000, V_jack(rnorm(n, 1), var))Teorijska disperzija ocene je 2/49 = 0.0408 Mi dobijamo
mean(jack_smp)[1] 0.04178114sto je vece od teorijske, sto i teorijski rezultati o dzeknajf oceni pokazuju. Ali pogledajmo ocenu gustine ovog uzorka:
plot(density(jack_smp))
abline(v = 2/49, col = "green")abline(v = mean(jack_smp), col = "red")
Vidimo da je ova raspodela asimetricna, sa zelenom vertikalnom linijom smo oznacili teorijsku disperziju, a sa crvenom ocekivanje \(\hat V_{jack}\) ocene. Uporedimo verovatnoce da dobijemo vrednosti manju od teorijske disperzije sa onim da dobijemo vece vrednosti od teorijske disperzije.
mean(jack_smp < 2/49)[1] 0.5809mean(jack_smp >= 2/49)[1] 0.4191Veca je verovatnoca da dobijemo ocenu disperzije koja je manja od stvarne disperzije, iako je ocekivanje ocene disperzije vece od stvarne.
Stoga se ne mozemo previse osloniti da nam Dzeknajf daje ocene koje premasuju teorijsku vrednost.
LS0tCnRpdGxlOiAiRHpla25hamYiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCkltcGxlbWVudGlyYW1vIGZ1bmtjaWp1IGtvamEgcmFjdW5hIER6ZWtuYWpmIHBzZXVkb3ZyZWRub3N0aS4KUHJvc2xlZGp1amVtbyBqb2ogdmVrdG9yIHBvZGF0YWthIGkgZnVua2NpanUga29qYSByYWN1bmEgb2NlbnUgbmEgb3Nub3Z1IHBvZGF0YWthLgpgYGB7cn0KcHNldWRvX2phY2sgPC0gZnVuY3Rpb24oeCwgaCkgewogIGlmIChpcy52ZWN0b3IoeCkpIHsKICAgIG4gPC0gbGVuZ3RoKHgpCiAgfSBlbHNlIHsKICAgIG4gPC0gbnJvdyh4KQogIH0KICBsZWF2ZV9vbmVfb3V0IDwtIHNhcHBseSgxOm4sIGZ1bmN0aW9uKGkpIGgoc3Vic2V0KHgsIDE6biAhPSBpKSkpCiAgcHNldWRvX3ZhbHVlcyA8LSBuKmgoeCkgLSAobi0xKSpsZWF2ZV9vbmVfb3V0CiAgcHNldWRvX3ZhbHVlcwp9CmBgYAoKS2FkYSBpbWFtbyBwc2V1ZG8gdnJlZG5vc3RpLCBsYWtvIHJhY3VuYW1vIHN2ZSBzdG8gbmFtIGplIHBvdHJlYm5vIHphIER6ZWtuYWpmLgoKRHpla25hamYgb2NlbnUgZnVua2Npb25hbGEgaCwgb2NlbnUgcHJpc3RyYXNub3N0aSBvY2VuZSAkXGhhdCBoJCBpIG9jZW51IGRpc3BlcnppamUgb2NlbmUgJFxoYXQgaCQgcmFjdW5hbW8gc2E6CmBgYHtyfQpoX2phY2sgPC0gZnVuY3Rpb24oeCwgaCkgbWVhbihwc2V1ZG9famFjayh4LCBoKSkKYl9qYWNrIDwtIGZ1bmN0aW9uKHgsIGgpIGgoeCkgLSBoX2phY2soeCwgaCkKVl9qYWNrIDwtIGZ1bmN0aW9uKHgsIGgpIHsKICBwaiA8LSBwc2V1ZG9famFjayh4LCBoKQogIHZhcihwaikgLyBsZW5ndGgocGopCn0KYGBgCgpHZW5lcmlzaW1vIHBvZGF0a2UgaXogcmFkYSBFZnJvbi1hIGkgU3RlaW4tYSBvIExTQVQgcmV6dWx0YXRpbWEgdSBhbWVyaWNpLgpgYGB7cn0KWF9lZnJvbiA8LSBkYXRhLmZyYW1lKExTQVQgPSBjKDU3NiwgNjM1LCA1NTgsIDU3OCwgNjY2LCA1ODAsIDU1NSwgNjYxLCA2NTEsIDYwNSwgNjUzLCA1NzUsIDU0NSwgNTcyLCA1OTQpLAogICAgICAgICAgICAgICAgR1BBID0gYygzLjM5LCAzLjMwLCAyLjgxLCAzLjAzLCAzLjQ0LCAzLjA3LCAzLjAwLCAzLjQzLCAzLjM2LCAzLjEzLCAzLjEyLCAyLjc0LCAyLjc2LCAyLjg4LCAyLjk2KSkKYGBgCgokUyQgY2UgYml0aSB1em9yYWNrYSBrb3JlbGFjaWphIGl6bWVkanUgTFNBVCBpIEdQQSBpIGl6cmFjdW5hY2VtbyBvY2VuZSBrdmFsaXRldGEgdGUgb2NlbmUuCmBgYHtyfQpTIDwtIGZ1bmN0aW9uKFgpIGNvcihYJExTQVQsIFgkR1BBKQpjYXQoIlMgPSAiLCBTKFhfZWZyb24pLCAnXG4nLAogICAgImhfamFjayA9ICIsIGhfamFjayhYX2Vmcm9uLCBTKSwgJ1xuJywKICAgICJiX2phY2sgPSAiLCBiX2phY2soWF9lZnJvbiwgUyksICdcbicsCiAgICAiVl9qYWNrID0gIiwgVl9qYWNrKFhfZWZyb24sIFMpLApzZXAgPSAiIikKYGBgCgpFZnJvbiBpIFN0ZWluIHUgc3ZvbSByYWR1IGRva2F6dWp1IGRhIGplIG9jZW5hIFZfamFjayBwb3ppdGl2bm8gcHJpc3RyYXNuYSBpIGRhanUgbW9ndWN1IGtvcmVrY2lqdSB0ZSBwcmlzdHJhc25vc3RpLCBrb2p1IGNlbW8gaW1wbGVtZW50aXJhdGkuCmBgYHtyfQpWX2phY2tfY29yciA8LSBmdW5jdGlvbih4LCBoLCB2amFjayA9IE5VTEwpIHsKICBpZiAoaXMudmVjdG9yKHgpKSB7CiAgICBuIDwtIGxlbmd0aCh4KQogIH0gZWxzZSB7CiAgICBuIDwtIG5yb3coeCkKICB9CiAgCiAgaWYgKGlzLm51bGwodmphY2spKQogICAgdmphY2sgPC0gVl9qYWNrKHgsIGgpCiAgCiAgUWlpIDwtIGFwcGx5KGNvbWJuKDE6biwgMiksIDIsIGZ1bmN0aW9uKGlqKSB7CiAgICAgICAgICAgICBuKmgoeCkgLSAobi0xKSooaChzdWJzZXQoeCwgMTpuICE9IGlqWzFdKSkgKyBoKHN1YnNldCh4LCAxOm4gIT0gaWpbMl0pKSkgKwogICAgICAgICAgICAgKG4tMikqKGgoc3Vic2V0KHgsICEoMTpuICVpbiUgaWopKSkpCiAgICAgICAgIH0pCiAgCiAgdmphY2sgLSAoKG4qKG4tMSkvMiAtIDEpIC8gKG4qKG4rMSkpKSp2YXIoUWlpKQp9CmBgYAoKS29yaWdvdmFuYSBvY2VuYSBkaXNwZXJ6aWplIG9jZW5lIFMgamUKCmBgYHtyfQpjYXQoIlZfamFja19jb3JyID0gIiwgVl9qYWNrX2NvcnIoWF9lZnJvbiwgUyksIHNlcCA9ICIiKQpgYGAKClRlc3RpcmFqbW8gcmV6dWx0YXRlIHphIG5la2kgdmVzdGFja2kgdXpvcmFrIHphIGtvamkgem5hbW8gcmFzcG9kZWx1CmBgYHtyfQpzZXQuc2VlZCg4MjUpCm4gPC0gNTAKWCA8LSBybm9ybShuKQpjYXQoIlZhciA9ICIsIHZhcihYKSwgJ1xuJywKICAgICJoX2phY2sgPSAiLCBoX2phY2soWCwgdmFyKSwgJ1xuJywKICAgICJiX2phY2sgPSAiLCBiX2phY2soWCwgdmFyKSwgJ1xuJywKICAgICJ2YXIodmFyKSA9ICIsIDIvKG4tMSksICdcbicsCiAgICAiVl9qYWNrID0gIiwgVl9qYWNrKFgsIHZhciksICdcbicsCiAgICAiVl9qYWNrX2NvcnIgPSAiLCAgVl9qYWNrX2NvcnIoWCwgdmFyKSwKc2VwID0gIiIpCmBgYAoKRm9ybWlyYW1vIGxpc3R1IG9jZW5hIHphIGtvamUgY2VtbyBwb3JlZGl0aSBkaXNwZXJ6aWplIGkgZHpla25hamYgb2NlbmUgaXN0aWgKYGBge3J9CmVzdGltYXRlcyA8LSBsaXN0KG1lYW4gPSBtZWFuLCB2YXIgPSB2YXIsIGludl9tZWFuID0gZnVuY3Rpb24oeCkgMS9tZWFuKHgpLCBjdiA9IGZ1bmN0aW9uKHgpIHNkKHgpL21lYW4oeCkpCmBgYAoKSXpyYWN1bmFtbyBvY2VuZSwgcHJhdmUgZGlzcGVyemlqZSBpIGR6ZWtuYWpmIG9jZW5lIHRpaCBkaXNwZXJ6aWphIHphIG5vcm1hbG51LCBla3Nwb25lbmNpamFsbnUgaSBnZW9tZXRyaWpza3UgcmFzcG9kZWx1LgpgYGB7cn0Kc2V0LnNlZWQoNDA2KQpuIDwtIDUwCk4gPC0gMTAwMApub3JtYWwgPC0gcm5vcm0obiwgMSkKbm9ybV9yZXN1bHRzIDwtIHNhcHBseShlc3RpbWF0ZXMsIGZ1bmN0aW9uKGYpIAogIGMoZXN0ID0gZihub3JtYWwpLAogICAgdmFyID0gdmFyKHJlcGxpY2F0ZShOLCBmKHJub3JtKG4sIDEpKSkpLAogICAgVl9qYWNrID0gVl9qYWNrKG5vcm1hbCwgZikpCikKCmV4cG9uZW50aWFsIDwtIHJleHAobikKCmV4cF9yZXN1bHRzIDwtIHNhcHBseShlc3RpbWF0ZXMsIGZ1bmN0aW9uKGYpIAogIGMoZXN0ID0gZihleHBvbmVudGlhbCksCiAgICB2YXIgPSB2YXIocmVwbGljYXRlKE4sIGYocmV4cChuKSkpKSwKICAgIFZfamFjayA9IFZfamFjayhleHBvbmVudGlhbCwgZikpCikKCmdlb21ldHJpYyA8LSByZ2VvbShuLCAwLjUpCgpnZW9tX3Jlc3VsdHMgPC0gc2FwcGx5KGVzdGltYXRlcywgZnVuY3Rpb24oZikgCiAgYyhlc3QgPSBmKGdlb21ldHJpYyksCiAgICB2YXIgPSB2YXIocmVwbGljYXRlKE4sIGYocmdlb20obiwgMC41KSkpKSwKICAgIFZfamFjayA9IFZfamFjayhnZW9tZXRyaWMsIGYpKQopCmBgYAppIHByaWthemVtbyByZXp1bHRhdGU6CmBgYHtyfQpyZXN1bHRzIDwtIHJvdW5kKHQoY2JpbmQobm9ybV9yZXN1bHRzLCBleHBfcmVzdWx0cywgZ2VvbV9yZXN1bHRzKSksIDMpCnJvd25hbWVzKHJlc3VsdHMpIDwtIGMoIk4oMSwxKSBtZWFuIiwKICAgICAgICAgICAgICAgICAgICAgICAiTigxLDEpIHZhciIsCiAgICAgICAgICAgICAgICAgICAgICAgIk4oMSwxKSAxL21lYW4iLAogICAgICAgICAgICAgICAgICAgICAgICJOKDEsMSkgc2QvbWVhbiIsCiAgICAgICAgICAgICAgICAgICAgICAgIkV4cCgxKSBtZWFuIiwKICAgICAgICAgICAgICAgICAgICAgICAiRXhwKDEpIHZhciIsCiAgICAgICAgICAgICAgICAgICAgICAgIkV4cCgxKSAxL21lYW4iLAogICAgICAgICAgICAgICAgICAgICAgICJFeHAoMSkgc2QvbWVhbiIsCiAgICAgICAgICAgICAgICAgICAgICAgIkdlb20oMC41KSBtZWFuIiwKICAgICAgICAgICAgICAgICAgICAgICAiR2VvbSgwLjUpIHZhciIsCiAgICAgICAgICAgICAgICAgICAgICAgIkdlb20oMC41KSAxL21lYW4iLAogICAgICAgICAgICAgICAgICAgICAgICJHZW9tKDAuNSkgc2QvbWVhbiIpCgpjb2xuYW1lcyhyZXN1bHRzKSA8LSBjKCJlc3RpbWF0ZSIsICJ0cnVlIHZhcmlhbmNlIiwgImphY2trbmlmZSB2YXJpYW5jZSIpCgpyZXN1bHRzCmBgYAoKLS0tLS0tLS0tLQoKIyMjIEludGVyZXNhbnRhbiByZXp1bHRhdCEKCkdlbmVyaXNhY2VtbyB1em9yYWsgb2QgMTAgaGlsamFkYSBkemVrbmFqZiBvY2VuYSBkaXNwZXJ6aWplIHphIHBvcHJhdmxqZW51IHV6b3JhY2t1IGRpc3BlcnppanUgbmFkIHV6b3Jrb20gb2JpbWEgNTAgaXogc3RhbmRhcmRuZSBub3JtYWxuZSByYXNwb2RlbGUgaSBhbmFsaXppcmF0aSByYXNwb2RlbHUgZHpla25hamYgb2NlbmUuCgpgYGB7cn0Kc2V0LnNlZWQoMTMpCiMgZ2VuZXJpc2VtbyB1em9yYWsKamFja19zbXAgPC0gcmVwbGljYXRlKDEwMDAwLCBWX2phY2socm5vcm0obiwgMSksIHZhcikpCmBgYApUZW9yaWpza2EgZGlzcGVyemlqYSBvY2VuZSBqZSAyLzQ5ID0gMC4wNDA4Ck1pIGRvYmlqYW1vCmBgYHtyfQptZWFuKGphY2tfc21wKQpgYGAKc3RvIGplIHZlY2Ugb2QgdGVvcmlqc2tlLCBzdG8gaSB0ZW9yaWpza2kgcmV6dWx0YXRpIG8gZHpla25hamYgb2NlbmkgcG9rYXp1anUuIEFsaSBwb2dsZWRham1vIG9jZW51IGd1c3RpbmUgb3ZvZyB1em9ya2E6CmBgYHtyfQpwbG90KGRlbnNpdHkoamFja19zbXApKQphYmxpbmUodiA9IDIvNDksIGNvbCA9ICJncmVlbiIpCmFibGluZSh2ID0gbWVhbihqYWNrX3NtcCksIGNvbCA9ICJyZWQiKQpgYGAKClZpZGltbyBkYSBqZSBvdmEgcmFzcG9kZWxhIGFzaW1ldHJpY25hLCBzYSB6ZWxlbm9tIHZlcnRpa2Fsbm9tIGxpbmlqb20gc21vIG96bmFjaWxpIHRlb3JpanNrdSBkaXNwZXJ6aWp1LCBhIHNhIGNydmVub20gb2Nla2l2YW5qZSAkXGhhdCBWX3tqYWNrfSQgb2NlbmUuIFVwb3JlZGltbyB2ZXJvdmF0bm9jZSBkYSBkb2JpamVtbyB2cmVkbm9zdGkgbWFuanUgb2QgdGVvcmlqc2tlIGRpc3BlcnppamUgc2Egb25pbSBkYSBkb2JpamVtbyB2ZWNlIHZyZWRub3N0aSBvZCB0ZW9yaWpza2UgZGlzcGVyemlqZS4KYGBge3J9Cm1lYW4oamFja19zbXAgPCAyLzQ5KQptZWFuKGphY2tfc21wID49IDIvNDkpCmBgYAoKVmVjYSBqZSB2ZXJvdmF0bm9jYSBkYSBkb2JpamVtbyBvY2VudSBkaXNwZXJ6aWplIGtvamEgamUgbWFuamEgb2Qgc3R2YXJuZSBkaXNwZXJ6aWplLCBpYWtvIGplIG9jZWtpdmFuamUgb2NlbmUgZGlzcGVyemlqZSB2ZWNlIG9kIHN0dmFybmUuCgpTdG9nYSBzZSBuZSBtb3plbW8gcHJldmlzZSBvc2xvbml0aSBkYSBuYW0gRHpla25hamYgZGFqZSBvY2VuZSBrb2plIHByZW1hc3VqdSB0ZW9yaWpza3UgdnJlZG5vc3QuCg==