向量化的 R 函数

Question

如下所示的a 和b 是相同的数量，但在R 中以两种不同的方式计算。它们大部分相同，但有几个很大的不同。我不明白为什么会这样。

theta0 <- c(-0.4, 10)

OS.mean <- function(shape, rank, n=100){
  term1 <- factorial(n)/(factorial(rank-1)*factorial(n-rank))
  term2 <- beta(n-rank+1, rank) - beta(n-rank+shape+1, rank)
  term1*term2/shape
}

OS.mean.theta0.100 <- OS.mean(theta0[1], rank=seq(1, 100, by=1))

Bias.MOP <- function(shape, scale, alpha){
  scale*shape*OS.mean.theta0.100[alpha*100]/(1-(1-alpha)^shape) - scale
}

a <- rep(0, 98)
for(i in 2:99){
  a[i-1] <- Bias.MOP(theta0[1], theta0[2], i/100)
}
plot(a)

b <- Bias.MOP(theta0[1], theta0[2], seq(0.02, 0.99, by=0.01))
plot(b)

a-b

另外一件奇怪的事情如下。

b[13] # -0.8185083
Bias.MOP(theta0[1], theta0[2], 0.14) # -0.03333929

它们应该是相同的。但他们显然不是。为什么？

Answer 1

问题是您在这一行中按数字 alpha*100 进行索引：

OS.mean.theta0.100[alpha*100]

当floating point error 导致seq(0.02, 0.99, by=0.01) 甚至略小于2:99 中的相应整数时，您最终会从theta0.100 中提取错误的数字。例如，参见：

x <- 1:10
x[5]
# [1] 5
x[6]
# [1] 6
x[5.99999999]
# [1] 5

一种快速的解决方案是将alpha*100 更改为round(alpha*100)，如下所示，以确保您始终选择最接近的整数。

Bias.MOP <- function(shape, scale, alpha){
  scale*shape*OS.mean.theta0.100[round(alpha*100)]/(1-(1-alpha)^shape) - scale
}