alpha
估计
beta和
beta distribution参数。我似乎无法在网上找到这样的例子。下面我尝试了三种方法：
betareg包、Ben Bolker的
bbmle包和
optim中的
R。
betareg
包返回
phi而不是
alpha和
beta。它还返回一个
warning。我认为
bbmle包可以满足我的需求，尽管我似乎必须偷工减料
data.frame才能获得估计值。我的
optim函数没有返回良好的估计值，我不知道如何纠正它。我希望有人可以验证我对
bbmle
包所做的事情是否正确，也许可以告诉我如何纠正我的
optim功能。下面我创建了一个概率数据集，并展示了所有三种方法的
R代码。创建概率数据集：
set.seed(1234)
# generate data
n.samples <- 1000
alpha <- 2
beta <- 2
my.probs <- rbeta(n.samples, alpha, beta)
mean(my.probs)
#[1] 0.4925867
这是我使用 R
包的
betareg代码。我知道
µ = p/(p + q)和
φ = p + q，我认为
p = alpha和
q = beta。我不确定如何从 phi 的估计中获得
alpha和
beta的估计。但鉴于我有两个带有两个未知数的方程，我可能可以定义
p = φ - q并将其代入
µ的方程中。
# install.packages('betareg')
library(betareg)
my.model <- betareg(formula = my.probs ~ 1)
summary(my.model)
# Call:
# betareg(formula = my.probs ~ 1)
#
# Standardized weighted residuals 2:
# Min 1Q Median 3Q Max
# -3.8992 -0.6445 -0.0096 0.6284 3.1924
#
# Coefficients (mean model with logit link):
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.02684 0.02740 -0.98 0.327
#
# Phi coefficients (precision model with identity link):
# Estimate Std. Error z value Pr(>|z|)
# (phi) 4.1774 0.1686 24.77 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Type of estimator: ML (maximum likelihood)
# Log-likelihood: 138.9 on 2 Df
# Number of iterations: 10 (BFGS) + 2 (Fisher scoring)
# Warning message:
# In deparse(x$call, width.cutoff = floor(getOption("width") * 0.85)) :
# invalid 'cutoff' value for 'deparse', using default
exp(-0.02684) / (1 + exp(-0.02684))
#[1] 0.4932904
这是我的 R
包的
bbmle代码。我定义了
data.frame和
x的占位符
y，其中每个总是
1。这些似乎有效，但我不确定是否可以这样做。我明白了
warnings。
β 二项式和 beta 分布的 alpha 和 beta 估计
# install.packages("bbmle")
library("bbmle")
x <- rep(1, length(my.probs))
y <- rep(1, length(my.probs))
my.model <- mle2(minuslogl = my.probs ~ dbeta(shape1, shape2), start = list(shape1 = 2, shape2 = 30), data = data.frame(x, y))
# There were 12 warnings (use warnings() to see them)
summary(my.model)
# Maximum likelihood estimation
#
# Call:
# mle2(minuslogl = my.probs ~ dbeta(shape1, shape2), start = list(shape1 = 2,
# shape2 = 30), data = data.frame(x, y))
#
# Coefficients:
# Estimate Std. Error z value Pr(z)
# shape1 2.060660 0.087702 23.496 < 2.2e-16 ***
# shape2 2.116721 0.090365 23.424 < 2.2e-16 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# -2 log L: -277.8204
这是我的 optim
功能。我承认我在这里有点猜测，因为我没有
beta估计的经验。估计似乎不太好，我得到
warnings。
beta.reg = function(betas, my.probs){
b0 = betas[1]
alpha = betas[2]
beta = betas[3]
n <- length(my.probs)
llh <- rep(0, n)
for(i in 1:n){
r = my.probs[i] - b0
llh[i] = llh[i] + r^2
}
-sum(dbeta(llh, alpha, beta, log = TRUE))
}
alpha.int <- alpha+1
beta.int <- beta-1
my.model <- optim(c(mean(my.probs), alpha.int, beta.int),
beta.reg, my.probs = my.probs, method = "BFGS", hessian=TRUE)
# There were 41 warnings (use warnings() to see them)
my.model$par
#[1] 0.4901505 0.5561770 11.0689391
mean(my.probs)
my.model$value
#[1] -2170.605
covmat <- solve(my.model$hessian)
mySE <- diag(covmat)^0.5
mySE
#[1] 0.0007311783 0.0207619704 0.6066361917
感谢您对此提供的任何帮助。
alpha
、
beta包
optim和
R包
bbmle获得了几乎相同的
R和
betareg估计。但是，对于
betareg，我可能需要使用
delta method来获取
SE和
alpha的
beta。首先我生成一些数据，这次使用与我原来的帖子中使用的不同的
alpha
和
beta值。
set.seed(4321)
n.samples <- 1000
alpha <- 3
beta <- 5
my.probs <- rbeta(n.samples, alpha, beta)
mean(my.probs)
#[1] 0.3808036
这是我的 R
代码，带有
optim。我不再尝试估计
b0。我仍然明白
warnings。
beta.reg = function(betas, my.probs){
alpha = betas[1]
beta = betas[2]
n <- length(my.probs)
llh <- rep(0, n)
for(i in 1:n){
llh[i] = llh[i] + dbeta(my.probs[i], alpha, beta, log = TRUE)
}
-sum(llh)
}
alpha.int <- alpha+1
beta.int <- beta-1
my.model <- optim(c(alpha.int, beta.int),
beta.reg, my.probs = my.probs, method = "BFGS", hessian=TRUE)
#There were 50 or more warnings (use warnings() to see the first 50)
my.model$par
#[1] 3.058483 4.979886
mean(my.probs)
my.model$value
#[1] -428.0194
2 * my.model$value
#[1] -856.0388
covmat <- solve(my.model$hessian)
mySE <- diag(covmat)^0.5
mySE
#[1] 0.1308876 0.2199231
这是 R
包的
bbmle代码。我仍然得到
warnings并且我仍然必须使用占位符
data.frame。然而，
point estimates、
SE和
log-likelihood实际上与上面
optim估计的相同。
# install.packages("bbmle")
library("bbmle")
x <- rep(1, length(my.probs))
y <- rep(1, length(my.probs))
my.model <- mle2(my.probs ~ dbeta(shape1, shape2),
start = list(shape1 = 4, shape2 = 4),
data = data.frame(x, y))
# Warning messages:
# 1: In dbeta(x = c(0.295375313127293, 0.332344598431757, 0.525281773592328, :
# NaNs produced
# 2: In dbeta(x = c(0.295375313127293, 0.332344598431757, 0.525281773592328, :
# NaNs produced
# 3: In dbeta(x = c(0.295375313127293, 0.332344598431757, 0.525281773592328, :
# NaNs produced
# 4: In dbeta(x = c(0.295375313127293, 0.332344598431757, 0.525281773592328, :
# NaNs produced
summary(my.model)
# Maximum likelihood estimation
#
# Call:
# mle2(minuslogl = my.probs ~ dbeta(shape1, shape2), start = list(shape1 = 2,
# shape2 = 30), data = data.frame(x, y))
#
# Coefficients:
# Estimate Std. Error z value Pr(z)
# shape1 3.05847 0.13089 23.367 < 2.2e-16 ***
# shape2 4.97987 0.21992 22.644 < 2.2e-16 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# -2 log L: -856.0388
这是 R
的
betareg代码。我仍然对这种方法感到满意。这次我实际上将
warnings和
φ转换为
µ和
alpha但我并不试图获得它们的
beta。
SE's