CODEWITHSUNDEEP
r
LaTeXFan
LaTeXFan

Reputation: 1231

Vectorised R functions

The a and b as shown below are the same quantities but are calculated in two different ways in R. They are mostly the same but with several big differences. I could not figure out why that was the case. 

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

One other strange thing is as follows.

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

They are supposed to be the same. But they clearly are not. Why?

Upvotes: 1

Views: 115

Answers (1)

David Robinson
David Robinson

Reputation: 78590

The problem is you are indexing by a numeric, alpha*100 in this line:

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

When floating point error causes seq(0.02, 0.99, by=0.01) to be even slightly less than the corresponding integer in 2:99, you end up extracting the wrong number from theta0.100. For example, see:

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

A quick solution is to change alpha*100 to round(alpha*100), as below, to ensure that you're always selecting by the nearest round number.

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

Upvotes: 5

Related Questions