为什么我在 R 的这个计算中得到 NA？

Question

在处理 Rcpp 程序时，我使用了 sample() 函数，它给了我以下错误：“NAs not allowed in probability.”我将此问题追溯到我使用的概率向量中包含 NA 值的事实。我不知道怎么做。下面是一些捕获错误的 R 代码：

n.0=20
n.1=20
n.reps=1
beta0.vals=rep(seq(-.3,.1,,n.0),n.reps)
beta1.vals=rep(seq(-7,0,,n.1),n.reps)
beta.grd=as.matrix(expand.grid(beta0.vals,beta1.vals))

n.rnd=200
beta.rnd.grd=cbind(runif(n.rnd,min(beta0.vals),max(beta0.vals)),runif(n.rnd,min(beta1.vals),max(beta1.vals)))
beta.grd=rbind(beta.grd,beta.rnd.grd)
  
N = 22670
count = 0

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];
    phi01 = exp(logit)/(1 + exp(logit))
      
    if(is.na(phi01)){ 
      count = count + 1
    }
  }
}

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

这里，$\beta_0 \in (-0.3, 0.1), \beta_1 \in (-7, 0), N = 22670, N_\text{sum} \in (1, N)$。注意 $N$ 和 $N_\text{sum}$ 是整数，而 beta 值可能不是。

由于在数学上，$\phi_{01} \in (0,1)$，我假设 NA 的出现是因为 R 不喜欢极小价值观。我也收到了大量的 NA 值。比数字更重要。为什么我会在这段代码中得到 NA？

Answer 1

在count = count + 1 旁边包含print(logit)，您会发现很多logit > 1000 值。 exp(1000) == Inf 所以你将Inf 除以Inf 得到NaN 和NaN 是NA：

> exp(500)
[1] 1.403592e+217
> Inf/Inf
[1] NaN
> is.na(NaN)
[1] TRUE

因此，您的问题不是太小，而是从exp(x) 的评估中首先出现的大量问题x 大于大约 700：

> exp(709)
[1] 8.218407e+307
> exp(710)
[1] Inf

Answer 2

Bernhard's answer 正确识别问题： 如果logit 很大，则exp(logit) = Inf。 这是一个解决方案：

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

我们可以使用更稳定的1 / (1 + exp(-logit) （为了让自己相信这一点，请将您的表达式乘以exp(-logit) / exp(-logit)）， 幸运的是，无论哪种方式，R 都有一个内置函数plogis()，可以快速准确地计算这些概率。 您可以从帮助文件 (?plogis) 中看到，此函数会计算我给出的表达式，但您也可以仔细检查以确保自己

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