Motivation
(Note) "이 내용은 Week 4: Unconstrained Optimization 1 의 응용에 대한 것으로 이차근사를 이용한 최적화문제에 대한 기초적인 지식이 필요합니다."
여기서는 부호제약 조건을 가지는 로지스틱 회귀모형을 적합하는 방법에 대해서 알아봅니다.
"3년동안" 확보한 재무제표 정보를 이용해서 어떤 기업의 신용을 평가하는 문제를 생각해봅시다. 재무제표 항복별 부도율을 높이는 것과 낮추는 것이 있을 것입니다. 사전적으로 변수가 가져야할 부호에 대해서 제약조건이 있는 경우가 있습니다. 이러한 제약조건이 있다면 로지스틱 회귀모형에 얻어진 회귀계수가 해당 부호제약조건을 만족하길 원할 것입니다.
한편 같은 항목에 대해서 최근 재무제표정보가 우선적으로 반영되길 바랄 수도 있습니다. 예를 들어 자본 대 부채비율의 크기는 부도율을 높인다고 알려져 있고, 최근 부채비율이 부도율에 더 큰 영향을 줘야한다고 믿는 것은 자연스럽습니다.
위에서 언급한 두 가지 사항을 자본 대 부채비율에 대응되는 회귀계수를 모형화하는데 수학적인 표현을 빌어 쓰면 다음과 같습니다.
- 부채비율의 회귀계수는 모두 0이하다.
- 전기의 회귀계수가 당기의 회귀계수보다 작을 수 없다. (음의 효과가 클 수 없다.)
- 전전기의 회귀계수가 전기 회귀계수보다 작을 수 없다. (음의 효과가 클 수 없다.)
우리는 이 제약조건을 반영하는 로지스틱 회귀모형을 적합해 보겠습니다.
모형 만들기
먼저 시뮬레이션 샘플을 만들어보겠습니다. $n$ 은 샘플 수, $p$ 는 당기에 사용할 변수의 수 입니다. 이 로지스틱 모형은 당기를 포함한 3년간 변수를 사용하므로 변수의 수는 $3p$ 개 입니다.
# 데이터의 생성
import numpy as np
# 샘플 수
n = 500
# 변수 수 (당기)
p = 10
# 로지스틱 모형 생성 (변수수 p*3: 당기, 전기, 전전기)
np.random.seed(1)
X = np.random.normal(size = (n,3*p))
true_beta = np.random.normal(size = 3*p)
Xb = X @ true_beta
phat = 1/(1+np.exp(-Xb))
y = np.random.binomial(1,phat)
우리가 만든 데이터의 형식은 다음과 같음을 가정합니다.
- 시점별 p개의 변수를 고려하고 있습니다.
- (당기, 전기, 전전기) 3년치의 데이터를 가지고 있습니다.
- 데이터는 당기 (p개), 전기 (p개), 전전기(p개) 의 순으로 정렬되야 합니다.
- 기간별로 변수의 순서는 같아야 합니다.
다시 말해 만약 실제 데이터를 사용하는 경우 위와 같이 변수들이 정렬되어 있어야 함을 의미합니다.
다음으로 변수의 부호를 정합니다. 이 부호는 사전시직에 의한 것으로 신중하게 결정되어야 합니다. 당기에 주어진 변수 중에서 양의 부호를 가지는 것의 index는 g1 음의 부호를 가지는 index 는 g2 에 리스트 형태로 입력합니다. 내부적으로 부호가 정해지지 않는 변수의 index 는 g3 로 계산됩니다.
# Group 의 지정
g1 = [0,2,4]
g2 = [6,7,9]
# g3 의 계산
g = set(range(0,p))
g3 = list((g - set(g1)) - set(g2))
g3.sort()
g1, g2, g3 에 있는 변수 이름을 편의상 변수1*,.., 변수p* 라 표시하겠습니다. 이어서 기간별 부등제약 조건을 쉽게 주기 위해 함수 내부에서 입력변수 X를 재정렬하고 마지막열은 상수항을 처리하기 위해 1 열벡터를 추가합니다. 정렬 후에 정렬된 입력변수의 순서는 다음과 같습니다
(변수1* 당기, 변수1* 전기, 변수1* 전전기, 변수2* 당기,..., 변수p* 전전기, 상수항)
g = g1 + g2 + g3
v = []
for i in range(p):
v.append(g[i])
v.append(g[i]+p)
v.append(g[i]+2*p)
print("v:", v)
print("g:", g)
x = X[:,v]
x = np.concatenate((x,np.ones((n,1))),axis = 1)
우리는 재정렬된 $i$번째 공변량을 다음과 같이 표시하겠습니다. $$x= (x_{1,1}, x_{1,2},x_{1,1}, x_{2,1}, \cdots, x_{p,3}, 1)\in \mathbb{R}^{3p+1}$$
다음은 제약조건을 줄 수 있는 행렬을 구성합니다. G1은 isotonic 양수제약 조건을 가지는 변수에 대한 부등식을 만드는 행렬입니다. $$\begin{eqnarray} \beta_{\cdot,2} &\leq& \beta_{\cdot,1} \\ \beta_{\cdot,3} &\leq & \beta_{\cdot,2}\\ \beta_{\cdot,3} &\leq &0 \end{eqnarray}$$ 이를 행렬과 벡터로 표시하면 다음과 같습니다. $$ G_1 = \left(\begin{array}{ccc} -1 & 1 & 0 \\ 0 & -1 & 1 \\ 0 & 0 & -1 \end{array}\right), ~~\beta = \left(\begin{array}{c} \beta_{\cdot,1} \\ \beta_{\cdot,2} \\ \beta_{\cdot,3} \end{array}\right), ~~ h = 0 \in \mathbb{R}^3 $$
$G_2 = -G_1$로 음수제약 조건을 정의하는 행렬입니다. 위쪽 Motivation에서 이미 설명하였습니다.
# block matrix 만들기
G1 = np.array([[-1,1,0], [0,-1,1],[0,0,-1],
[-1,0,0], [0,-1,0], [0,0,-1]], dtype=float)
G2 = -G1.copy()
상수항을 포함하여 제약조건이 없는 변수를 처리하기 위해 $G_1$, $G_2$를 이용해 block matrix 를 만든 후에 영행렬을 추가하여 행렬 $G$를 만듭니다.
# 블록 행렬 함수
matrices = []
for i in g1:
matrices.append(G1)
for i in g2:
matrices.append(G2)
def block_diag(*arrs):
# Calculate the total size of the resulting matrix
rows = sum(arr.shape[0] for arr in arrs)
cols = sum(arr.shape[1] for arr in arrs)
# Create a zero matrix with the calculated size
result = np.zeros((rows, cols))
# Fill in the block diagonal
r, c = 0, 0
for arr in arrs:
rows, cols = arr.shape
result[r:r+rows, c:c+cols] = arr
r += rows
c += cols
return result
# Create the block diagonal matrix
G = block_diag(*matrices)
G.shape
G = np.concatenate((G, np.zeros((G.shape[0],3*p-G.shape[1] + 1) )),
axis = 1)
다음으로 로지스틱 회귀모형에 대한 미분, 헤시안을 정의합니다. 여기서는 음로그우도함수를 최소화하는 목적함수로 설정하였습니다. 로지스틱 회귀모형의 미분과 헤시안은 Week 04: Unconstrained Optimization (slide 1) 의 Example 3 을 확인하세요
# 초기값 설정
init_b = np.zeros(G.shape[1])
# 시그모이드 함수 정의
def sigmoid(z):
return 1 / (1 + np.exp(-z))
# 도함수 정의
def neg_log_likelihood_gradient(X, y, beta):
m = X.shape[0]
predictions = sigmoid(np.dot(X, beta))
gradient = np.dot(X.T, (predictions - y)) / m
return gradient
# 헤시안 함수 정의
def neg_log_likelihood_hessian(X, beta):
m = X.shape[0]
predictions = sigmoid(np.dot(X, beta))
D = np.diag(predictions * (1 - predictions))
hessian = np.dot(np.dot(X.T, D), X) / m
return hessian
# 미분 함수 계산
gradient = neg_log_likelihood_gradient(x, y, init_b)
# 헤시안 행렬 계산
hessian = neg_log_likelihood_hessian(x, init_b)
다음으로 QP 식에 사용할 수 있는 형태로 P, q, G, h, 를 설정합니다. QP 함수는 다음과 같은 이차함수의 최적화 문제를 풀 수 있습니다. $$ \begin{eqnarray} L(\beta) &=& \frac{1}{2}\beta^\top P \beta + q^\top \beta\\ s.t& & G\beta \leq h \end{eqnarray} $$
로지스틱 음우도함수의 테일러 전개를 이용하여 QP 의 목적함수의 형태로 만들 수 있습니다. 아래 식에서 어떤 값 $\beta_0$ 음우도함수의 테일러전개를 확인해보세요.
$$\begin{eqnarray} l(\beta) &=& l(\beta_0) + \nabla l(\beta_0)^\top (\beta - \beta_0) + \frac{1}{2}(\beta - \beta_0)^\top \nabla^2 l(\beta_0) (\beta - \beta_0) + const.\\ &=&\frac{1}{2} \beta^\top \nabla^2 l(\beta_0) \beta + (\nabla l(\beta_0) - \nabla^2 l(\beta_0)\beta_0)^\top \beta + const.. \end{eqnarray} $$ 여기서 $\nabla l(\beta_0)$이 $\beta_0$ 에서 gradient, $\nabla^2 l(\beta_0)$ 이 Hessian 행렬입니다.
# QP form
P = hessian
q = gradient - hessian @init_b
h = np.zeros(G.shape[0])
cvxopt 에서 정의한 matrix 형태로 변환합니다.
# 형식 변환
from cvxopt import matrix, solvers
# 메시지 끄기
solvers.options['show_progress'] = False
cP = matrix(P)
cq = matrix(q)
cG = matrix(G)
ch = matrix(h)
QP 문제플 풀고 해를 얻습니다.
# solution
sol = solvers.qp(cP,cq,cG,ch)
sol['primal objective']
위 과정은 목적함수의 이사근차 후 해를 한번 구한 것이기 때문에 이 과정을 반복합니다.
# 미분 함수 계산
b = init_b.copy()
b = b.squeeze()
for i in range(500):
gradient = neg_log_likelihood_gradient(x, y, b)
# 헤시안 행렬 계산
hessian = neg_log_likelihood_hessian(x, b)
P = hessian
q = gradient - hessian @init_b
h = np.zeros(G.shape[0])
cP = matrix(P)
cq = matrix(q)
cG = matrix(G)
ch = matrix(h)
sol = solvers.qp(cP,cq,cG,ch)
b = np.array(sol['x']).squeeze()
if i%10==9:
print("primal objective:", round(sol['primal objective'],5))
마지막으로 원래 입력변수의 위치로 돌려줍니다. 원본 데이터를 재정렬할 때 사용한 변수 v 를 활용합니다.
beta_est = np.zeros(3*p)
intercept = b[-1]
beta_est
# v 는 변수의 위치 정보가지고 있음
beta_est[v] = b[0:-1]