给定 R 中的分位数,在 qbeta() 中求解 shape1 和 shape2
Solving for shape1 and shape2 in qbeta() given its quantiles in R
背景:
qbeta(p, shape1, shape2)
是内置的基本 R 函数。在此函数中,shape1
和 shape2
始终 > 0(出于我的目的,始终 > 1)。此外,p
是概率,因此 0 <= p <= 1.
问题:
假设我知道 qbeta(p = c(.025, .975), shape1 = x, shape2 = y ) = c(.3, .8 )
,那么,有没有一种方法可以找出 x 和 y 可以合理是(即在一定误差范围内)?
注意:我相信HERE中提供了一些相关的注释和代码。因此,如果其中一位同事能够理解或使用这些代码和注释并将其应用到此处,我将不胜感激。
最后,我的最终目标是用这个 x & y[=37= 包装一个函数 ] 寻找过程。
这是一个暴力破解的方法。我希望其他人可以想出更聪明的方法来做到这一点。
d = expand.grid(shape1 = 1:100/10, shape2 = 1:100/10)
for (i in 1:NROW(d)){
x = round(qbeta(p = c(0.025, 0.975), shape1 = d[i,1], shape2 = d[i,2]),2)
if(x[1] > 0.28 & x[1] < 0.32 & x[2] > 0.78 & x[2] < 0.82){
print(d[i,])
break
}
}
# shape1 shape2
#5368 6.8 5.4
qbeta(p = c(.025, .975), shape1 = 6.8, shape2 = 5.4)
#[1] 0.2860268 0.8109389
我在下面提供whuber's solution from here
################################
#REQUIRED FUNCTIONS
# Logistic transformation of the Beta CDF.
f.beta <- function(alpha, beta, x, lower=0, upper=1) {
p <- pbeta((x-lower)/(upper-lower), alpha, beta)
log(p/(1-p))
}
# Sums of squares.
delta <- function(fit, actual){
sum((fit-actual)^2)
}
# The objective function handles the transformed parameters `theta` and
# uses `f.beta` and `delta` to fit the values and measure their discrepancies.
objective <- function(theta, x, prob, ...) {
ab <- exp(theta) # Parameters are the *logs* of alpha and beta
fit <- f.beta(ab[1], ab[2], x, ...)
return (delta(fit, prob))
}
################################
#INPUT
x <- c(.3, .8)
x.p <- (function(p) log(p/(1-p)))(c(0.025, 0.975))
#SOLVE
start <- log(c(1e1, 1e1))
sol <- nlm(objective, start, x=x, prob=x.p, lower=0, upper=1,
typsize=c(1,1), fscale=1e-12, gradtol=1e-12)
parms <- exp(sol$estimate)
#ANSWER
parms
#[1] 7.597429 6.016240
#TEST
qbeta(p = c(.025, .975), parms[1], parms[2])
#[1] 0.3 0.8
背景:
qbeta(p, shape1, shape2)
是内置的基本 R 函数。在此函数中,shape1
和 shape2
始终 > 0(出于我的目的,始终 > 1)。此外,p
是概率,因此 0 <= p <= 1.
问题:
假设我知道 qbeta(p = c(.025, .975), shape1 = x, shape2 = y ) = c(.3, .8 )
,那么,有没有一种方法可以找出 x 和 y 可以合理是(即在一定误差范围内)?
注意:我相信HERE中提供了一些相关的注释和代码。因此,如果其中一位同事能够理解或使用这些代码和注释并将其应用到此处,我将不胜感激。
最后,我的最终目标是用这个 x & y[=37= 包装一个函数 ] 寻找过程。
这是一个暴力破解的方法。我希望其他人可以想出更聪明的方法来做到这一点。
d = expand.grid(shape1 = 1:100/10, shape2 = 1:100/10)
for (i in 1:NROW(d)){
x = round(qbeta(p = c(0.025, 0.975), shape1 = d[i,1], shape2 = d[i,2]),2)
if(x[1] > 0.28 & x[1] < 0.32 & x[2] > 0.78 & x[2] < 0.82){
print(d[i,])
break
}
}
# shape1 shape2
#5368 6.8 5.4
qbeta(p = c(.025, .975), shape1 = 6.8, shape2 = 5.4)
#[1] 0.2860268 0.8109389
我在下面提供whuber's solution from here
################################
#REQUIRED FUNCTIONS
# Logistic transformation of the Beta CDF.
f.beta <- function(alpha, beta, x, lower=0, upper=1) {
p <- pbeta((x-lower)/(upper-lower), alpha, beta)
log(p/(1-p))
}
# Sums of squares.
delta <- function(fit, actual){
sum((fit-actual)^2)
}
# The objective function handles the transformed parameters `theta` and
# uses `f.beta` and `delta` to fit the values and measure their discrepancies.
objective <- function(theta, x, prob, ...) {
ab <- exp(theta) # Parameters are the *logs* of alpha and beta
fit <- f.beta(ab[1], ab[2], x, ...)
return (delta(fit, prob))
}
################################
#INPUT
x <- c(.3, .8)
x.p <- (function(p) log(p/(1-p)))(c(0.025, 0.975))
#SOLVE
start <- log(c(1e1, 1e1))
sol <- nlm(objective, start, x=x, prob=x.p, lower=0, upper=1,
typsize=c(1,1), fscale=1e-12, gradtol=1e-12)
parms <- exp(sol$estimate)
#ANSWER
parms
#[1] 7.597429 6.016240
#TEST
qbeta(p = c(.025, .975), parms[1], parms[2])
#[1] 0.3 0.8