从现有向量生成 beta-binomial 分布

Question

是否有可能 to/how 我可以从现有向量生成 beta-binomial 分布吗？

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

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

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

从现有向量生成 beta-binomial 分布

Generate beta-binomial distribution from existing vector

r

distribution