将概率模型拟合到 3 个标准正态协变量的 200 个观测值

Question

背景和任务：考虑一个大小为 n 的随机样本，其结果为二元结果 Y_i。假设 Y_i ~ Bern(pi_i)。假设 probit link 函数 pi_i=Phi(X_i^T beta).

创建一个 X 矩阵，其中包含 200 个关于三个协变量的观测值，每个协变量都具有标准正态分布。假设 probit 模型对数据是正确的，生成一个 y 向量，即选择一些 beta 值并将这些值与 X 值和 probit link 函数一起使用以生成一组结果 y.

在 R 中编写一个函数来拟合 probit 模型，采用矢量响应 y 和协变量矢量以及截距 X。运行基于这些数据的模型并将系数估计值与真实值进行比较.

代码:

X=rnorm(200*3) # generate 200x3=600 random standard normal values
dim(X)=c(200, 3) # set the dimensions of X to be 200x3
X=cbind(1, X) # add a column of 1's for the intercept
X # print X
beta=c(1,4,2,3) # choose some values of beta
pi_i=pnorm(X%*%beta)
for (i in 1:200) { # generate y vector
  y[i]=rbinom(1, 1, pi_i[i])
}

loglik=function(par, X, y) {
  pi_est=pnorm(X%*%par) # probit link function
  ll=sum(y*log(pi_est)+(1-y)*log(1-pi_est)) # log likelihood for bernoulli sample
  return(ll)
}
opt.out=optim(par=c(0,0,0,0), fn=loglik, X=X, y=y, method="BFGS", control=list(fnscale=-1), hessian=TRUE) # error in this line

问题：我遇到了错误

Error in optim(par = c(0, 0, 0, 0), fn = loglik, X = X, y = y, method = "BFGS", : non-finite finite-difference value [3]

有人知道这是为什么吗？

Answer 1

当运行问题中的loglik 函数可以return NaN 值时使用下面的数据。这是由于 pi_est 在数值上接近 1，因此术语 log(1 - pi_est) 等同于 log(0) 导致无限值。

par <- c(1, 4, 2, 3)
pi_est <- pnorm(X %*% par) 
ll <- sum(y* log(pi_est) + (1 - y)* log(1 - pi_est)) 
ll
# [1] NaN\

特别是 pi_est 的值被评估为 1 -- 数值准确性问题。

idx <- which(is.infinite(log(1 - pi_est)))
print(pi_est[idx], digits=21)
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

您可以通过计算 pnorm 内的积分 log 来降低这种情况发生的可能性。还要注意正常的 CDF(-x) = 1 - CDF(x)。将 log(pi_est) 替换为 pnorm(X %*% par, log.p=TRUE) 并将 log(1 - pi_est) 替换为 pnorm(X %*% par, log.p=TRUE, lower.tail = FALSE)（等于 pnorm(-X %*% par, log.p=TRUE)）会导致更稳定的计算。

loglik <- function(par, X, y) {
   lp = X %*% par
   pi_est1 = pnorm(lp, log.p=TRUE) 
   pi_est2 = pnorm(lp, log.p=TRUE, lower.tail=FALSE) 
   ll = -sum(y*pi_est1 + (1-y)* pi_est2)
   return(ll)
 }
opt.out <- optim(par=c(1,1,1,1), 
                 fn=loglik, X=X, y=y, 
                 method="BFGS",
                 hessian=TRUE) 
opt.out$par
# [1] 1.207355 4.585248 2.064004 3.430316

# Using `glm`
m = glm(y ~ X-1, family=binomial(link="probit"))
# Warning message:
# glm.fit: fitted probabilities numerically 0 or 1 occurred 
coef(m)
#       X1       X2       X3       X4 
# 1.207346 4.585221 2.063990 3.430295

可能有一种方法可以避免计算两次积分。

数据

set.seed(65819138)
X <- matrix(rnorm(200*3), ncol=3)
X <- cbind(1, X) 
beta <- c(1,4,2,3) 
pi_i <- pnorm(X%*%beta)
y <- rbinom(200, 1, pi_i)

将概率模型拟合到 3 个标准正态协变量的 200 个观测值

Fitting a probit model to 200 observations of 3 standard normal covariates

optimization

r

glm