一个相关的问题是"使用二项式分布的fitdistplus中的fitdist".fitdistrplus::fitdist是一个函数，它采用单变量数据并开始猜测参数。要拟合二项式和β二项式数据，虽然是单变量数据，但也需要大小。如果每个数据的大小都是固定的，那么前面提到的链接就有所需的修复。然而，如果大小不同，需要传递向量，我不确定如何获得正常运行的调用。
下面代码中的opt_one是前面提到的链接文章中提供的解决方案，也就是说，集群大小是已知的并且是固定的。对于opt_one，我错误地指定了fix.arg=list(size=125)(本质上使N的每个元素都为125(，这已经足够接近了，并且代码运行了。然而，N中的簇大小实际上是不同的。我尝试在opt_two中指定此项，但得到一个错误。任何想法都将不胜感激。
library(fitdistrplus)
library(VGAM)
set.seed(123)
N <- 100 + rbinom(1000,25,0.9)
Y <- rbetabinom.ab(rep(1,length(N)), N, 1, 2)
head(cbind(Y,N))
opt_one <-
fitdist(data=Y,
distr=pbetabinom.ab,
fix.arg=list(size=125),
start=list(shape1=1,shape2=1)
)
opt_one
哪个给出：
> head(cbind(Y,N))
Y N
[1,] 67 123
[2,] 14 121
[3,] 15 123
[4,] 42 121
[5,] 86 120
[6,] 28 125
> opt_one <-
+ fitdist(data=Y,
+ distr=pbetabinom.ab,
+ fix.arg=list(size=125),
+ start=list(shape1=1,shape2=1)
+ )
Warning messages:
1: In fitdist(data = Y, distr = pbetabinom.ab, fix.arg = list(size = 125), :
The dbetabinom.ab function should return a zero-length vector when input has length zero
2: In fitdist(data = Y, distr = pbetabinom.ab, fix.arg = list(size = 125), :
The pbetabinom.ab function should return a zero-length vector when input has length zero
> opt_one
Fitting of the distribution ' betabinom.ab ' by maximum likelihood
Parameters:
estimate Std. Error
shape1 0.9694054 0.04132912
shape2 2.1337839 0.10108720
Fixed parameters:
value
size 125
不错，很糟糕，因为shape1和shape2参数分别为1和2，正如我们创建Y时指定的那样。这是选项2:
opt_two <-
fitdist(data=Y,
distr=pbetabinom.ab,
fix.arg=list(size=N),
start=list(shape1=1,shape2=1)
)
这给出了一个错误：
> opt_two <-
+ fitdist(data=Y,
+ distr=pbetabinom.ab,
+ fix.arg=list(size=N),
+ start=list(shape1=1,shape2=1)
+ )
Error in checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, :
'fix.arg' must specify names which are arguments to 'distr'.
首次发布后的尝试(感谢Dean Follmann(
我知道我可以编写自己的二项式似然(opt_three，如下所示(，但我真的很想使用带有fitdist对象的工具——也就是说，让opt_two工作。
library(Rfast)
loglik <-function(parm){
A<-parm[1];B<-parm[2]
-sum( Lgamma(A+B) - Lgamma(A)- Lgamma(B) + Lgamma(Y+A) + Lgamma(N-Y+B) - Lgamma(N+A+B) )
}
opt_three <- optim(c(1,1),loglik, method = "L-BFGS-B", lower=c(0,0))
opt_three
哪个给出：
> opt_three
$par
[1] 0.9525161 2.0262342
$value
[1] 61805.54
$counts
function gradient
7 7
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
还与Ben Bolker使用mle2的回答有关。fitdist解决方案仍然逍遥法外。
查看?fitdistrplus::fitdist()帮助页面的示例4：
# (4) defining your own distribution functions, here for the Gumbel distribution
# for other distributions, see the CRAN task view
# dedicated to probability distributions
#
dgumbel <- function(x, a, b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b))
pgumbel <- function(q, a, b) exp(-exp((a-q)/b))
qgumbel <- function(p, a, b) a-b*log(-log(p))
fitgumbel <- fitdist(serving, "gumbel", start=list(a=10, b=10))
summary(fitgumbel)
plot(fitgumbel)
然后——因为你实际上RTM而感到受到启发和知情——用指定的N来制作你自己的[dpq]函数：
dbbspecifiedsize <- function(x, a, b) dbetabinom.ab(x, size=N, shape1=a, shape2=b)
pbbspecifiedsize <- function(q, a, b) pbetabinom.ab(q, size=N, shape1=a, shape2=b)
qbbspecifiedsize <- function(p, a, b) qbetabinom.ab(p, size=N, shape1=a, shape2=b)
opt_four <-
fitdist(data=Y,
distr="bbspecifiedsize",
start=list(a=1,b=1)
)
opt_four
它给出：
> opt_four
Fitting of the distribution ' bbspecifiedsize ' by maximum likelihood
Parameters:
estimate Std. Error
a 0.9526875 0.04058396
b 2.0261339 0.09576709
其与CCD_ 14的估计非常相似并且是CCD_。