矢量化R函数

Question

在a与b如下所示是相同的量,但是在河两种不同的方式进行计算,他们大多是相同的,但有几个大的差异.我无法弄清楚为什么会这样.

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]

当浮点错误导致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
}