Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

선형 회귀

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

1최적화

Data Generation Process

y=y0+ϵ\begin{aligned} y = y_0 + \epsilon \end{aligned}

Model

y^=f(x;θ)\hat{y} = f(x; \theta)

최적 매개변수

Loss=1nin(yiy^i)2θ=arg minθLoss(y,y^)\begin{aligned} &\text{Loss} = \frac{1}{n}\sum_{i}^{n}\left({y^i - \hat{y}^{i}}\right)^2 \\ &\theta^{*} = \argmin_{\theta}{\text{Loss}(y, \hat{y})} \end{aligned}

특성수가 dd인 데이터포인트를 xiRdx^i \in \mathbb{R}^{d}라고 할 때, 데이터포인트의 집합인 데이터셋 {xi}i=1n\{x^i\}_{i=1}^{n}은 행렬 XRn×dX \in \mathbb{R}^{n \times d}로 표현할 수 있습니다. 각 데이터포인트의 목표값은 스칼라값 yiRy^i \in \mathbb{R}이며, 목표값 집합 yRn y \in \mathbb{R}^{n}입니다. 선형 모형의 가중치는 wRdw \in \mathbb{R}^{d}입니다.

목표는 손실을 최소화하는 매개변수를 찾는 것입니다. 즉, 손실은 매개변수의 함수이기 때문에 목표 (objective) 함수는 다음과 같습니다.

J(w)=i=1n(yiwxi)2J(w) = \sum_{i=1}^{n}(y^i - w^{\top}x^i)^{2}
J(w)=i=1n(yiwxi)2=yXw2=(yXw)(yXw)=(ywX)(yXw)=yyyXwwXy+wXXw=yy2yXw+wXXw\begin{align} J(w) &= \sum_{i=1}^{n}(y^i - w^{\top}x^i)^{2} \\ &= \| y - Xw \|^{2} \\ &= (y - Xw)^{\top}(y - Xw) \\ &= (y^{\top} - w^{\top}X^{\top})(y - Xw) \\ &= y^{\top}y - y^{\top}Xw - w^{\top}X^{\top}y + w^{\top}X^{\top}Xw \\ &= y^{\top}y - 2y^{\top}Xw + w^{\top}X^{\top}Xw \end{align}

손실은 2차항의 볼록 함수이기 때문에 미분이 0이 되는 지점이 최소값입니다.

wJ(w)=(Jw1,Jw2,Jwd)0=2Xy+2XXwXXw=Xy\begin{align} \nabla_{w}{J}(w) &= \left( \frac{\partial{J}}{\partial{w_1}}, \frac{\partial{J}}{\partial{w_2}}, \dots \frac{\partial{J}}{\partial{w_d}} \right) \\ 0 &=-2X^{\top}y+2X^{\top}Xw \\ X^{\top}Xw &= X^{\top}y \end{align}

따라서, 미분이 0이 되는 지점에서 최적해 ww는 다음과 같이 계산됩니다.

wJ(w)=0XXw=Xy\nabla_w J(w) = 0 \quad \Longrightarrow \quad X^{\top}Xw = X^{\top}y

XXX^{\top}X의 역행렬이 존재하는 경우, 아래와 같이 ww를 구할 수 있습니다.

w=(XX)1Xy\begin{align} w & = (X^{\top}X)^{-1}X^{\top} y \end{align}

XXX^{\top}X의 역행렬이 존재하지 않는 특이 행렬(singular matrix)의 경우(즉, det(XX)=0\det(X^{\top}X)=0)를 포함해, 선형회귀의 최적해는 일반적으로 무어-펜로즈 유사역행렬 (Moore-Penrose pseudoinverse)로 구할 수 있습니다.

w=X+yw = X^{+}y
# DGP: Data Generating Process. "진짜 원인"
xs = np.linspace(0, 10, 100)
noise = np.random.default_rng(seed=4).normal(0, 1, size=xs.shape)
ys = 1.2 * xs + 3.4 + noise

model = lambda x, w, b: w * x + b # 선형 모델
params = [{'w': -1.0, 'b': 15.0}] # 초기 추정

# 최적해 방정식: 정규 방정식: X.T @ X @ w = X.T @ y
X = np.stack([np.ones_like(xs), xs], axis=1) # [x0=1, x1=xs]
best_params = np.linalg.inv(X.T @ X) @ X.T @ ys
params.append({'w': best_params[1], 'b': best_params[0]}) # 최적해

outputs = []
for 매개변수 in params:
    y_pred = model(xs, **매개변수)
    outputs.append(y_pred)

평균제곱오차 = lambda y_true, y_pred: np.mean((y_true - y_pred) ** 2)
손실 = [평균제곱오차(ys, y_pred) for y_pred in outputs]
display(pd.DataFrame(params).assign(손실=손실).round(3))

plt.scatter(xs, ys)
plt.plot(xs, ys - noise, color='red', linestyle='--', label='DGP')
for (w_and_b, outputs) in zip(params, outputs):
    라벨 = f"w={w_and_b['w']}, b={w_and_b['b']}"
    plt.plot(xs, outputs, label=라벨, color='black')
plt.legend()
plt.show()

2최대우도 추정

선형 회귀의 최적 매개변수는 최대우도추정(Maximum Likelihood Estimation; MLE) 관점에서 정당화 될 수 있습니다.

yi=wxi+ϵi\begin{aligned} y^{i} = w^{\top}x^{i} + \epsilon^{i} \end{aligned}

관측 오류 ϵ\epsilon에 대해 표준 정규 분포 ϵN(0,σ2)\epsilon \sim \mathcal{N}(0, \sigma^2) 를 가정하면, 선형 회귀의 출력은 정규 분포 yN(wxi,σ2) y \sim \mathcal{N}(w^{\top}x^{i}, \sigma^{2})의 확률 분포를 가집니다. 즉, 기대값(평균) μ=E[yx]=wx\mu = \mathbb{E}[y \mid x] = w^{\top}x

$$

N(yμ,σ)=12πσ2exp(12σ2(yμ)2)=12πσ2exp(12σ2(ywx)2)\begin{align} \mathcal{N}(y | \mu, \sigma) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{1}{2\sigma^2} (y - \mu)^2 \right) \\ &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{1}{2\sigma^2} (y - w^{\top}x)^2 \right) \end{align}

$$

따라서, 선형 회귀는 입력 xix^i와 매개변수 ww가 주어졌을 때 출력 yiy^i의 조건부 확률 분포 p(yixi,w,σ2)p(y^i \mid x^i, w, \sigma^2)를 모델링합니다.

전체 데이터에 대한 우도(likelihood)는 다음과 같습니다. 우도(likelihood)는 매개변수 ww에 의해 결정되기 때문에 이 문제는 매개변수 ww에 대한 최대우도추정의 문제가 됩니다.

$$ \begin{aligned} \mathcal{L}(w) &= \prod_{i=1}^n p(y^i \mid x^i,w, \sigma^2)

\ &= \prod_{i=1}^n \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{1}{2\sigma^2} (y^i - w^{\top}x^i)^2 \right) \end{aligned} $$

확률값은 곱해지면 작아지기 때문에 계산의 편의를 위해 로그를 취해 로그 우도 (log likelihood)로 계산합니다. 로그로 변환하면, 곱하기는 더하기가 됩니다.

logL(w)=i=1nlogp(yixi,w,σ2)=i=1n[12log(2πσ2)12σ2(yiwxi)2]=n2log(2πσ2)12σ2i=1n(yiwxi)2\begin{aligned} \log \mathcal{L}(w) &= \sum_{i=1}^n \log p(y^i \mid x^i, w, \sigma^2) \\ &= \sum_{i=1}^n \left[ -\frac{1}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2}(y^i - w^\top x^i)^2 \right] \\ &= -\frac{n}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2} \sum_{i=1}^n (y^i - w^\top x^i)^2 \end{aligned}

상수항을 제외하면 로그 우도는 다음과 같이 정리할 수 있습니다.

logL(w)i=1n(yiwxi)2\log \mathcal{L}(w) \propto - \sum_{i=1}^n \left(y^i - w^{\top}x^i \right)^2

따라서, 로그우도를 최대화하는 것은 오차제곱합을 최소화하는 것과 같습니다.

arg maxlogL(w)=arg mini=1n(yiwxi)2\argmax \log \mathcal{L}(w) = \argmin \sum_{i=1}^{n} \left(y^i - w^{\top}x^i \right)^2