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

如何解决为什么我在 R 的这个计算中得到 NAs？

在处理 Rcpp 程序时，我使用了 sample() 函数，这给了我以下错误：“概率不允许 NA。”我将这个问题追溯到我使用的概率向量中包含 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,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 值。比数字更重要。为什么我会在这段代码中得到 NAs？

解决方法

在 print(logit) 旁边包含 count = count + 1，您会发现很多 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

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