النسب المئوية للفصل بين النسب المئوية

لنفترض أنك ترغب في أن تكون العينة كافية لبعض النسب المئوية حتى لا تتخطى بعض النسب المئوية. إليك الرمز:

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"))
  
}

الطريقة الأولى تطبق طريقة هان وميكر 1991، والطريقة الثانية هي مجرد تمهيد غير بارامترى.

ن.س<-مجموعة تسلسل تسلسلي (50,200,5)

س.90 <- 0.90 # أقل صرامة
q.99<-0.99
results.90<-مصفوفة(NA,2,length(n.s))
results.99<-مصفوفة(NA,2,length(n.s))

بالنسبة إلى (i في 1:الطول (n.s)) {
  
  lu.90 <- Quantile.CI(n.s[i]، q.90)$Interval
  results.90[1,i]<-lu.90[1]
  results.90[2,i]<-lu.90[2]
  
  lu.99 <- Quantile.CI(n.s[i]، q.99)$Interval
  results.99[1,i]<-lu.99[1]
  results.99[2,i]<-lu.99[2]
}


رسم بياني(results.90[2,],ylab="ترتيب حد CI التقديري",xlab="n",type="l",col="blue",axes=FALSE,pch=19)
المحور(2)
محور(1, at=seq_along(results.90[1,]),labels=as.character(n.s), las=2)
صندوق()
خطوط(results.99[1,],col="أرجواني",type="l",pch=18)

أبلاين(v = الذي(n.s==95)، اللون="أسود"، lwd=3، lty=2)
# إضافة وسيلة إيضاح
وسيلة الإيضاح(1, 175, وسيلة الإيضاح=ج("الحد الأدنى 99"، "الحد الأعلى 90"),
       اللون=c("أرجواني"، "أزرق")، lty=1، cex=0.8)

كما نرى أن العينة n=95 كافية وفقاً لطريقة التمهيد.

نتائج التمهيد:

ن.س <-التسلسل الهرمي(50,200,5)
q.95<-0.90
q.99<-0.99
النتائج.95<-مصفوفة(NA,2,2,length(n.s))
results.99<-مصفوفة(NA,2,length(n.s))


بالنسبة إلى (i في 1:طول (n.s)) {

  x <- rnorm(n.s[i]، المتوسط = which.mean.max، sd = which.sd.max) 1TP5هناك أكبر سيرة ذاتية
  lu.95 <-quantile.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 <- Quantile.CI.via.bootstrap(x, 0.99, 0.1)$bca
  results.99[1,i]<-lu.99[4]
  results.99[2,i]<-lu.99[5]
}



رسم بياني(results.95[2,],ylab="حد CI المقدر",xlab="n",type="l",col="blue",axes=FALSE,pch=19,ylim=c(300,400))
محور(2)
المحور(1, at=seq_along(results.95[1,]),labels=as.character(n.s), las=2)
صندوق()
خطوط(results.99[1,],col="أرجواني",type="l",pch=18)

أبلاين(v = الذي(n.s==100)، اللون="أسود"، lwd=3، lty=2)
# إضافة وسيلة إيضاح
وسيلة الإيضاح(1,5، وسيلة الإيضاح=ج("الحد الأدنى 99"، "الحد الأعلى 90"),
       اللون=c("أرجواني"، "أزرق")، lty=1، cex=0.8)

وفقًا لنتائج boostrap 100 كافية. لاحظ أن استخدم هنا نوع BCA من التقديرات في عينات التمهيد المعززة وهي نسخة قوية.

لا يعني ذلك أن كلتا الطريقتين خاليتين من التوزيع، وغير بارامترية!

يمكنك استخدامها لنسب مئوية أخرى أيضاً. استمتع بوقتك!