-1

Nether 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。

4

1 回答 1

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 函数和这些替换函数的数学运算应该是准确的。希望您能接受这项工作并使其适应您的需求。祝你好运！

于 2013-05-24T19:27:33.613 回答