提问人：Aku-Ville Lehtimäki 提问时间：10/3/2023 更新时间：10/4/2023 访问量：69
在 R 中，一种巧妙的方法来获得平均值以及给定字符串格式的概率分布的 5 个和 95 个百分位数？
A neat way in R to get mean and 5 and 95 percentiles given a probability distribution in a string format?
问：
我想编写一个将字符串作为输入的函数。确切地说，此输入表示概率分布，输出为 5 个百分位数、平均值和 95 个百分位数作为向量或列表。
例如，如果我输入：“beta（0.2,0.3）” 我会得到：2.793073e-06 0.4 0.9992341
首先，要做的是从字符串中读取分布，这可以通过正则表达式来完成：
dist <- gsub("[^A-z]","","beta(0.2,0.3)")
params <- gsub("[^0-9.,]","","beta(0.2,0.3)")
可以使用 strsplit 将参数放入数值向量中，例如
params <- as.numeric(unlist(strsplit(gsub("[^0-9.,]","","beta(0.2,0.3)"),split=",")))
现在，我需要声明一个均值函数（因为根据我的理解，R 中不存在随机分布均值函数）。对于 beta 分发，这将是：
beta_mean <- function(alpha, beta) {
return(alpha / (alpha + beta))
}
我可以从 qbeta 函数中获得百分位数，即：
qbeta(c(0.05,0.95), params[1], params[2])
既然，我想处理许多不同的发行版，有没有比以下更优雅的方法：
meanvalue <- NA
percentiles <- NA
if(dist == "beta") {
meanvalue <- beta_mean(params[1], params[2])
percentiles <- qbeta(c(0.05,0.95), params[1], params[2])
} else if (dist == "gamma") {
meanvalue <- gamma_mean(params[1], params[2])
percentiles <- qgamma(c(0.05,0.95), params[1], params[2])
} #and so on and so on and so on...
return(c(percentiles[1],meanvalue,percentiles[2])
所以，我想做的是将发行版名称字符串（例如“beta”或“gamma”或其他）链接到相应的函数（DISTRIBUTIONNAME_mean 和 qDISTRIBUTIONNAME），这样我就不必使用那个非常长的 if-else-结构，它包含太多（不必要的？）重复。
我怎样才能做到这一点？
答：
2赞
Michael M
10/3/2023
#1
我只能给出涵盖分位数的部分答案：
a_string <- "beta(0.2, 0.3)"
eval(parse(text = paste0("q", sub("\\(", "(c(0.05, 0.95), ", a_string))))
# 2.793073e-06 9.992341e-01
为了得到平均值，你需要参数化所有分布的期望值，这很痛苦。
2赞
jpsmith
10/3/2023
#2
您可以在这里一起使用：eval(parse(...))paste0
distfun <- function(x){
par <- unlist(stringr::str_extract_all(a_string, "\\d+\\.\\d"))
m <- paste0(gsub("[^a-zA-Z]", "", x), "_mean(", par[1], ",", par[2], ")")
setNames(c(
eval(parse(text = m)),
eval(parse(text = paste0("q", sub("\\(", "(c(0.05, 0.95), ", x))))
), c("mean", "q05", "q95"))
}
strng <- "beta(0.2, 0.3)"
distfun(strng)
# mean q05 q95
# 4.000000e-01 2.793073e-06 9.992341e-01
或者，如果您希望按照问题的确切顺序进行处理：
distfun <- function(x){
par <- unlist(stringr::str_extract_all(a_string, "\\d+\\.\\d"))
res <- setNames(c(
eval(parse(text = paste0(gsub("[^a-zA-Z]", "", x), "_mean(", par[1], ",", par[2], ")"))),
eval(parse(text = paste0("q", sub("\\(", "(c(0.05, 0.95), ", x))))
), c("mean", "q05", "q95"))
res[c(2,1,3)]
}
strng <- "beta(0.2, 0.3)"
distfun(strng)
# q05 mean q95
# 2.793073e-06 4.000000e-01 9.992341e-01
0赞
Aku-Ville Lehtimäki
10/4/2023
#3
我今天想出了这个解决方案。首先，我们可以获取统计包中的每一个分布：
library(dplyr)
fnnames <- ls("package:stats") %>% #Get the list of function names in stats package
grep("^d|^p|^q|^r", ., value = TRUE) %>% #Exclude names NOT starting with "d", "p", "q" or "r"
substring(2) %>% #Remove first letter (e.g. rnorm, pnorm, qnorm, and rnorm become all just norm)
table %>% #Create a frequency table, a vector values are frequencies
.[. == 4] %>% #Since distribution related functions come in quadruplets, we keep only quadruplets
names(.) #Finally, we are interested in the function names, not the frequencies any more.
现在，我们知道有哪些功能可用。让我们创建一个查找字典（这是 R，所以它实际上是一个命名列表，尽管具有与字典相同的功能）：
distlist <- lapply(fnnames, function(fnname){
list( #a list inside of dist list
"q" = get(paste0("q",fnname)), #these sub-lists contain "key" q for quantile function. However, I am not sure whether or not "get" function is doing the same as eval(parse(text=...)) under the hood.
"mean" = NULL #I'll return to this later.
)
}) %>% setNames(fnnames) #Let's set names for the sub-lists
现在，我们像这样使用我们的查找字典：
distlist[["beta"]][["q"]](c(0.05, 0.95), params[1], params[2])
[1] 2.793073e-06 9.992341e-01
现在，均值函数。大多数随机分布都有均值解析解，这更容易使用。不幸的是（据我所知），不存在这样的函数（meannorm、meanbeta、meangamma 等，将是一个很好的补充）。但是，我们当然可以集成：
integrate(function(x)(x*dbeta(x, params[1], params[2])), lower = 0, upper = 1)
0.4 with absolute error < 1.9e-05
这很慢，并且 d 函数不喜欢（即它们返回错误），如果提供了超出其支持范围的值。但是，这可以通过 try-catch 来规避，但由于期望值的分析形式大多非常简单，我更喜欢它们。
也许，我会为不同的随机分布的不同时刻（mean 是第一个时刻）编写和发布一个包，如果没有人这样做的话。
评论