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

Question

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

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

PS。我的代码：

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

当您使用 [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))

甚至可以在 α 和 β 的估计值上使用参数 bootstrap 以获得 B = 1000 个不同的估计值损失分布的第 95 个百分位数。并使用这些估计值构建一个 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))

为了进行非参数化 bootstrap，我们首先“扩展”数据以反映每个单独的观察结果。然后我们从行号中进行替换抽样，计算频率 table，估计模型及其 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))

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

How can we find percentile or quantile of gamma distribution in MATLAB?

matlab

distribution

gamma-distribution