将 beta 正态分布拟合到数据

Question

我想将 beta 正态分布拟合到我的数据中

X=c(5.20 , 6.80, 11.00, 21.00 ,25.50, 28.50, 30.90 ,30.90, 27.20, 17.70 ,10.50 , 6.70, 5.00,  8.00 ,14.30, 20.90 ,24.10 ,28.40 ,29.80, 30.80 ,26.80 ,20.50, 12.50 , 9.30, 19.20 , 5.60,  9.20  ,1.80 ,15.80 ,22.20 ,27.90, 30.60 ,31.10 ,28.80 ,23.30, 13.40, 4.30 , 6.80  ,7.20 ,10.30 ,17.00, 21.20 ,27.40 ,32.10, 30.20 ,25.50, 22.20, 11.30, 6.00  ,6.60 , 9.80 ,14.20 ,16.10 ,22.00 ,29.30 ,29.30 ,31.30, 26.20, 20.30 ,13.30, 5.30  ,5.00,  4.80 ,13.00 ,18.90 ,22.40, 28.30, 32.40 ,30.20 ,27.20 ,21.30, 12.00, 10.20  ,3.00 , 9.80 ,14.10 ,19.30, 24.80, 30.30 ,31.40 ,31.60, 26.40, 31.60, 11.50, 4.00 , 4.70  ,7.60, 10.00 ,17.20 ,24.20, 29.10 ,30.90 ,30.40 ,26.50, 20.00, 14.60, 5.80, -2.00 , 5.00 ,17.00, 20.70 ,23.90, 28.80 ,31.60 ,30.20 ,27.30 ,20.70 ,11.70, 7.40 , 5.40  ,9.10 ,13.80 ,14.20 ,23.70 ,26.80 ,31.70 ,29.80 ,24.40 ,19.40, 11.90, 7.40 , 8.90 , 8.90 ,14.50 ,18.01 ,23.20 ,30.00 ,32.00 ,29.10 ,25.80 ,23.30 ,13.00, 9.80)

我使用 fitdistplus 分布，但它没有这个分布，想定义它的密度和分布。我定义密度如下

mu=mean(data)
sigma=sd(data)
Phi=pnorm((x-mu)/sigma)
phi=dnorm((x-mu)/sigma)
dbetanorm=function(a,b){
((gamma(a+b))/(sigma*gamma(a)*gamma(b)))*phi*(Phi)^(a-1)((1-Phi)^(b-1))
}

但我无法定义 pbetanorm。它的公式是

F(x)=I_Φ((x-μ)/σ)  (a,b)    

有谁知道如何在 R 中做到这一点？ 感谢您的帮助

Answer 1

最简单的方法是在VGAM 包中使用pbetanorm(...)。文档here。

library(fitdistrplus)
library(VGAM)

set.seed(1)       # for reproducible example
X <- rbetanorm(1000,2,6)
params <- fitdist(X,distr=dbetanorm,start=list(shape1=1,shape2=1),
                  fix.arg=list(mean=0,sd=1))
params
# Fitting of the distribution ' betanorm ' by maximum likelihood 
# Parameters:
#        estimate Std. Error
# shape1 2.010134 0.08380233
# shape2 5.980691 0.27065835

# always plot the results...
hist(X,freq=F,xlim=c(-3,1),ylim=c(0,1),breaks=20)
x <- seq(-3,1,.1)
lines(x,dbetanorm(x,params$estimate[1],params$estimate[2]),col="red",lty=2)

编辑对 OP 包含数据的回应。

这是一个很好的对象课程，在问题中包含您的数据的重要性！！！

所以这里发生了几件事。首先也是最重要的，拟合可能很差，因为您的数据不是 beta 正态分布的。你为什么相信它？

其次，beta 正态分布有 4 个参数，shape1、shape2、mean 和 sd。在我的示例中，我修复了后两个以显示它是如何完成的，并且因为您似乎在您的问题中假设了值。在实际情况下，拟合所有 4 个参数可能很重要。

第三，fitdistr(...) 使用optim(...) 函数进行实际拟合。 optim(...) 使用局部最小化，这意味着如果您的初始参数估计值不合理地接近真实值，您可能会得到一个“拟合”（局部最小值），这是对您的数据的不良表示。经过一番实验，很明显，对于您的数据集，shape1

params <- fitdist(X,distr=dbetanorm,
                  start=list(shape1=.001,shape2=.01,mean=mean(X),sd=1))
params
# Fitting of the distribution ' betanorm ' by maximum likelihood 
# Parameters:
#           estimate  Std. Error
# shape1  0.01898767 0.010817437
# shape2  0.01562771 0.009071755
# mean   17.36844300 0.870978734
# sd      0.89752769 0.265374668


hist(X,freq=F,xlim=c(-5,35),ylim=c(0,.08),breaks=20)
x <- seq(min(X),max(X),.1)
with(params,lines(x,dbetanorm(x,estimate[1],estimate[2],estimate[3],estimate[4]),
                  col="red",lty=2))

给出这个：

这仍然不是很好。

如果没有更多的工作，我倾向于说您的数据不是 beta 正态分布的。可能是两个分布的总和？