probability - Why am I getting NAs in this calculation in R?

Question

question from:https://stackoverflow.com/questions/65876352/why-am-i-getting-nas-in-this-calculation-in-r

Answer 1

Bernhard's answer correctly identifies the problem: If logit is large, exp(logit) = Inf. Here is a solution:

for(i in 1:dim(beta.grd)[1]){ # iterate through 600 possible beta values in beta grid
    
    beta.ind = 0 # indicator for current pair of beta values
    
    for(j in 1:N){ # iterate through all possible Nsums
        logit = beta.grd[i,1]/N*(j - .1*N)^2 + beta.grd[i,2];
        ## This one isn't great because exp(logit) can be very large
        # phi01 = exp(logit)/(1 + exp(logit))
        ## So, we say instead
        ## phi01 = 1 / ( 1 + exp(-logit) )
        phi01 = plogis(logit)
        
        
        if(is.na(phi01)){ 
            count = count + 1
        }
    }
}

cat("Total number of invalid probabilities: ", count)
# Total number of invalid probabilities:  0

We can use the more stable 1 / (1 + exp(-logit) (to convince yourself of this, multiply your expression with exp(-logit) / exp(-logit)), and luckily either way, R has a builtin function plogis() that can calculate these probabilities quickly and accurately. You can see from the help file (?plogis) that this function evaluates the expression I gave, but you can also double check to assure yourself

x = rnorm(1000)
y = 1 / (1 + exp(-x))
z = plogis(x)
all.equal(y, z)
[1] TRUE