从现有向量生成 beta 二项分布

Question

我是否可以/如何从现有向量生成 beta 二项分布？

我的最终目标是从以下数据生成 beta 二项分布，然后获得该分布的 95% 置信区间。

我的数据是兽医记录的身体状况评分。身体状况的值范围为 0-5，增量为 0.5。有人建议我 here 我的数据遵循 beta 二项分布，离散值具有有限的范围。

set1 <- as.data.frame(c(3,3,2.5,2.5,4.5,3,2,4,3,3.5,3.5,2.5,3,3,3.5,3,3,4,3.5,3.5,4,3.5,3.5,4,3.5))
colnames(set1) <- "numbers"

我看到似乎有多个函数可以做到这一点，VGAM 中的betabinomial() 和emdbook 中的rbetabinom()，但我的统计数据和编码知识还不足以理解并实施功能帮助页面上提供的说明，至少还没有以对我的预期目的有所帮助的方式。

Answer 1

我们可以看看你的变量分布，y轴是概率：

x1 = set1$numbers*2
h = hist(x1,breaks=seq(0,10))
bp = barplot(h$counts/length(x1),names.arg=(h$mids+0.5)/2,ylim=c(0,0.35))

您可以尝试拟合它，但您的数据点太少，无法估计 beta 二项式所需的 3 个参数。因此，我修复了概率，以便平均值是您分数的平均值，并且查看上面的分布似乎没问题：

library(bbmle)
library(emdbook)
library(MASS)

mtmp <- function(prob,size,theta) {
-sum(dbetabinom(x1,prob,size,theta,log=TRUE))
}

m0 <- mle2(mtmp,start=list(theta=100),
data=list(size=10,prob=mean(x1)/10),control=list(maxit=1000))
THETA=coef(m0)[1]

我们也可以使用正态分布：

normal_fit = fitdistr(x1,"normal")
MEAN=normal_fit$estimate[1]
SD=normal_fit$estimate[2]

绘制它们两个：

lines(bp[,1],dbetabinom(1:10,size=10,prob=mean(x1)/10,theta=THETA),
col="blue",lwd=2)
lines(bp[,1],dnorm(1:10,MEAN,SD),col="orange",lwd=2)
legend("topleft",c("normal","betabinomial"),fill=c("orange","blue"))

我认为您实际上可以使用正常估计，在这种情况下它将是：

normal_fit$estimate
    mean       sd 
6.560000 1.134196