使用R计算与excel中的公式相同[关闭]

Question

下界 R 和一般统计数据都是我的强项。因此，作为客户构建的一部分，他们要求在一些图表中绘制二项式数据，并为我提供了一些公式，用于获得他们期望的结果。

这些公式在 excel 中，因为我使用的是 CentOS VPS，所以我已经安装了 R 来为我做这件事，但我一直找不到正确的函数来使用。

当前已安装，最新的 R 实验室和 Binom 包，我使用 Rscript 进行计算，使用 PHP 和 Pchart 生成实际图形。

要绘制的数据是4条二项式曲线，alpha分别为0.9995、0.0005、0.995和0.005，n为X轴上的位置

在这种情况下（将单元格编号替换为 sudo 变量以使其更易于阅读）：

start = 1
xaxis = 0 (increments)
p = 0.01
alpha = 0.005
B = 1.7

excel公式：

n = start+ROUNDDOWN(xaxis*EXP(1)^B,0)
critB = CRITBINOM(n,p,alpha)
Low-adj = critB-(BINOMDIST(critB,n,p,TRUE)-alpha)/(BINOMDIST(critB,n,p,TRUE)-BINOMDIST(critB-1,n,p,TRUE))
Low Alert = IF(ISERROR(100*Low-adj/n)=TRUE,"",Low-adj/n)

上面应该什么都不返回并且应该继续这样做，直到 xaxis = 14，其中预期的结果是：Low Alert = 0.000156

有没有人可以帮助编写一个 Rscript 来处理这个问题？我使用的是 binom.confinit() 但客户端现在返回了上述内容，因此我需要将 binom.confinit 替换为可以实现此目的的脚本。

值得注意的是，这是网站的一部分，而不是一次性的，因此是 php / pchart 而不是 gnumeric。

Answer 1

我没有得到 OP 从公式中报告的确切答案，因此我为 Excel 等效项编写了一系列直接替换函数，并编写了如何在 R 中使用它们的示例。这不是R 最有效的方法，但它可能是 OP 走向实现的最方便的方法。

CRITBINOM：这本质上是一个计数函数。它需要二项式试验的规模、成功的概率和 alpha 值。它返回累积概率大于给定 alpha 值的最小样本量。

CRITBINOM <- function(.trials, .probability_s, .alpha){
    i <- 0
    while(sum(dbinom(0:i, .trials, .probability_s)) < .alpha){
        i <- i + 1
        print(i)
    }
    return(i)
}

BINOMDIST：在 Excel 中，这实际上是两个带有布尔开关的函数。如果开关为 TRUE，则函数返回给定数量的二项式成功的左尾累积分布值给定试验规模和成功概率。如果开关为假，则该函数返回给定相同信息的概率质量函数（成功次数的概率）。

BINOMDIST <- function(.number_s, .trials, .probability_s, .cumulative){
    if(.cumulative){
        return(sum(dbinom(0:.number_s, .trials, .probability_s)))
    }else{
        return(choose(.trials,.number_s)*.probability_s^.number_s*(1-.probability_s)^(.trials-.number_s))
    }
}

ISERROR：在这种情况下，函数实际上只是检查函数的结果是否为无限（未定义）。我不会为这种特定用途复制所有 Excel 功能。

ISERROR <- function(.value){
    return(is.infinite(.value))
}

ROUNDDOWN：这是那些奇怪的 Excel 小函数之一。它四舍五入，但它只向下四舍五入。在这种情况下，我们实际上并没有四舍五入，而是通过乘以 10^digits、删除任何余数然后除以 10^digits 来截断数字。

ROUNDDOWN <- function(.number, .num_digits){
    num_digits <- as.integer(.num_digits)
    return(as.integer(.number*10^num_digits)/(10^num_digits))
}

示例 R 代码：现在我将展示如何使用所有这些来复制 OP 的 Excel 任务。首先，我将定义一个矢量化函数来一次计算所有的“n”值。

n <- function(.start, .increments, .B){
    return(.start + ROUNDDOWN(.increments * exp(1)^.B, 0))
}

接下来，我创建一个函数来确定单个低警报值。这是 OP 的大部分工作都包含在其中的地方。函数看起来应该几乎完全相同。

generate_Low_Alert <- function(.n, .probability_s, .alpha){
    critB <- CRITBINOM(.n, .probability_s, .alpha)
    Low_adj <- critB-(BINOMDIST(critB, .n, .probability_s,TRUE)-.alpha)/(BINOMDIST(critB, .n, .probability_s,TRUE)-BINOMDIST(critB-1, .n, .probability_s,TRUE))
    if(ISERROR(100 * Low_adj / .n)){
        return("")
    }else{
        return(Low_adj/.n)
    }
}

最后，我做了一个包装来填满整个混乱。

generate_data <- function(.B, .probability_s, .alpha, .start, .increments){
    Low_Alerts <- integer(length(.increments))
    n_values <- n(.start, .increments, .B)

    for(i in 1:length(n_values)){
        Low_Alerts[i] <- generate_Low_Alert(n_values[i], .probability_s, .alpha)
    }
    return(Low_Alerts)
}

我基本上只是循环遍历每个“n”值并生成警报（“”或 Low_adj/n 的值）。所有这些都存储在一个数组中，并作为函数的结果返回。

要全部使用，我会调用包装函数，如下所示：

generate_data(1.7, 0.01, 0.005, 1, 0:100)

现在，这与 Excel 方法略有不同，因为我们在开始时定义了整套步骤 (0:100)，而不是一次一个。否则，这将完全复制我构建的 Excel 版本。

免责声明：我无法获得与 OP 相同的结果（在 x_axis=97 处看到警报，而不是 14），但 Excel 函数和这些替换函数的数学运算应该是准确的.希望您能接受这项工作并使其适应您的需求。祝你好运！