29. Piecewise Polynomial(Constant, Linear) Regression에 대해서 알아보자 with Python

---

이번 포스팅에서는 Piecewise Polynomial(Constant, Linear) Regression에 대한 개념과 이를 파이썬으로 구현하는 방법에 대해서 알아보려고 한다.

 

- 목차 -

1. Piecewise Polynomial Regression 이란 무엇인가?

2. Piecewise Polynomial Regression 알고리즘

3. 장단점

4. 파이썬 구현

 

---

   1. Piecewise Polynomial Regression 이란 무엇인가?

1) 정의

Piecewise Polynomial Regression(PPR)은 주어진 매듭점으로 나누어진 영역에서 설명 변수의 다항 함수로 회귀 함수를 추정하는 방법을 말한다.

---

2) 파헤치기

앞에서 살펴본 정의를 하나하나 살펴보자. 이때 설명과 구현의 편의를 위하여 설명 변수는 하나인 회귀 문제만을 고려한다.

a. PPR은 먼저 주어진 매듭점으로 데이터 영역을 나눈다. 

먼저 데이터 $(x_i, y_i), i=1, \ldots, n$이 있다고 하자. 그리고 $K$개의 매듭점 $(\tau_1, \ldots, \tau_K)$이 주어졌다. PPR은 먼저 주어진 매듭점을 이용하여 $K+1$개의 영역을 나누게 된다. 아래 그림은 $K=3$인 경우에 데이터 영역을 나눈 것을 나타낸 것이다.

 

여기서 각 매듭점 사이의 거리가 반드시 등간격($\tau_2-\tau_1 = \tau_3-\tau_2$)일 필요는 없다. 또한 매듭점은 데이터 범위, 즉 설명 변수의 최소값과 최대값 사이에 있어야 하며 매듭점이 반드시 데이터에 포함될 필요는 없다. 즉, 모든 $k$에 대하여

$$\tau_k = x_i \;\; \text{for some} \;\; i=1, \ldots, n$$

일 필요는 없다는 것이다.

b. PPR은 각 영역에서 다항 함수로 적합하여 회귀 함수를 적합한다.

말 그대로 $K+1$개의 영역에 대해서 다항 회귀 함수를 적합한다. 아래 그림은 $K=3$인 경우 4($=K+1$) 개의 영역에 대하여 일차 함수를 적합한 것을 나타낸다.

이제 이를 엄밀하게 수식으로 표현해보자. PPR은 $K+1$개의 영역에 대하여 다음과 같이 $p$차 다항식을 적합한다. 즉, $k=1, \ldots, K+1$에 대하여

$$f(x) = b_{0k} + b_{1k}x + b_{2k}x^2 + \cdots + b_{pk}x^p, \;\; \tau_{k-1} \leq x < \tau_{k}\tag{1}$$

을 적합하게 된다. 이때 $\tau_0 = -\infty, \tau_{K+1} = \infty$이다. 

 

그러면 모든 영역에 대해서 매번 다항 함수를 적합해야 하는가?...

 

아니다. 지시 함수(Indicator Function)을 이용하면 한 번의 적합만으로 마치 각 영역에 대해서 적합한 것과 같은 효과를 얻게 된다. 이를 위해 $f_k$를 다음과 같이 정의하자.

$$f_k(x) = \left( b_{0k} + b_{1k}x + b_{2k}x^2 + \cdots + b_{pk}x^p \right )I(\tau_{k-1}\leq x < \tau_k), k=1, \ldots, K+1$$

그러면 (1)에서 $f$는 다음과 같이 표현된다.

$$f(x) = f_1(x) + f_2(x) + \cdots + f_{K+1}(x)\tag{2}$$

 

이를 이용하면 다음과 같이 오차 제곱합을 최소화하는 한 번에 회귀 계수 $b = (b_{01}, \ldots,b_{p1}, \ldots,  b_{1,K+1}, \ldots, b_{p,K+1})$ 를 추정할 수 있게 된다.

$$\sum_{i=1}^n(y_i-f(x_i))^2 = \|y-X\beta \|^2 \tag{3}$$

이때 $X$는 $n\times (p+1)(K+1)$인 설계 행렬이며 다음과 같은 Block 행렬로 이루어져 있다.

$$X = \begin{pmatrix} X^1 & X^2 & \cdots & X^{K+1} \end{pmatrix} $$

여기서 $X^k, k=1,\ldots, K+1$는 $n\times (p+1)$ 행렬로 $X^k_{ij} = x_i^{j-1}I(\tau_{k-1} \leq x_i < \tau_k)$이다.

 

이때 회귀 계수에 대한 해석은 구간 별로 일반적인 회귀 계수의 해석과 동일하게 할 수 있다. 예를 들어 구간 별로 일차 모형을 적합했을 경우 $b_{1k}$는 구간 $[\tau_{k-1}, \tau_k)$에서 $x$가 한 단위 변할 때 $y$의 평균 변화율로 해석할 수 있다는 것이다.

---

3) 함수의 연속 조건

일반적으로 아무런 조건 없이 최소 제곱법으로 추정하게 되면 매듭점에서 불연속 함수가 된다. 불연속이 되면 매듭점에서의 추정치를 매우 신뢰할 수 없게 된다.

 

따라서 매듭점에서 연속이 되도록 (3)에 다음과 같은 제약식을 주어야 한다.

$$f_k(\tau_k) = f_{k+1}(\tau_k), k=1, \ldots, K \tag{4}$$

 

아래 그림은 PPR 방법으로 추정한 회귀 함수가 연속이 되기 위한 조건을 나타낸 그림이다.

연속이 되기 위한 조건

(3)과 (4)를 종합하면 우리가 풀고자 하는 문제는 결국 등호 제약식이 포함된 오차 제곱합을 최소화하는 $b$를 찾는 것이다.

$$\DeclareMathOperator*{\argminA}{arg\,min} \argminA_b \|y-Xb \|^2 \;\;\text{with constraint} \;\; Ab = 0 \tag{5}$$

이때 $A$는 $K\times (p+1)(K+1)$ 행렬이며 $A$의 $k$ 번째 행은 $(p+1)(k-1)+1$번째 부터 $(p+1)(k+1)$번째 원소가 $(1, \tau_k, \ldots, \tau_k^p, -1, -\tau_k, \ldots, -\tau_k^p)$이고 다른 곳에서는 모두 0이다.

---

4) 차수 결정

차수는 데이터의 분포를 시각적으로 확인한 뒤에 주관적으로 결정할 수도 있고 후보군 $p\in \{1, 2, \ldots, 5 \}$에서 교차 검증을 이용하여 결정할 수 있다. 하지만 보통 3을 넘지는 않는 것 같다. 

---

5) 매듭점 결정

지금까지는 매듭점이 주어진 경우를 생각했다. 하지만 매듭점을 자동으로 결정하는 방법도 많이 연구되어 있다. 매듭점 결정하는 방법은 다음과 같은 경우로 나누어 생각할 수 있다.

---

(1) 매듭점의 개수를 알지만 위치를 모르는 경우

(2) 매듭점의 개수와 위치를 모두 모르는 경우

---

이제 각 경우에 대해서 매듭점을 결정하는 방법을 알아보자.

(1) 매듭점의 개수를 알지만 위치를 모르는 경우

만약 매듭점의 개수 $K$를 아는 경우 전진 선택법(Forward Selection) 또는 후진 소거법(Back Elimination)을 이용할 수 있다.

 

전진 선택법은 먼저 매듭점의 후보군을 데이터의 유니크한 값에서 (다항 차수에 따라서) 왼쪽 부분과 오른쪽 데이터를 적절히 제외한 나머지로 설정한다. 다음으로 가장 적합도가 좋은 매듭점을 하나씩 추가하는 방식으로 $K$개까지 매듭점을 지정하게 된다. 적합도를 측정하는 측도로는 보통 잔차 제곱합을 사용하게 된다.

 

예를 들어 $K=2$이고 아래 그림과 같이 유니크한 값의 개수가 10개라고 해보자. 이 경우 Piecewise Constant, 즉 다항 차수가 0인 경우에는 $x_2$부터 $x_9$까지 후보군이 된다(왼쪽 그림). 왜냐하면 Constant를 추정하기 위해선 최소 한 개 데이터가 필요하기 때문이다. 만약 Piecewise Linear, 즉 다항 차수가 1인 경우에는 $x_3$부터 $x_8$까지 후보군이 된다(오른쪽 그림). 이는 직선을 추정하기 위해서 최소 2개의 점이 필요하기 때문이다.

 

 

Piecewise Linear로 적합한다고 가정하고 $x_3$부터 $x_8$까지 하나씩 매듭점으로 정하여 모형을 적합한 후 가장 적합력이 좋은 매듭점을 설정한다. 이때 다음과 같이 $x_4$에서 가장 적합력이 좋다면 $\tau_1 = x_4$로 정한다.

이제 $x_4$를 제외한 나머지 점에 대하여 가장 적합도가 좋은 매듭점을 찾는다. 예를 들어 $x_7$이 가장 좋은 적합력을 가졌다면 $\tau_2 = x_7$로 설정한다. $K=2$이므로 $\tau_1, \tau_2$를 최종 매듭점으로 정한다. 이때 각 구간의 최소 샘플 수를 지정할 수도 있다.

 

후진 소거법은 후보군을 모두 매듭점으로 하여 적합한 다음 적합도가 가장 안 좋은 점을 제외하게 된다. 이를 매듭점 개수가 $K$개가 될 때까지 반복하는 것이다.

(2) 매듭점의 개수와 위치를 모두 모르는 경우

매듭점의 개수를 모르는 경우 매듭점의 개수 후보를 정하고 각 개수에 대하여 (1)의 방법을 이용한다. 이 과정에서 최종 매듭점의 위치와 이 매듭점에 대한 적합도가 저장되어 있다. 다음으로 각 개수를 x축, 그에 따른 적합력을 y축으로 하여 Elbow Plot(Scree Plot)을 그려보고 적합력의 변화가 크지 않은 점을 기준으로 매듭점의 개수 결정하고 이때 저장되어 있던 위치를 최종 매듭점으로 한다.

 

예를 들어 아래 그림에서 3, 4 즈음에서 적합력의 감소 정도가 낮으므로 3 또는 4를 매듭점의 개수로 설정한다.

---

---

여기서 궁금한 점 ...

매듭점의 개수를 지정할 때 다항 차수와 최종 모형을 적합할 때 다항 차수가 꼭 같아야 하는지에 대한 자료를 찾지 못했다. 하지만 꼭 같을 필요는 없을 듯하다. 매듭점 개수를 정할 때에는 직선 즉, 다항 차수를 1로 하고 모형 적합시에는 다항 차수를 3으로 할 수도 있다는 것이다.

---

---

   2. Piecewise Polynomial Regression 알고리즘

PPR 알고리즘은 다음과 같다.

 

Algorithm

 

1 단계)

주어진 매듭점 $\tau_1, \ldots, \tau_K$, 다항 차수$p$에 대하여 $n\times (p+1)(K+1)$ 설계 행렬

$$X = \begin{pmatrix} X^1 & X^2 & \cdots & X^{K+1} \end{pmatrix} $$

을 구성한다. 이때 $X^k_{ij} = x_i^{j-1}I(\tau_{k-1} \leq x_i < \tau_k)$이다.

 

2 단계)

(3)을 최소화하는 회귀 계수 벡터 $b$를 추정한다. 만약 연속을 가정한다면 (5)를 만족하는 $b$를 찾는다.

 

3 단계)

$x$가 주어진 경우 $x$를 포함하는 매듭점 구간이 $[\tau_{k-1}, \tau_k)$이라면 $b_{0k} + b_{1k}x + \cdots + b_{pk}x^p$로 예측한다.

 

---

   3. 장단점

- 장점 -

a. 특정 구간별로 패턴을 잡아 낼 수 있다.

 

b. 구간별로 회귀 계수의 해석이 가능하다.

구간별로 회귀 계수의 해석이 어려운 경우도 있다(B-Spline)

 

c. 구현이 어렵지 않다.

PPR 추정 문제를 제약식이 있는 최소 제곱법으로 바꾸면 어렵지 않게 모형을 추정할 수 있다.

---

- 단점 -

a. 매듭점을 전후로 급격한 변화로 인하여 해당 점에서의 추정치를 신뢰할 수 없게 된다.

 

b. 매듭점 개수가 많다면 과적합이 일어날 수도 있다.

 

반응형

---

   4. 파이썬 구현

이제 개념을 알았으니 파이썬으로 구현해보아야 한다.

1) 클래스 정의

먼저 주어진 매듭점에 대해서 지시 함수를 만드는 데 필요한 StepFunction 클래스를 정의해주었다. 이 클래스는 get_value를 통하여 주어진 값이 구간에 들어오면 1 아니면 0을 반환한다.

 

class StepFunction():
    def __init__(self, lower_knot, upper_knot):
        assert lower_knot < upper_knot
        self.lower_knot = lower_knot
        self.upper_knot = upper_knot
        
    def get_value(self, x):
        lower_knot = self.lower_knot
        upper_knot = self.upper_knot
        lower_step = np.piecewise(x, [x<=lower_knot, x>lower_knot], [1, 0])
        upper_step = np.piecewise(x, [x<=upper_knot, x>upper_knot], [1, 0])
        return int(upper_step - lower_step)

 

다음으로 우리의 주인공 PiesewiseRegression 클래스를 구현했다. fit 메서드 설계 행렬 X와 반응 변수 벡터 y를 기본적으로 받고 constraint가 False이면 불연속 회귀 함수를, True면 연속 회귀 함수를 피팅한다. 그리고 (5)에서 행렬 $A$를 만들어주는 make_constraint_matrix, 매듭점으로 나누어지는 구간 개수만큼 지시 함수를 만드는 make_bases, 지시 함수와 다항 차수를 이용하여 실질적으로 회귀 모형 피팅에 사용되는 설계 행렬을 만드는 make_design_matrix 도 정의했다. 그리고 예측을 위한 predict 메서드도 구현하였다.

 

class PiecewiseRegression():    
    def __init__(self, knots, order=1):
        assert order<=3 and order>=0
        self.order = order
        self.knots = knots
        self.coef = None
        self.basis = None
#         self.model = None
        self.coef_ = None
        self.intercept_ = None
        self.fitted_values = None
        self.new_X = None
        self.constraint = None
        
    def fit(self, X, y, constraint=False):
        from sklearn.linear_model import LinearRegression
        from scipy.linalg import lapack
        self.constraint = constraint
        self.make_bases()
        
        new_X = self.make_design_matrix(X)
        self.new_X = new_X
        if not constraint:
            model = LinearRegression(fit_intercept=False).fit(new_X, y)
            self.coef_ = model.coef_
            self.fitted_values = model.predict(new_X)
        else:
            A = self.make_constraint_matrix()
            self.coef_ = lapack.dgglse(new_X, A, y, np.zeros(len(self.knots)))[3]
            self.fitted_values = new_X.dot(self.coef_)
        return self
    def make_constraint_matrix(self):
        row = []
        knots = self.knots
        order = self.order
        for k in range(len(knots)):
            knot = knots[k]
            temp_rows = np.zeros((order+1)*(len(knots)+1))
            temp_knot_poly = np.array([knot**o for o in range(order+1)])
            temp_rows[(order+1)*k:(order+1)*(k+2)] = np.concatenate([temp_knot_poly, -temp_knot_poly])
            row.append(temp_rows)

        return np.row_stack(row)
    def make_bases(self):
        knots = self.knots     
        basis = []
        num_knots = len(knots)
        for i in range(num_knots+1):
            if i == 0:
                basis.append(StepFunction(-np.infty, knots[i]))
            elif i == num_knots:
                basis.append(StepFunction(knots[i-1], np.infty))
            else:
                basis.append(StepFunction(knots[i-1], knots[i]))
                
        self.basis = basis
        
    def make_design_matrix(self, X):
        num_intervals = len(self.basis)
        order = self.order
        col = []
        for b in self.basis:
            for o in range(order+1):
                col.append(np.array([(x**(o))*b.get_value(x) for x in X.flatten()]))

        return np.column_stack(col)
    
    def predict(self, X):
        return np.array([self._predict(x) for x in X])
    
    def _predict(self, x):
        order = self.order
        x = x[0]
        temp_array = [b.get_value(x) for b in self.basis]
        temp_idx = temp_array.index(1)
        x_array = np.array([x**(o) for o in range(order+1)])
        target_coef = self.coef_[(order+1)*temp_idx:(order+1)*(temp_idx+1)]
        return x_array.dot(target_coef)

---

2) 실험

여기서는 실제 데이터가 아닌 시뮬레이션 데이터를 생성하여 앞에서 구현한 클래스가 잘 동작하는지 확인해본다. 시뮬레이션 데이터는 사인 함수에 노이즈를 추가하여 생성하였다. 그리고 매듭점을 3, 7로 설정하고 다항 차수를 1로 적합해보았다.

 

np.random.seed(100)
## 시뮬레이션 데이터
x = np.arange(0, 10, 0.1)
y = np.sin(x) + np.random.randn(len(x))
X = np.expand_dims(x, axis=1)

## PPR 적합
knots = [3, 7] ## 매듭점
pr = PiecewiseRegression(knots, order=1).fit(X, y) ## Piecewise Linear

## 시각화
fig = plt.figure(figsize=(10,6))
fig.set_facecolor('white')
ax = fig.add_subplot()
ax.scatter(x, y, color='k')
# ax.scatter(x, pr.fitted_values)
for i in range(len(knots)+1):
    add_vline = False
    if i == 0:
        idx = x<=knots[i]
        add_vline=True
    elif i == len(knots):
        idx = x>knots[len(knots)-1]
    else:
        idx = (x<=knots[i])&(x>knots[i-1])
        add_vline=True
    if add_vline:    
        ax.axvline(knots[i], linestyle='--', color='gray')
    ax.plot(x[idx], pr.predict(X)[idx], color='r')
plt.show()

 

보는 바와 같이 각 구간에서 독립적으로 선형 함수를 적합한 결과를 보여주었다. 이때 각 매듭점에서 불연속이다.

 

이번엔 각 매듭점에서 연속이 되도록 회귀 모형을 적합해보자. 위 코드에서 9번째 줄을 아래와 같이 바꿔준다.

 

pr = PiecewiseRegression(knots, order=1).fit(X, y, constraint=True)

 

 

보는 바와 같이 연속적으로 잘 피팅된 것을 알 수 있다.

 

이번엔 매듭점 개수만 알고 그 위치를 모른다고 했을 때 다항 차수를 2로 하여 최적 모형을 찾고자 한다. 여기서는 전진 선택법을 이용했으며 구간별 최소 샘플 수를 20으로 맞추었다.

 

from tqdm import tqdm

np.random.seed(100)
x = np.arange(0, 10, 0.1)
y = np.sin(x) + np.random.randn(len(x))
X = np.expand_dims(x, axis=1)

order = 2 ## 다항 차수
min_sample = 20 ## 구간별 최소 샘플 수
max_knot = 2 ## 찾아야할 매듭점 개수

## 최적 매듭점 찾기
num_knot = 0
final_knots = []
candidates = x[order:-order]
while num_knot < max_knot:
    obj = np.infty
    for i, knot in tqdm(enumerate(candidates), total=len(candidates)):
        knots = final_knots+[knot]
        knots = sorted(knots)
        ## 구간 사이의 최소 샘플 수를 만족 해야함.
        if len(knots) == 1:
            if not(np.sum(x <= knots[0]) >= min_sample and np.sum(x > knots[0]) >= min_sample):
                continue
        else:
            first = np.sum(x <= knots[0]) >= min_sample
            last = np.sum(x > knots[-1]) >= min_sample
            mid = [np.sum((x > knots[i])&(x<=knots[i+1])) >= min_sample for i in range(len(knots)-1)]
            if not(all(mid+[first]+[last])):
                continue
        pr = PiecewiseRegression(knots, order=order).fit(X, y, constraint=True)
        temp_val = np.linalg.norm(y-pr.predict(X))
        if temp_val < obj:
            obj = temp_val
            opt_knot_idx = i
            
    final_knots.append(candidates[opt_knot_idx])
    candidates = np.delete(candidates, opt_knot_idx)
    num_knot += 1
    
## 최종 매듭점
final_knots = sorted(final_knots)
pr = PiecewiseRegression(final_knots, order=order).fit(X, y, constraint=True)

fig = plt.figure(figsize=(10,6))
fig.set_facecolor('white')
ax = fig.add_subplot()
ax.scatter(x, y, color='k')
# ax.scatter(x, pr.fitted_values)
for i in range(len(final_knots)+1):
    add_vline = False
    if i == 0:
        idx = x<=final_knots[i]
        add_vline=True
    elif i == len(final_knots):
        idx = x>final_knots[len(final_knots)-1]
    else:
        idx = (x<=final_knots[i])&(x>final_knots[i-1])
        add_vline=True
    if add_vline:    
        ax.axvline(final_knots[i], linestyle='--', color='gray')
    ax.plot(x[idx], pr.predict(X)[idx], color='r')
plt.show()

 

코드를 돌려보면 아래와 같이 최종 모형(빨간 선)이 선택된 것을 알 수 있으며 사인곡선과 형태가 조금 다르지만 적합은 잘된 것 같다.

 

 

매듭점의 개수와 위치를 모르는 상황은 생략한다. 귀찮아서..

---

 

원래는 pyearth를 이용하여 내가 구현한 것과 비교해보고 싶었으나 설치가 안되어 못했다. 이번에는 머리로 대충 알고 있던 Piecewise Polynomial Regression에 대해서 공부하고 포스팅해보았다. 확실히 글로 쓰고 직접 구현해봐야 그 내용을 더 잘 알 수 있다. 매일 공부하면서 느끼는 거지만 세상에는 공부해야 할게 너무 많다. 그래도 하나씩 공부해서 포스팅하면 뭔가 쌓이는 느낌이 나에겐 너무 좋다~!!

 

- 참고자료 -

The Elements of Statistical Learning

---