- 공유 링크 만들기
- 이메일
- 기타 앱
내용
로지스틱회귀(Logistic Regression)
로지스틱회귀
로지스틱 회귀분석은 독립변수(특성, 설명변수)에 대해 반응변수(라벨)를 로짓변수(logit variavble, 반응변수의 발생 확률의 자연로그)로 변환한 후 인스턴스(smaple)가 두 개의 클래스로 구분된 라벨(반응변수) 중에 특정한 클래스에 속하는 확률(최대우도)을 추정하기 위해 사용됩니다. 그러므로 이 회귀 모델은 binary classificatier를 구축하게 됩니다.
선형회귀와 유사하게 로지스틱회귀 역시 입력 변수들과 가중치들의 곱이 계산됩니다. 그러나 선형회귀의 경우 이 곱의 결과가 직접적으로 사용되는 것에 비해 로지스틱 회귀의 경우는 시그모이드 함수(sigmoid function, σ(·))을 통해 [0, 1] 사이의 값으로 변환시킵니다. 즉, 식 1과 같은 연산에 의해 인스턴스의 결과(확률)가 계산됩니다.
$$\begin{align}\tag{1} &\hat{p}=h_\beta(x)=\sigma\left(x^T\beta\right)\\ &t=x^T\beta\\ \tag{2}&\sigma(t)=\frac{1}{1+\exp(-t)}\end{align}$$시그모이드 함수(식 2)는 변수와 가중치의 곱의 결과를 0과 1사이의 값으로 변환합니다.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import FinanceDataReader as fdr
from mpl_toolkits.mplot3d import Axes3D
font1={'size':11, 'weight':'bold'}
def sigmoid(x):
return 1/(1+np.exp(-x))
t=np.linspace(-10, 10, 100)
y=sigmoid(t)
plt.figure(figsize=(7, 3), dpi=100)
plt.plot(t, y, label="sigmoid")
plt.legend(loc='best', prop=font1)
plt.axhline(y=0.5, linestyle="--", color='k')
plt.axvline(x=0, linestyle="--", color='k')
plt.xlabel('t', fontdict=font1)
plt.ylabel('y', fontdict=font1)
plt.show()
위 결과와 같이 t<0인 경우 σ<0.5, t≥0인 경우 σ≥0.5가 됩니다. 두 클래스 0, 1을 각각 negative, positive라고한다면 식 3과 같이 나타낼 수 있습니다.
$$\tag{3}\hat{y}=\begin{cases}0&\; \text{if}\; \hat{p}\lt 0.5\\1&\; \text{if}\; \hat{p}\geq 0.5\end{cases} $$로짓함수(logit function)는 식 4와 같이 정의됩니다.
$$\begin{align}\tag{4}\text{logit}(p)&=\ln\left(\frac{p}{1-p}\right)\\&= \ln(\text{odds}) \end{align}$$그러므로 데이터가 확률(p)이라면 식 4와 같이 odds ratio를 사용합니다. odds ratio는 두 사건 사이의 상관 강도를 정량화하는 통계량으로 로그화하면 다음 표와 같이 변수의 범위를 ∞로 확장할 수 있습니다.
| 항목 | 범위 |
|---|---|
| p | [0, 1] |
| odds | [0, ∞] |
| log(odds) | [-∞, ∞] |
결과적으로 log(odds)를 적용하면 범위는 무한대로 확장되므로 범위의 제한으로 발생하는 문제들을 방지할 수 있으며 로그화에 의해 비선형 관계를 포함할 수 있습니다.
만약 확률 0, 1에 대해 선형모델을 적용하면 다음과 같은 문제가 발생하며 이것은 로지스틱 회기모델로 보상됩니다.
- 오차항은 X 독립변수의 중간 값에서 크고 극단값에서 작은 경향이 있습니다. 이는 오차가 평균이 0이어야 하고 정규 분포를 따라야 한다는 선형 회귀 가정을 위반하는 것입니다.
- X의 끝 값에서 1보다 크고 0보다 작은 무의미한 예측을 생성합니다.
- 일반 최소 자승(OLS) 추정은 비효율적이며 표준 오차는 편향됩니다.
- X의 중간 값에서 높은 오차 분산 및 끝에서 낮은 분산이 발생됩니다.
비용함수
인스턴스 중에 확률(p)이 1 또는 0과 같이 극단적인 경우의 비용 즉, 실측치와의 차이는 다음과 같습니다.
- p → 0 ⇒ log(p)=-∞, -log(p)=∞인 경우 positive에 해당하는 인스턴스(y=1)가 0에 접근할 경우의 비용은 ∞에 접근
- p → 1 ⇒ log(1-p)=-∞, -log(1-p)=∞인 경우 negative에 해당하는 인스턴스(y=0)가 1에 접근할 경우의 비용은 ∞에 접근
위 결과는 식 5와 같이 매개변수 β에 대한 하나의 인스턴스에 대한 비용함수(c())로 나타낼 수 있습니다.
$$\tag{5}c(\beta)=\begin{cases}-\ln(\hat{p})&\; \text{if y=1}\\-\ln(1-\hat{p})&\; \text{if y=0}\end{cases}$$식 5는 단일 인스턴스에 대한 비용을 나타내므로 전체 훈련세트에 대한 비용은 모든 비용의 합에 대한 평균을 사용합니다. 식 6과 같이 하나의 식으로 표현할 수 있습니다.
$$\tag{6}J(\beta)=-\frac{1}{n}\sum^n_{i=1}\left[y^{(i)}\log\left(\hat{p}^{(i)}\right)+\left(1-y^{(i)}\right)\log\left(1-\hat{p}^{(i)}\right)\right]$$식 6의 비용함수를 최소로 하는 β를 직접적으로 찿을 수 있는 식은 없습니다. 그러나 경사하강법(gradient descent)을 사용하여 비용함수를 최소로하는 β를 발견할 수 있습니다. 식 6의 미분은 식 7과 같습니다.
$$\tag{7}\frac{\partial}{\partial \beta_j}J(\beta)=\frac{1}{n}\sum^n_{i=1}\left(\sigma\left(\beta^Tx^{(i)}\right)-y^{(i)}\right)x^{(i)}_j$$식 7은 j 번째 회귀계수(β)에 대해 훈련 데이터 셋트의 비용함수 변화량의 평균을 계산한 것입니다.
LogisticRegression()
로지스틱 회귀 모델은 sklearn.linear_model.LogisticRegression() 클래스를 사용할 수 있습니다. 이 클래스로 생성한 모델의 독립변수에 대한 예측확률은 .predict_proba() 메서드로 확인할 수 있습니다. 이 메서드는 각 인스턴스가 클래스 0에 속할 확률과 1에 속할 확률 2개를 반환합니다.
다음은 iris 데이터 셋으로 부터 petal with를 독립변수, target을 반응변수로 하고 반응변수의 클래스를 'virginica'(1)와 'not virginica'(0)으로 분류합니다.
from sklearn import datasets
iris=datasets.load_iris() iris.keys()
dict_keys(['data', 'target', 'frame', 'target_names', 'DESCR', 'feature_names', 'filename', 'data_module'])
iris['feature_names'], iris['target_names']
(['sepal length (cm)',
'sepal width (cm)',
'petal length (cm)',
'petal width (cm)'],
array(['setosa', 'versicolor', 'virginica'], dtype='<U10'))
np.unique(iris["target"], return_counts=True)
(array([0, 1, 2]), array([50, 50, 50]))
X=iris["data"][:, 3:] y=(iris["target"]==2).astype(int)
np.unique(X, return_counts=True)
(array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1. , 1.1, 1.2, 1.3, 1.4, 1.5, 1.6,
1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3, 2.4, 2.5]),
array([ 5, 29, 7, 7, 1, 1, 7, 3, 5, 13, 8, 12, 4, 2, 12, 5, 6,
6, 3, 8, 3, 3]))
iris 데이터로 부터 설정한 변수 X, y에 대해 로지스틱 모델을 생성합니다. 모델 생성을 위해 LogisticRegression()클래스를 적용하며 각 인스턴스에 대한 예측 확률을 산출해 봅니다.
from sklearn.linear_model import LogisticRegression
logit=LogisticRegression().fit(X, y)
prePro=logit.predict_proba(X) prePro[:3, :]
array([[0.99821774, 0.00178226],
[0.99821774, 0.00178226],
[0.99821774, 0.00178226]])
결정경계(decision boundary)
다음은 데이터의 petal width의 임의 값에 따라 모델로 부터 예측한 확률과 실제 데이터의 진위 여부를 나타낸 것입니다. 실제 데이터에서 petal width는 1.4이상에서 virginica임을 나타내지만 모델에 의한 분류는 약 2.0 이상에서 확신할 수 있습니다. 이 모델은 예측확률 50%를 기준으로 그 이상에서 virginica로 결정합니다. 즉, 예측확률 50%가 결정경계(decision boundary)가 됩니다.
plt.figure(dpi=100)
X1=np.linspace(0, 3, 1000).reshape(-1,1)
x1_prob=logit.predict_proba(X1)
y1=logit.predict(X1)
plt.scatter(X, y, label="data")
plt.plot(X1, x1_prob[:,0], label="not Virginica")
plt.plot(X1, x1_prob[:,1], label="Virginica")
plt.plot(X1, y1, color='red', label="model")
plt.axhline(y=0.5, linestyle="--", color='k', label="50%", alpha=0.5)
plt.legend(bbox_to_anchor=(1,1), prop=font1)
plt.xlabel("petal width (cm)", fontdict=font1)
plt.ylabel("Probability", fontdict=font1)
plt.grid(True)
plt.show()
LogisticRegression() 클래스의 decision_function(feature) 메서드는 각 샘플(인스턴스)의 신뢰도 점수를 반환합니다. 신뢰도 점수는 model에서 판단의 기준이 되는 결정경계와의 거리를 나타냅니다. 다음의 결과와 같이 결정경계에 가장 근접한 점수는 약 0.1715이며 이에 가장 근접한 두 개의 X값은 모두 1.7이지만 이에 대응하는 y값은 0과 1입니다. 그러나 이 값에 대한 예측은 모두 1입니다. 즉, 이 값 이상에서는 모두 1이 됨을 알 수 있습니다.
np.abs(logit.decision_function(X)).min()
0.17153551269203593
np.where((logit.decision_function(X) > 0.1715355)&(logit.decision_function(X) < 0.17153552))
(array([ 77, 106]),)
X[[77,106]], y[[77,106]]
(array([[1.7],
[1.7]]),
array([0, 1]))
logit.predict(X)[[77,106]]
array([1, 1])
petal length를 독립변수에 추가하여 결정경계를 확인해 봅니다.
X_f2=iris["data"][:, [2,3]] logit2=LogisticRegression().fit(X_f2, y)
fig=plt.figure(dpi=100)
x1=np.linspace(1, 7, 100)
y1=np.linspace(0, 2.5, 100)
ax=fig.add_subplot(111, projection='3d')
ax.scatter(X_f2[:,0], X_f2[:,1], y, s=5, color=['blue' if i==1 else 'red' for i in y],label="data")
ax.plot(x1, y1, logit2.predict(np.c_[x1, y1]), label="model")
ax.legend(bbox_to_anchor=(1,1), prop=font1)
plt.ylabel("petal width (cm)", fontdict=font1)
plt.xlabel("petal length (cm)", fontdict=font1)
plt.savefig('img/logit3D.png')
plt.show()
위 결과의 결정경계를 계산해 보면 다음과 같습니다. 즉, decision_function()에 의해 각 샘플의 점수를 계산하고 각 절대값의 최소(이 지점이 모델의 경계에 가장 근접한 지점)에 대응하는 X값과 그 지점의 실측 라벨과 예측 라벨을 비교해 봅니다.
decision=logit2.decision_function(X_f2) pre=logit2.predict(X_f2) scoreMin=np.abs(decision).min() scoreMin
0.07842675185569803
idxMin=np.where((decision > scoreMin-1e-8)&(decision < scoreMin+1e-8))[0]; idxMin
array([ 70, 126, 138])
X_f2[idxMin]
array([[4.8, 1.8],
[4.8, 1.8],
[4.8, 1.8]])
y[idxMin], pre[idxMin]
(array([0, 1, 1]), array([1, 1, 1]))
logit2.score(X_f2, y)
0.9666666666666667
위 결과에 의하면 독립변수가 [4.8, 1.8] 이상일 경우 1이 되므로 이 값이 결정경계가 됩니다.
주가 자료에 적용
다음은 FinanceDataReader.DataReader() 함수를 호출한 환율과 미국의 주요지수, 코스피 지수의 일간변화율에 대해 로지스틱회귀를 적용하여 봅니다. 이 데이터 세트의 변수 중 코스피 "Close"의 상승(1)과 하락(0)을 반응변수로 합니다.
st=pd.Timestamp(2010,3, 1)
et=pd.Timestamp(2022, 6, 8)
nme={'exchg':'USD/KRW','dj':"DJI",'nasd':'IXIC','soxx':"SOXX", 'vix':'VIX'}
kos=fdr.DataReader('KS11', st, et)
stock={}
for i, j in zip(nme.keys(), nme.values()):
stock[i]=fdr.DataReader(j, st, et)['Close']
stock=pd.DataFrame(stock.values(), index=stock.keys()).T
호출한 자료는 결측값을 포함할 수 있습니다. 특히나 DataReader()함수에 의해 호출한 자료는 결측값을 0으로 처리된 형태입니다. 이 결측값을 처리하기 위해 Pandas 라이브러리의 .mask(), .ffill() 메서드를 적용하였습니다.
Stock=pd.concat([kos, stock], join='outer', axis=1) Stock=Stock.mask(Stock==0) Stock=Stock.ffill().dropna() Stock.tail(3)
| Close | Open | High | Low | Volume | Change | exchg | dj | nasd | soxx | vix | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Date | |||||||||||
| 2022-06-06 | 2670.65 | 2679.57 | 2681.51 | 2663.00 | 562470000.0 | 0.0044 | 1254.45 | 32914.20 | 12061.37 | 421.74 | 25.07 |
| 2022-06-07 | 2626.34 | 2659.27 | 2662.04 | 2622.78 | 524710000.0 | -0.0166 | 1254.36 | 33179.81 | 12175.23 | 425.11 | 24.02 |
| 2022-06-08 | 2626.15 | 2633.53 | 2639.52 | 2621.96 | 560660.0 | -0.0001 | 1257.28 | 32910.38 | 12086.27 | 416.06 | 23.96 |
data=Stock.pct_change().dropna() data=data.mask(data==0).ffill() data.head()
| Close | Open | High | Low | Volume | Change | exchg | dj | nasd | soxx | vix | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Date | |||||||||||
| 2010-03-03 | 0.004532 | 0.001749 | 0.001531 | 0.001479 | -0.028017 | -0.651163 | -0.002402 | -0.000886 | -0.000048 | -0.004841 | -0.012067 |
| 2010-03-04 | -0.002613 | 0.006253 | 0.005220 | 0.000695 | -0.050636 | -1.577778 | 0.001532 | 0.004557 | 0.005099 | -0.000846 | -0.005842 |
| 2010-03-05 | 0.010116 | -0.000338 | 0.002127 | 0.007609 | -0.101713 | -4.884615 | -0.003059 | 0.011687 | 0.014850 | 0.011219 | -0.069444 |
| 2010-03-08 | 0.015582 | 0.017688 | 0.015729 | 0.015344 | 0.068967 | 0.544554 | -0.007014 | -0.001295 | 0.002519 | -0.000419 | 0.021240 |
| 2010-03-09 | 0.000476 | 0.002661 | 0.000590 | 0.002613 | 0.725646 | -0.967949 | 0.002349 | 0.001124 | 0.003632 | -0.002723 | 0.007307 |
np.where(data==0)
(array([], dtype=int64), array([], dtype=int64))
위 객체 data에서 반응변수를 정의합니다.
y=data["Close"][1:] y=(y>=0).astype(int) y.head(2)
Date
2010-03-04 0
2010-03-05 1
Name: Close, dtype: int64
위 객체 Stock은 변수별로 데이터의 scale에 차이를 보입니다. 그러므로 표준화를 실시하고 정제된 데이터와 반응변수를 훈련과 검정 세트로 분리합니다.
from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split
ind=data.iloc[:, 1:].values scalerF=StandardScaler().fit(ind) feat=scalerF.transform(ind) x=feat[:-1,:] new=feat[-1,:]
x.shape, y.shape
((3201, 10), (3201,))
xtr, xte, ytr, yte=train_test_split(x, y, test_size=0.3, random_state=1)
로지스틱 모델을 구축하고 훈련과 검정 세트에서의 score(R2)를 산출하였습니다.
logit_stock=LogisticRegression(random_state=2).fit(xtr, ytr) logit_stock.score(xtr, ytr),logit_stock.score(xte, yte)
(0.6200892857142857, 0.6066597294484911)
preXtr=logit_stock.predict(xtr)
결정경계를 설정하기 위해 decision_function()메서드를 적용하였으며 그 결과 중 절대값이 최소가 되는 위치의 독립변수와 관찰치의 값, 그리고 예측치의 값을 확인합니다. 그 지점의 관찰치와 예측치 모두 1이며 구현된 모델은 이 지점에서의 독립변수 값 보다 증가한 부분에서는 상승(1)을 추정합니다.
decis=logit_stock.decision_function(xtr) decisMin=np.abs(decis).min() decisIdx=np.where((decis > decisMin-1e-5)&(decis < decisMin+1e-5))[0] decisIdx
array([1061])
xtr[decisIdx], ytr[decisIdx], preXtr[decisIdx]
(array([[ 3.88273832, 1.34537194, 2.64276259, -0.03847125, 0.05802377,
-4.4652264 , -1.15398721, -0.77372962, -0.97823835, 0.92170431]]),
Date
2010-05-31 1
Name: Close, dtype: int64,
array([1]))
deciBoundary=pd.DataFrame(scalerF.inverse_transform(xtr[decisIdx]), columns=data.columns[1:]) deciBoundary
| Open | High | Low | Volume | Change | exchg | dj | nasd | soxx | vix | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.04123 | 0.011535 | 0.026421 | -0.116512 | -0.40625 | -0.023653 | -0.011927 | -0.009062 | -0.016079 | 0.080526 |
예측을 위해 작성된 변수 new에 대한 예측을 실시합니다.
new
array([-0.94171816, -1.03808245, -0.06060143, -0.06166036, 0.00678508,
0.43284885, -0.79877753, -0.63374162, -1.27859125, -0.06452137])
logit_stock.predict(new.reshape(1,-1))
array([0])