伽马分布参数的广义置信区间及其覆盖率的R代码

R代码，关于伽马分布中参数的广义置信区间及其覆盖率，论文不好复制这里，广义置信区间的代码都可以参考一下。

GCIab<- function(a,b,n,nr,alpha){
  y <- rgamma(n,shape = a,scale = b)
  x1<- y^(1/3)
  x1bar <- mean(x1)
  x1sd <- sd(x1)
  chi <- rchisq(nr,n-1)
  z <- rnorm(nr,0,1)
  Gu <- x1bar+z*sqrt(n-1)*x1sd/(sqrt(chi)*sqrt(n))
  Gsigma_2 <- (n-1)*(x1sd^2)/chi
  Ga=((1+0.5*Gu^2/Gsigma_2)+((1+0.5*Gu^2/Gsigma_2)^2-1)^(1/2))/9
  Ga<-sort(Ga)
  ll<-nr*(alpha/2)
  ul<-nr*(1-alpha/2) 
  Gall<-Ga[ll]
  Gaul<-Ga[ul]
  GCIa <- c(Gall,Gaul)
  Gb <- 27*sqrt(Ga)*Gsigma_2^(3/2)
  Gb<-sort(Gb)
  ll<-n*(alpha/2)
  ul<-n*(1-alpha/2) 
  Gbll<-Gb[ll]
  Gbul<-Gb[ul]
  GCIb <- c(Gbll,Gbul)
  return(c(GCIa,GCIb))
}
GCIab(2,1,1000,1000,0.05)
#覆盖率
coverage <-function(a,b,n,nr,alpha,m1,m2){
  p=0
  p1=0
  p2=0
  q=0
  q1=0
  q2=0
  #a
  for(i in 1:m1){
    CIa=GCIab(a,b,n,nr,alpha)
    l <- CIa[1]
    u <- CIa[2]
    if(l<a&a<u){
      p=p+1
    } else
      if(a<l){
        p1=p1+1
      } else{
        p2=p2+1
      }
  }
  acoverage <-c(p/m1,p1/m1,p2/m1)
  #b
  for(i in 1:m2){
    CIb=GCIab(a,b,n,nr,alpha)
    h <- CIb[3]
    j <- CIb[4]
    if(l<=h&b<=j){
      q=q+1
    } else
      if(b<l){
        q1=q1+1
      } else{
        q2=q2+1
      }
  }
  bcoverage <- c(q/m2,q1/m2,q2/m2)
  return (c(acoverage, bcoverage))
}
coverage(2,1,1000,1000,0.05,10000,10000)