Procentielscheiding

Stel dat je de steekproef genoeg wilt hebben om sommige percentielen niet te overschrijden. Hier is de code:

quantile.CI <- function(n, q, alpha=0.1) {
  #
  # Search over a small range of upper and lower order statistics for the 
  # closest coverage to 1-alpha (but not less than it, if possible).
  #
  u <- qbinom(1-alpha/2, n, q) + (-2:2) + 1
  l <- qbinom(alpha/2, n, q) + (-2:2)
  u[u > n] <- Inf
  l[l < 0] <- -Inf
  coverage <- outer(l, u, function(a,b) pbinom(b-1,n,q) - pbinom(a-1,n,q))
  if (max(coverage) < 1-alpha) i <- which(coverage==max(coverage)) else
    i <- which(coverage == min(coverage[coverage >= 1-alpha]))
  i <- i[1]
  #
  # Return the order statistics and the actual coverage.
  #
  u <- rep(u, each=5)[i]
  l <- rep(l, 5)[i]
  return(list(Interval=c(l,u), Coverage=coverage[i]))
}


###bootstrap alternative method############

quantile.CI.via.bootstrap <- function(x, p, alpha = 0.1) {
  ## Purpose:
  ##   Calculate a two-sided confidence interval with confidence level of (1 - alpha) for
  ##   a quantile, based on the (computing intensive) bootstrap resampling method. 
  ##
  ## Arguments:
  ##   - x: a vector of values, representing a data sample. 
  ##   - p: probability cutpoint for the quantile (between 0 and 1).
  ##   - alpha: type I error level (default to 0.1 so a 90% CI is calculated)
  ##
  ## Return:
  ##   - CI: the lower and upper limits of the two-sided CI. 
  
  q <- quantile(x, probs = p)         
  message("Quantile Point Estimate = ", q, " (Probability Cutpoint = ", p, ")\n")
  
  ## Bootstrap resampling with 2000 replications
  library(boot)
  set.seed(1)
  b <- boot(x, function(x, i) quantile(x[i], probs = p), R = 2000)
  
  boot.ci(b, conf = 1 - alpha, type = c("norm", "basic", "perc", "bca"))
  
}

De eerste methode implementeert de methode van Hahn & Meeker 1991 en de tweede is gewoon een niet-parametrische bootstrap.

n.s<-seq(50,200,5)

q.90<-0.90 #less strikt
q.99<-0.99
resultaten.90<-matrix(NA,2,lengte(n.s))
resultaten.99<-matrix(NA,2,lengte(n.s))

for (i in 1:length(n.s)){
  
  lu.90 <- kwantiel.CI(n.s[i], q.90)$Interval
  results.90[1,i]<-lu.90[1]
  results.90[2,i]<-lu.90[2]
  
  lu.99 <- kwantiel.CI(n.s[i], q.99)$Interval
  results.99[1,i]<-lu.99[1]
  results.99[2,i]<-lu.99[2]
}


plot(results.90[2,],ylab="Geschatte CI-limietorder",xlab="n",type="l",col="blue",axes=FALSE,pch=19)
as(2)
axis(1, at=seq_along(results.90[1,]),labels=as.character(n.s), las=2)
vak()
lijnen(resultaten.99[1,],col="paars",type="l",pch=18)

abline(v = welke(n.s==95), col="black", lwd=3, lty=2)
# Een legende toevoegen
legenda(1, 175, legenda=c("Ondergrens 99e", "Bovengrens 90e"),
       col=c("paars", "blauw"), lty=1, cex=0,8)

Zoals we kunnen zien is een steekproef n=95 voldoende volgens de bootstrapmethode.

Bootstrap resultaten:

n.s<-seq(50,200,5)
q.95<-0.90
q.99<-0.99
resultaten.95<-matrix(NA,2,lengte(n.s))
resultaten.99<-matrix(NA,2,lengte(n.s))


voor (i in 1:lengte(n.s)){

  x <- rnorm(n.s[i], mean = welke.mean.max, sd = welke.sd.max) 1TP5Daar de grootste CV
  lu.95 <-kwantiel.CI.via.bootstrap(x, 0,90, 0,1)$bca
  results.95[1,i]<-lu.95[4]
  results.95[2,i]<-lu.95[5]
  
  lu.99 <- kwantiel.CI.via.bootstrap(x, 0.99, 0.1)$bca
  results.99[1,i]<-lu.99[4]
  results.99[2,i]<-lu.99[5]
}



plot(results.95[2,],ylab="Geschatte CI-limiet",xlab="n",type="l",col="blue",axes=FALSE,pch=19,ylim=c(300,400))
as(2)
as(1, at=seq_along(results.95[1,]),labels=as.character(n.s), las=2)
vak()
lijnen(resultaten.99[1,],col="paars",type="l",pch=18)

abline(v = welke(n.s==100), col="black", lwd=3, lty=2)
# Een legende toevoegen
legenda(1,5 , legenda=c("Ondergrens 99e", "Bovengrens 90e"),
       col=c("paars", "blauw"), lty=1, cex=0,8)

Volgens de resultaten van boostrap is 100 voldoende. Merk op dat se hier BCA-type schattingen gebruiken in een bootrstrapmonster, wat een robuuste versie is.

Niet dat beide methoden verdelingsvrij, niet-parametrisch zijn!

Je kunt ze ook voor andere percentielen gebruiken. Veel plezier!