<jupyter_start><jupyter_text>로지스틱 회귀 분석 로지스틱 회귀(Logistic Regression) 분석은 회귀 분석이라는 명칭을 가지고 있지만 분류(classsification) 방법의 일종이다.로지스틱 회귀 모형에서는 베르누이 확률 변수(Bernoilli random variable)의 모수(parameter) $\theta$가 독립 변수 $x$에 의존한다고 가정한다.$$ p(y \mid x, \theta) = \text{Ber} (y \mid  \theta(x) )$$여기에서 모수 $\theta$ 는 0과 1사이의 실수이며 다음과 같이 $x$의 값에 의존하는 함수이다.$$\theta = f(w^Tx)$$  시그모이드 함수 모수 $\theta$는 일반적인 회귀 분석의 종속 변수와 달리 0 부터 1까지의 실수값만 가질 수 있기 때문에 시그모이드 함수(sigmoid function)이라 불리는 특별한 형태의 함수 $f$를 사용해야 한다.시그모이드 함수는 종속 변수의 모든 실수 값에 대해 유한한 구간 $(a,b)$ 사이의 한정된(bounded) 값과 양의 기울기를 가지는 함수를 말하며 다음과 같은 함수들이 주로 사용된다. * 로지스틱 함수 (Logistic Function)  $$ \text{logitstic}(z) = \dfrac{1}{1+\exp{(-z)}} $$  * 오차 함수 (Error Function)  $$ \text{erf}(z) = \frac{2}{\sqrt\pi}\int_0^z e^{-t^2}\,dt $$  * 하이퍼볼릭 탄젠트 함수 (Hyperbolic tangent)  $$ \tanh(z) = \frac{\sinh z}{\cosh z} = \frac {e^z - e^{-z}} {e^z + e^{-z}} $$  * 역 탄젠트 함수 (Arc-tangent)  $$ \arctan(z) = \tan^{-1}(z) $$<jupyter_code>xx = np.linspace(-10, 10, 1000)
plt.plot(xx, (1/(1+np.exp(-xx)))*2-1, label="logistic (scaled)")
plt.plot(xx, sp.special.erf(0.5*np.sqrt(np.pi)*xx), label="erf (scaled)")
plt.plot(xx, np.tanh(xx), label="tanh")
plt.ylim([-1.1, 1.1])
plt.legend(loc=2)
plt.show()<jupyter_output><empty_output><jupyter_text>로지스틱 함수 여러가지 시그모이드 중 로지스틱 함수는 다음과 같은 물리적인 의미를 부여할 수 있기 때문에 많이 사용된다.우선 Bernoulli 시도에서 1이 나올 확률 $\theta$ 과 0이 나올 확률 $1-\theta$ 의 비(ratio)는 다음과 같은 수식이 되며 odds ratio 라고 한다.$$ \text{odds ratio} = \dfrac{\theta}{1-\theta} $$이 odds ratio 를 로그 변환한 것이 로지트 함수(Logit function)이다.$$ z = \text{logit}(\text{odds ratio}) = \log \left(\dfrac{\theta}{1-\theta}\right) $$ 로지스틱 함수(Logistic function) 는 이 로지트 함수의 역함수이다.$$ \text{logitstic}(z) = \theta(z) = \dfrac{1}{1+\exp{(-z)}} $$    로지스틱 모형의 모수 추정 로지스틱 모형은 일종의 비선형 회귀 모형이지만 다음과 같이 MLE(Maximum Likelihood Estimation) 방법으로 모수 $w$를 추정할 수 있다.여기에서는 종속 변수 $y$가 베르누이 확률 변수라고 가정한다.$$ p(y \mid x, \theta) = \text{Ber} (y \mid  \theta(x) )$$ 데이터 표본이 $\{ x_i, y_i \}$일 경우 Log Likelihood $\text{LL}$ 를 구하면 다음과 같다.$$\begin{eqnarray}\text{LL} &=& \log \prod_{i=1}^N \theta_i(x_i)^{y_i} (1-\theta_i(x_i))^{1-y_i} \\&=& \sum_{i=1}^N \left( y_i \log\theta_i(x_i) +  (1-y_i)\log(1-\theta_i(x_i)) \right) \\\end{eqnarray}$$$\theta$가 로지스틱 함수 형태로 표현된다면$$\log \left(\dfrac{\theta(x)}{1-\theta(x)}\right) = w^T x$$$$\theta(x) = \dfrac{1}{1 + \exp{(-w^Tx)}}$$가 되고 이를  Log Likelihood 에 적용하면 다음과 같다.$$\begin{eqnarray}\text{LL} &=& \sum_{i=1}^N \left( y_i \log\theta_i(x_i) +  (1-y_i)\log(1-\theta_i(x_i)) \right) \\&=& \sum_{i=1}^N \left( y_i \log\left(\dfrac{1}{1 + \exp{(-w^Tx_i)}}\right) - (1-y_i)\log\left(\dfrac{\exp{(-w^Tx_i)}}{1 + \exp{(-w^Tx_i)}}\right) \right) \\\end{eqnarray}$$ 이 값의 최대화하는 값을 구하기 위해 chain rule를 사용하여 $w$로 미분해야 한다.우선 $\theta$를 $w$로 미분하면$$ \dfrac{\partial \theta}{\partial w} = \dfrac{\partial}{\partial w} \dfrac{1}{1 + \exp{(-w^Tx)}} \ = \dfrac{\exp{(-w^Tx)}}{(1 + \exp{(-w^Tx)})^2} x \ = \theta(1-\theta) x $$chain rule를 적용하면 $$ \begin{eqnarray}\dfrac{\partial \text{LL}}{\partial w} &=& \sum_{i=1}^N \left( y_i \dfrac{1}{\theta_i(x_i;w)} - (1-y_i)\dfrac{1}{1-\theta_i(x_i;w)} \right) \dfrac{\partial \theta}{\partial w}  \\&=& \sum_{i=1}^N \big( y_i (1-\theta_i(x_i;w)) - (1-y_i)\theta_i(x_i;w)  \big)  x_i \\&=& \sum_{i=1}^N \big( y_i  - \theta_i(x_i;w) \big) x_i \\\end{eqnarray}$$이 값은 $w$에 대한 비선형 함수이므로 선형 모형과 같이 간단하게 그레디언트가 0이 되는 모수 $w$ 값에 대한 수식을 구할 수 없으며 수치적인 최적화 방법(numerical optimization)을 통해 최적 모수 $w$의 값을 구해야 한다.  수치적 최적화 단순한 Steepest Gradient 방법을 사용한다면 최적화 알고리즘은 다음과 같다.그레디언트 벡터는$$g_k = \dfrac{d}{dw}(-LL)$$이 방향으로 step size $\eta_k$ 만큼 움직이면 다음과 같이 반복적으로 최적 모수값을 구할 수 있다.$$\begin{eqnarray}w_{k+1} &=& w_{k} - \eta_k g_k \\&=& w_{k} + \eta_k \sum_{i=1}^N \big( y_i  - \theta_i(x_i) \big) x_i\\\end{eqnarray}$$  Scikit-Learn 패키지의 로지스틱 회귀 Scikit-Learn 패키지는 로지스틱 회귀 모형 `LogisticRegression` 를 제공한다.<jupyter_code>from sklearn.datasets import make_classification
X0, y = make_classification(n_features=1, n_redundant=0, n_informative=1, n_clusters_per_class=1, random_state=4)
X = sm.add_constant(X0)
from sklearn.linear_model import LogisticRegression
model = LogisticRegression().fit(X0, y)
xx = np.linspace(-3, 3, 100)
sigm = 1.0/(1 + np.exp(-model.coef_[0][0]*xx - model.intercept_[0]))
plt.plot(xx, sigm)
plt.scatter(X0, y, marker='o', c=y, s=100)
plt.scatter(X0, model.predict(X0), marker='x', c=y, s=200, lw=2, alpha=0.5, cmap=mpl.cm.jet)
plt.xlim(-3, 3)
plt.show()<jupyter_output><empty_output><jupyter_text>statsmodels 패키지의 로지스틱 회귀 statsmodels 패키지는 로지스틱 회귀 모형 `Logit` 를 제공한다. 사용방법은 `OLS` 와 동일하다. Scikit-Learn 패키지와 달리 `Logit` 클래스는 classification 되기 전의 값을 출력한다<jupyter_code>logit_mod = sm.Logit(y, X)
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())
xx = np.linspace(-3, 3, 100)
sigmoid = logit_res.predict(sm.add_constant(xx))
plt.plot(xx, sigmoid, lw=5, alpha=0.5)
plt.scatter(X0, y, marker='o', c=y, s=100)
plt.scatter(X0, logit_res.predict(X), marker='x', c=y, s=200, lw=2, alpha=0.5, cmap=mpl.cm.jet)
plt.xlim(-3, 3)
plt.show()<jupyter_output><empty_output><jupyter_text>예제 1: Michelin and Zagat 가이드 비교 다음 데이터는 뉴욕시의 레스토랑에 대한 두 개의 가이드북에서 발취한 것이다.* Food: Zagat Survey 2006 의 고객 평가 점수* InMichelin: 해당 고객 평가 점수를 받은 레스토랑 중 2006 Michelin Guide New York City 에 실린 레스토랑의 수* NotInMichelin: 해당 고객 평가 점수를 받은 레스토랑 중 2006 Michelin Guide New York City 에 실리지 않은 레스토랑의 수* mi: 해당 고객 평가 점수를 받은 레스토랑의 수* proportion: 해당 고객 평가 점수를 받은 레스토랑 중 2006 Michelin Guide New York City 에 실린 레스토랑의 비율<jupyter_code>df = pd.read_table("~/data/sheather/MichelinFood.txt")
df
df.plot(kind="scatter", x="Food", y="proportion", s=100)
plt.show()
X = sm.add_constant(df.Food)
y = df.proportion
model = sm.Logit(y, X)
result = model.fit()
print(result.summary())
df.plot(kind="scatter", x="Food", y="proportion", s=50, alpha=0.5)
xx = np.linspace(10, 35, 100)
plt.plot(xx, result.predict(sm.add_constant(xx)), "r", lw=4)
plt.xlim(10, 35)
plt.show()<jupyter_output><empty_output><jupyter_text>예제 2: Michelin 가이드 예측 다음 데이터는 뉴욕시의 개별 레스토랑의 고객 평가 점수와 Michelin 가이드 수록 여부를 보인 것이다.* InMichelin: Michelin 가이드 수록 여부* Restaurant Name: 레스토랑 이름* Food: 식사에 대한 고객 평가 점수 (1~30)* Decor: 인테리어에 대한 고객 평가 점수 (1~30)* Service: 서비스에 대한 고객 평가 점수 (1~30)* Price: 저녁 식사 가격 (US$)<jupyter_code>df = pd.read_csv("~/data/sheather/MichelinNY.csv")
df.tail()
sns.stripplot(x="Food", y="InMichelin", data=df, jitter=True, orient='h', order=[1, 0])
plt.grid(True)
plt.show()
X = sm.add_constant(df.Food)
y = df.InMichelin
model = sm.Logit(y, X)
result = model.fit()
print(result.summary())
xx = np.linspace(10, 35, 100)
pred = result.predict(sm.add_constant(xx))
decision_value = xx[np.argmax(pred > 0.5)]
print(decision_value)
plt.plot(xx, pred, "r", lw=4)
plt.axvline(decision_value)
plt.xlim(10, 35)
plt.show()<jupyter_output>21.8686868687<jupyter_text>예제 3: Fair's Affair Dataset<jupyter_code>print(sm.datasets.fair.SOURCE)
print(sm.datasets.fair.NOTE)
df = sm.datasets.fair.load_pandas().data
df.head()
sns.factorplot(x="affairs", y="children", row="yrs_married", data=df,
               orient="h", size=2, aspect=5, kind="box")
plt.show()
df['affair'] = (df['affairs'] > 0).astype(float)
modoel = smf.logit("affair ~ rate_marriage + religious + yrs_married + age + educ + children", df).fit()
print(modoel.summary())<jupyter_output>Optimization terminated successfully.
         Current function value: 0.547174
         Iterations 6
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                 affair   No. Observations:                 6366
Model:                          Logit   Df Residuals:                     6359
Method:                           MLE   Df Model:                            6
Date:                Mon, 13 Jun 2016   Pseudo R-squ.:                  0.1297
Time:                        08:05:46   Log-Likelihood:                -3483.3
converged:                       True   LL-Null:                       -4002.5
                                        LLR p-value:                4.345e-221
=================================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
-----------------------[...]
{'fraction_non_alphanumeric': 0.2022224935889608, 'fraction_numerical': 0.0319941384784467, 'mean_word_length': 3.507429829389103, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '252388', 'license': 'MIT', 'n_tokens_mistral': 4468, 'n_tokens_neox': 4868, 'n_words': 845, 'path': '17. 로지스틱 회귀 분석/01. 로지스틱 회귀 분석.ipynb', 'repo_name': 'zzsza/Datascience_School', 'size': '8189'}