2

这是我关于堆栈溢出的第一篇文章，所以请放纵！

我正在尝试为多项式（或条件）logit 模型开发代码（如http://data.princeton.edu/wws509/notes/c6s3.html中所述）这种类型的模型可以看作是经典的二元逻辑回归。多项式 logit 模型经常用于离散选择建模。

在下面的代码中，我描述了一个模型，解释了当所有参与者 (N) 面临相同的选择题 (T) 时在 2 个选项 (J) 之间的选择。我已经为 R 开发了一段类似的代码，它工作得很好，但速度很慢（因为“优化”例程）。由于溢出问题，使用 PYTHON 开发这种类型的代码似乎更加困难。第一个操作（即 numpsy.dot(X,beta)）可能导致非常大/小的值，然后变得与 numpsy.exp()“不兼容”。如https://lingpipe-blog.com/2012/02/16/howprevent-overflow-underflow-logistic-regression/中所示, 溢出问题往往发生在：abs(np.dot(X,beta)) > 500 时。初始溢出问题具有重要的后果，因为它们将导致其余代码中的 0、NaN 和 Inf 值，从而生成其他问题（例如，np.log(0) 或 number/Inf）。我花了很多时间试图用十进制模块、bigfloat 模块等解决这个问题，但似乎没有用（使用 numpy 处理 exp 中的溢出）。

第二个最佳选项似乎通过包含一些舍入操作（例如 num[num > 400] = 400）来修改似然函数，然后使用“Nelder-Mead”算法来最小化对数似然。理想情况下，我会使用“BFGS”算法，但它与舍入操作不兼容。

关于如何修复代码中的溢出问题的任何想法？

备注：非常大的 np.dot(X,beta) 值在 R、MATLAB 和 STATA 中似乎不是问题。

非常感谢您的帮助！

蟒蛇代码

N = len(set(data[:,0]))
T = 20
J = 2
Y = np.reshape(data[:,12],(N*T,J))
X = np.matrix(data[:,[33,10,5,6,7,8,9]])

def fn_mnl(beta):
    num = np.dot(X,beta)
    num[num >  400] =  400 # to avoid overflow  with np.exp()
    num[num < -400] = -400 # to avoid underflow with np.exp()
    num = np.exp(num)
    num = num.reshape(N*T,J)
    den = np.sum(num, axis=1)
    prbj = num/den
    prbi = prbj[Y==1]
    prbi[prbi <= 0] = 0.0000001 # to avoid issue with np.log()
    llik = np.sum(np.log(prbi))
    return(-llik)

sv = [0]*X.shape[1]
res = sc.minimize(fn_mnl, sv, method='Nelder-Mead')

代码

LOGIT = function(Beta, DATA, X, Y){
  num  = exp(as.matrix(DATA[,X]) %*% as.vector(Beta))
  mat1 = matrix(num, ncol=2, byrow=T)
  prbj = mat1/rowSums(mat1)
  prbi = prbj[Y==1]
  logl = sum(log(prbi))
  return(-logl)}
4

0 回答 0