matlab - 我们如何在 MATLAB 中找到 gamma 分布的百分位数或分位数？

Question

假设我们在 MATLAB 中有这个伽马分布：

我希望这部分分布具有更高的密度（x 轴范围）。如何在 MATLAB 中提取它？histfit我使用函数拟合这个分布。

附言。我的代码：

figure;
histfit(Data,20,'gamma');
    [phat, pci] = gamfit(Data);

phat =

   11.3360    4.2276

pci =

    8.4434    3.1281
   15.2196    5.7136

Answer 1

当您使用 将 gamma 分布拟合到数据时[phat, pci] = gamfit(Data);，phat包含 MLE 参数。您可以将其插入gaminv：

x = gaminv(p, phat(1), phat(2));

其中p是百分比向量，例如p = [.2, .8]。

Answer 2

我总是将我的代码基于以下内容。这是曾经制作的代码，我现在在必要时进行更改。也许你会发现它也很有用。

table2.10 <- cbind(lower=c(0,2.5,7.5,12.5,17.5,22.5,32.5,
                           47.5,67.5,87.5,125,225,300),upper=c(2.5,7.5,12.5,17.5,22.5,32.5,
                                                               47.5,67.5,87.5,125,225,300,Inf),freq=c(41,48,24,18,15,14,16,12,6,11,5,4,3))
loglik <-function(p,d){
  upper <- d[,2]
  lower <- d[,1]
  n <- d[,3]
  ll<-n*log(ifelse(upper<Inf,pgamma(upper,p[1],p[2]),1)-
              pgamma(lower,p[1],p[2]))
  sum( (ll) )
}

p0 <- c(alpha=0.47,beta=0.014)
m <- optim(p0,loglik,hessian=T,control=list(fnscale=-1),
           d=table2.10)

theta <- qgamma(0.95,m$par[1],m$par[2])
theta

也可以使用 Delta 方法创建 95% 的置信区间。为此，我们需要区分声明 F^(-1)_{X} (0.95; α, β) 相对于 α 和 β 的分布函数。

p <- m$par
eps <- c(1e-5,1e-6,1e-7)
d.alpha <- 0*eps
d.beta <- 0*eps
for (i in 1:3){
d.alpha[i] <- (qgamma(0.95,p[1]+eps[i],p[2])-qgamma(0.95,p[1]-eps[i],p[2]))/(2*eps[i])
d.beta[i] <- (qgamma(0.95,p[1],p[2]+eps[i])-qgamma(0.95,p[1],p[2]-eps[i]))/(2*eps[i])
}
d.alpha

d.beta

var.p <- solve(-m$hessian)
var.q95 <- t(c(d.alpha[2],d.beta[2])) %*% var.p %*% c(d.alpha[2],d.beta[2])
qgamma(0.95,p[1],p[2]) + qnorm(c(0.025,0.975))*sqrt(c(var.q95))

甚至可以对 α 和 β 的估计值使用参数引导程序，以获得损失分布的第 95 个百分位数的 B = 1000 个不同的估计值。并使用这些估计来构建 95% 的置信区间

library(mvtnorm)
B <- 10000
q.b <- rep(NA,B)
for (b in 1:B){
p.b <- rmvnorm(1,p,var.p)
if (!any(p.b<0)) q.b[b] <- qgamma(0.95,p.b[1],p.b[2])
}

quantile(q.b,c(0.025,0.975))

为了进行非参数引导，我们首先“扩展”数据以反映每个单独的观察结果。然后我们从行号中进行替换抽样，计算频率表，估计模型及其 95% 百分位数。

line.numbers <- rep(1:13,table2.10[,"freq"])
q.b <- rep(NA,B)
table2.10b <- table2.10
for (b in 1:B){
line.numbers.b <- sample(line.numbers,size=217,replace=TRUE)
table2.10b[,"freq"] <- table(factor(line.numbers.b,levels=1:13))
m.b <- optim(m$par,loglik,hessian=T,control=list(fnscale=-1),
d=table2.10b)
q.b[b] <- qgamma(0.95,m.b$par[1],m.b$par[2])
}
q.npb <- q.b
quantile(q.b,c(0.025,0.975))