Hồi quy logistic bằng cách sử dụng SciPy

Question

Tôi đang cố gắng mã hóa hồi quy logistic bằng Python bằng hàm SciPy fmin_bfgs, nhưng đang chạy vào một số vấn đề. Tôi đã viết các hàm cho hàm chuyển đổi logistic (sigmoid), và hàm chi phí, và các công việc đó tốt (tôi đã sử dụng các giá trị được tối ưu hóa của vector tham số được tìm thấy thông qua phần mềm đóng hộp để kiểm tra các hàm, và các phần tử phù hợp). Tôi không chắc chắn về việc triển khai chức năng gradient của tôi, nhưng có vẻ hợp lý.

Đây là mã:

# purpose: logistic regression 
import numpy as np 
import scipy.optimize 

# prepare the data 
data = np.loadtxt('data.csv', delimiter=',', skiprows=1) 
vY = data[:, 0] 
mX = data[:, 1:] 
intercept = np.ones(mX.shape[0]).reshape(mX.shape[0], 1) 
mX = np.concatenate((intercept, mX), axis = 1) 
iK = mX.shape[1] 
iN = mX.shape[0] 

# logistic transformation 
def logit(mX, vBeta): 
    return((1/(1.0 + np.exp(-np.dot(mX, vBeta))))) 

# test function call 
vBeta0 = np.array([-.10296645, -.0332327, -.01209484, .44626211, .92554137, .53973828, 
    1.7993371, .7148045 ]) 
logit(mX, vBeta0) 

# cost function 
def logLikelihoodLogit(vBeta, mX, vY): 
    return(-(np.sum(vY*np.log(logit(mX, vBeta)) + (1-vY)*(np.log(1-logit(mX, vBeta)))))) 
logLikelihoodLogit(vBeta0, mX, vY) # test function call 

# gradient function 
def likelihoodScore(vBeta, mX, vY): 
    return(np.dot(mX.T, 
        ((np.dot(mX, vBeta) - vY)/ 
        np.dot(mX, vBeta)).reshape(iN, 1)).reshape(iK, 1)) 

likelihoodScore(vBeta0, mX, vY).shape # test function call 

# optimize the function (without gradient) 
optimLogit = scipy.optimize.fmin_bfgs(logLikelihoodLogit, 
            x0 = np.array([-.1, -.03, -.01, .44, .92, .53, 
              1.8, .71]), 
            args = (mX, vY), gtol = 1e-3) 

# optimize the function (with gradient) 
optimLogit = scipy.optimize.fmin_bfgs(logLikelihoodLogit, 
            x0 = np.array([-.1, -.03, -.01, .44, .92, .53, 
              1.8, .71]), fprime = likelihoodScore, 
            args = (mX, vY), gtol = 1e-3) 

• Việc tối ưu hóa đầu tiên (không có gradient) kết thúc với một tổng thể rất nhiều thứ về phép chia cho không.

• Tối ưu hóa thứ hai (có độ dốc) kết thúc bằng ma trận không phù hợp với lỗi, điều này có thể có nghĩa là tôi đã nhận được sai số trả về sai.

Bất kỳ trợ giúp nào được đánh giá cao. Nếu bất cứ ai muốn thử điều này, dữ liệu được bao gồm dưới đây.

low,age,lwt,race,smoke,ptl,ht,ui 
0,19,182,2,0,0,0,1 
0,33,155,3,0,0,0,0 
0,20,105,1,1,0,0,0 
0,21,108,1,1,0,0,1 
0,18,107,1,1,0,0,1 
0,21,124,3,0,0,0,0 
0,22,118,1,0,0,0,0 
0,17,103,3,0,0,0,0 
0,29,123,1,1,0,0,0 
0,26,113,1,1,0,0,0 
0,19,95,3,0,0,0,0 
0,19,150,3,0,0,0,0 
0,22,95,3,0,0,1,0 
0,30,107,3,0,1,0,1 
0,18,100,1,1,0,0,0 
0,18,100,1,1,0,0,0 
0,15,98,2,0,0,0,0 
0,25,118,1,1,0,0,0 
0,20,120,3,0,0,0,1 
0,28,120,1,1,0,0,0 
0,32,121,3,0,0,0,0 
0,31,100,1,0,0,0,1 
0,36,202,1,0,0,0,0 
0,28,120,3,0,0,0,0 
0,25,120,3,0,0,0,1 
0,28,167,1,0,0,0,0 
0,17,122,1,1,0,0,0 
0,29,150,1,0,0,0,0 
0,26,168,2,1,0,0,0 
0,17,113,2,0,0,0,0 
0,17,113,2,0,0,0,0 
0,24,90,1,1,1,0,0 
0,35,121,2,1,1,0,0 
0,25,155,1,0,0,0,0 
0,25,125,2,0,0,0,0 
0,29,140,1,1,0,0,0 
0,19,138,1,1,0,0,0 
0,27,124,1,1,0,0,0 
0,31,215,1,1,0,0,0 
0,33,109,1,1,0,0,0 
0,21,185,2,1,0,0,0 
0,19,189,1,0,0,0,0 
0,23,130,2,0,0,0,0 
0,21,160,1,0,0,0,0 
0,18,90,1,1,0,0,1 
0,18,90,1,1,0,0,1 
0,32,132,1,0,0,0,0 
0,19,132,3,0,0,0,0 
0,24,115,1,0,0,0,0 
0,22,85,3,1,0,0,0 
0,22,120,1,0,0,1,0 
0,23,128,3,0,0,0,0 
0,22,130,1,1,0,0,0 
0,30,95,1,1,0,0,0 
0,19,115,3,0,0,0,0 
0,16,110,3,0,0,0,0 
0,21,110,3,1,0,0,1 
0,30,153,3,0,0,0,0 
0,20,103,3,0,0,0,0 
0,17,119,3,0,0,0,0 
0,17,119,3,0,0,0,0 
0,23,119,3,0,0,0,0 
0,24,110,3,0,0,0,0 
0,28,140,1,0,0,0,0 
0,26,133,3,1,2,0,0 
0,20,169,3,0,1,0,1 
0,24,115,3,0,0,0,0 
0,28,250,3,1,0,0,0 
0,20,141,1,0,2,0,1 
0,22,158,2,0,1,0,0 
0,22,112,1,1,2,0,0 
0,31,150,3,1,0,0,0 
0,23,115,3,1,0,0,0 
0,16,112,2,0,0,0,0 
0,16,135,1,1,0,0,0 
0,18,229,2,0,0,0,0 
0,25,140,1,0,0,0,0 
0,32,134,1,1,1,0,0 
0,20,121,2,1,0,0,0 
0,23,190,1,0,0,0,0 
0,22,131,1,0,0,0,0 
0,32,170,1,0,0,0,0 
0,30,110,3,0,0,0,0 
0,20,127,3,0,0,0,0 
0,23,123,3,0,0,0,0 
0,17,120,3,1,0,0,0 
0,19,105,3,0,0,0,0 
0,23,130,1,0,0,0,0 
0,36,175,1,0,0,0,0 
0,22,125,1,0,0,0,0 
0,24,133,1,0,0,0,0 
0,21,134,3,0,0,0,0 
0,19,235,1,1,0,1,0 
0,25,95,1,1,3,0,1 
0,16,135,1,1,0,0,0 
0,29,135,1,0,0,0,0 
0,29,154,1,0,0,0,0 
0,19,147,1,1,0,0,0 
0,19,147,1,1,0,0,0 
0,30,137,1,0,0,0,0 
0,24,110,1,0,0,0,0 
0,19,184,1,1,0,1,0 
0,24,110,3,0,1,0,0 
0,23,110,1,0,0,0,0 
0,20,120,3,0,0,0,0 
0,25,241,2,0,0,1,0 
0,30,112,1,0,0,0,0 
0,22,169,1,0,0,0,0 
0,18,120,1,1,0,0,0 
0,16,170,2,0,0,0,0 
0,32,186,1,0,0,0,0 
0,18,120,3,0,0,0,0 
0,29,130,1,1,0,0,0 
0,33,117,1,0,0,0,1 
0,20,170,1,1,0,0,0 
0,28,134,3,0,0,0,0 
0,14,135,1,0,0,0,0 
0,28,130,3,0,0,0,0 
0,25,120,1,0,0,0,0 
0,16,95,3,0,0,0,0 
0,20,158,1,0,0,0,0 
0,26,160,3,0,0,0,0 
0,21,115,1,0,0,0,0 
0,22,129,1,0,0,0,0 
0,25,130,1,0,0,0,0 
0,31,120,1,0,0,0,0 
0,35,170,1,0,1,0,0 
0,19,120,1,1,0,0,0 
0,24,116,1,0,0,0,0 
0,45,123,1,0,0,0,0 
1,28,120,3,1,1,0,1 
1,29,130,1,0,0,0,1 
1,34,187,2,1,0,1,0 
1,25,105,3,0,1,1,0 
1,25,85,3,0,0,0,1 
1,27,150,3,0,0,0,0 
1,23,97,3,0,0,0,1 
1,24,128,2,0,1,0,0 
1,24,132,3,0,0,1,0 
1,21,165,1,1,0,1,0 
1,32,105,1,1,0,0,0 
1,19,91,1,1,2,0,1 
1,25,115,3,0,0,0,0 
1,16,130,3,0,0,0,0 
1,25,92,1,1,0,0,0 
1,20,150,1,1,0,0,0 
1,21,200,2,0,0,0,1 
1,24,155,1,1,1,0,0 
1,21,103,3,0,0,0,0 
1,20,125,3,0,0,0,1 
1,25,89,3,0,2,0,0 
1,19,102,1,0,0,0,0 
1,19,112,1,1,0,0,1 
1,26,117,1,1,1,0,0 
1,24,138,1,0,0,0,0 
1,17,130,3,1,1,0,1 
1,20,120,2,1,0,0,0 
1,22,130,1,1,1,0,1 
1,27,130,2,0,0,0,1 
1,20,80,3,1,0,0,1 
1,17,110,1,1,0,0,0 
1,25,105,3,0,1,0,0 
1,20,109,3,0,0,0,0 
1,18,148,3,0,0,0,0 
1,18,110,2,1,1,0,0 
1,20,121,1,1,1,0,1 
1,21,100,3,0,1,0,0 
1,26,96,3,0,0,0,0 
1,31,102,1,1,1,0,0 
1,15,110,1,0,0,0,0 
1,23,187,2,1,0,0,0 
1,20,122,2,1,0,0,0 
1,24,105,2,1,0,0,0 
1,15,115,3,0,0,0,1 
1,23,120,3,0,0,0,0 
1,30,142,1,1,1,0,0 
1,22,130,1,1,0,0,0 
1,17,120,1,1,0,0,0 
1,23,110,1,1,1,0,0 
1,17,120,2,0,0,0,0 
1,26,154,3,0,1,1,0 
1,20,106,3,0,0,0,0 
1,26,190,1,1,0,0,0 
1,14,101,3,1,1,0,0 
1,28,95,1,1,0,0,0 
1,14,100,3,0,0,0,0 
1,23,94,3,1,0,0,0 
1,17,142,2,0,0,1,0 
1,21,130,1,1,0,1,0 

Answer 1

Đây là the answer I sent back to the SciPy list nơi câu hỏi này được đăng chéo. Cảm ơn @tiago vì câu trả lời của anh ấy. Về cơ bản, tôi reparametrized chức năng khả năng. Ngoài ra, thêm một cuộc gọi đến chức năng check_grad.

#===================================================== 
# purpose: logistic regression 
import numpy as np 
import scipy as sp 
import scipy.optimize 

import matplotlib as mpl 
import os 

# prepare the data 
data = np.loadtxt('data.csv', delimiter=',', skiprows=1) 
vY = data[:, 0] 
mX = data[:, 1:] 
# mX = (mX - np.mean(mX))/np.std(mX) # standardize the data; if required 

intercept = np.ones(mX.shape[0]).reshape(mX.shape[0], 1) 
mX = np.concatenate((intercept, mX), axis = 1) 
iK = mX.shape[1] 
iN = mX.shape[0] 

# logistic transformation 
def logit(mX, vBeta): 
    return((np.exp(np.dot(mX, vBeta))/(1.0 + np.exp(np.dot(mX, vBeta))))) 

# test function call 
vBeta0 = np.array([-.10296645, -.0332327, -.01209484, .44626211, .92554137, .53973828, 
    1.7993371, .7148045 ]) 
logit(mX, vBeta0) 

# cost function 
def logLikelihoodLogit(vBeta, mX, vY): 
    return(-(np.sum(vY*np.log(logit(mX, vBeta)) + (1-vY)*(np.log(1-logit(mX, vBeta)))))) 
logLikelihoodLogit(vBeta0, mX, vY) # test function call 

# different parametrization of the cost function 
def logLikelihoodLogitVerbose(vBeta, mX, vY): 
    return(-(np.sum(vY*(np.dot(mX, vBeta) - np.log((1.0 + np.exp(np.dot(mX, vBeta))))) + 
        (1-vY)*(-np.log((1.0 + np.exp(np.dot(mX, vBeta)))))))) 
logLikelihoodLogitVerbose(vBeta0, mX, vY) # test function call 

# gradient function 
def likelihoodScore(vBeta, mX, vY): 
    return(np.dot(mX.T, 
        (logit(mX, vBeta) - vY))) 
likelihoodScore(vBeta0, mX, vY).shape # test function call 
sp.optimize.check_grad(logLikelihoodLogitVerbose, likelihoodScore, 
         vBeta0, mX, vY) # check that the analytical gradient is close to 
               # numerical gradient 

# optimize the function (without gradient) 
optimLogit = scipy.optimize.fmin_bfgs(logLikelihoodLogitVerbose, 
            x0 = np.array([-.1, -.03, -.01, .44, .92, .53, 
              1.8, .71]), 
            args = (mX, vY), gtol = 1e-3) 

# optimize the function (with gradient) 
optimLogit = scipy.optimize.fmin_bfgs(logLikelihoodLogitVerbose, 
            x0 = np.array([-.1, -.03, -.01, .44, .92, .53, 
              1.8, .71]), fprime = likelihoodScore, 
            args = (mX, vY), gtol = 1e-3) 
#===================================================== 

Answer 2

Vấn đề của bạn là chức năng bạn đang cố gắng để giảm thiểu, logLikelihoodLogit, sẽ trở lại NaN với giá trị rất gần với ước tính ban đầu của bạn. Và nó cũng sẽ cố gắng để đánh giá logarit tiêu cực và gặp phải các vấn đề khác. fmin_bfgs không biết về điều này, sẽ cố gắng đánh giá chức năng cho các giá trị như vậy và gặp rắc rối.

Tôi khuyên bạn nên sử dụng tối ưu hóa được giới hạn thay thế. Bạn có thể sử dụng số optimize.fmin_l_bfgs_b của scipy cho việc này. Nó sử dụng một thuật toán tương tự như fmin_bfgs, nhưng nó hỗ trợ giới hạn trong không gian tham số. Bạn gọi nó tương tự, chỉ cần thêm một từ khóa giới hạn. Dưới đây là một ví dụ đơn giản về cách bạn muốn gọi fmin_l_bfgs_b:

from scipy.optimize import fmin_bfgs, fmin_l_bfgs_b 

# list of bounds: each item is a tuple with the (lower, upper) bounds 
bd = [(0, 1.), ...] 

test = fmin_l_bfgs_b(logLikelihoodLogit, x0=x0, args=(mX, vY), bounds=bd, 
         approx_grad=True) 

Ở đây tôi đang sử dụng một gradient gần đúng (dường như làm việc tốt với dữ liệu của bạn), nhưng bạn có thể vượt qua fprime như trong ví dụ của bạn (I don' t có thời gian để kiểm tra tính chính xác của nó). Bạn sẽ biết không gian tham số của mình tốt hơn tôi, chỉ cần đảm bảo xây dựng mảng giới hạn cho tất cả các giá trị có ý nghĩa mà các tham số của bạn có thể thực hiện.

Answer 3

Tôi đang đối mặt với cùng một vấn đề. Khi tôi thử nghiệm với các thuật toán thực hiện khác nhau trong scipy.optimize.minimize, tôi thấy rằng để tìm các tham số hồi quy logistic tối ưu cho tập dữ liệu của tôi, Newton Conjugate Gradient tỏ ra hữu ích. Cuộc gọi có thể được thực hiện theo cách như sau:

Result = scipy.optimize.minimize(fun = logLikelihoodLogit, 
           x0 = np.array([-.1, -.03, -.01, .44, .92, .53,1.8, .71]), 
           args = (mX, vY), 
           method = 'TNC', 
           jac = likelihoodScore); 
optimLogit = Result.x;