拟合β正态分布到数据

Question

我想将betanormal分配到我的数据中

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.the的公式是

F(x)=I_?((x-?)/?)  (a,b)    

有谁知道如何在R中这样做？谢谢你的帮助

Answer 1

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

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)

EDIT对OP包含数据的回应.

对于将数据包含在问题中的重要性,这是一个很好的对象课程!

所以这里有几件事情.首先也是最重要的,因为您的数据不是正常分布的,所以拟合可能很差.为什么你相信它？

其次,betanormal分布有4个参数,shape1,shape2,mean,和sd.在我的例子中,我修复后两个以显示它是如何完成的,因为你似乎在你的问题中假设值.在实际情况中,适合所有4个参数可能很重要.

三,fitdistr(...)使用optim(...)功能进行实际拟合.optim(...)使用局部最小化,这意味着如果您的初始参数估计值不合理地接近真实值,您可能会得到一个"拟合"(局部最小值),这对您的数据来说代表性很差.经过实验后,很明显,对于您的数据集,shape1 << 1和shape2 <1.因此,重新拟合更好的初始估计:

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

给出这个:

这仍然不是很好.

如果没有更多的工作,我倾向于说你的数据不是正常分布的.它可能是两个发行版的总和吗？