[머신러닝] 파이썬 단순 선형 회귀분석 & 비선형 회귀분석 예제 소스코드

예전 글

[머신러닝] 파이썬 단순 선형 회귀분석 & 비선형 회귀분석 예제 소스코드

gyuuuna 2021. 8. 20. 23:02

※ 전체 소스코드는 글 최하단에 있습니다

※ Levenberg Method 및 Levenberg-Marquardt Method의 경우 damping factor를 계속 update해야 하는데 update하지 않도록 코드를 짰다.. 수정할 예정입니다!!!!

간단 이론 정리

선형 회귀는 주어진 데이터에 가장 잘 fit하는 선형 모델을 구하는 것이다. 모델이 주어진 데이터에 가장 잘 fit할 때 잔차(=추정값-실제값)의 합은 최소가 되며, 이때 likelihood/우도/가능도는 최대가 된다. likelihood/우도/가능도는 모델이 주어졌을 때 데이터가 관측될 확률을 말한다. 이를 응용하여, 우도/가능도를 최대화하여 선형 회귀를 하는 방법을 MLE(Maximum Likelihood Estimation)/최대 우도 추정이라고 한다. 위 이미지의 ①, ②는 우도를 최대하는 방법으로부터 모델이 fit할 때의 a, b 값을 구하는 방법을 나타내는 수식들이다. ②는 ①의 양변에 log를 씌운 결과다.

선형 회귀의 모수 추정 방법으로는 크게 두 가지가 있다. 바로 경사하강법(Gradient Descent)과 최소제곱법(Least Square)이다. 경사하강법은 iteration(반복)을 이용한 방법이다. iteration을 이용한 방법은 초깃값, step size을 잘 정하는 것이 중요하다. 위 그림에서 경사하강법 부분에 써놓은 식은 반복할 때마다 a, b의 값을 어떻게 갱신할 것인지를 나타낸다. 최소제곱법에는 크게 두 가지 방법이 있다. 해석적(analytic)인 방법은 오차 함수가 최솟값을 가질 때 이는 극값일 것이므로 오차 함수의 a, b 각각에 대한 편미분이 0일 것이라고 간주하고 그때의 a, b값을 구하는 것이다. 이때 역행렬을 이용하여 연립방정식을 풀며, inverse를 이용한다. 대수적(algebric)인 방법은 모든 실험데이터의 (x, y)쌍을 이용하여 연립방정식을 세우고, pseudo inverse를 이용해 연립방정식을 풀이하는 것이다.

Nonlinear regression에서는 circle fitting을 예시로 살펴보겠다. iteration을 이용해 파라미터를 업데이트해가며 fit하는 원의 중심=(a, b), 반지름=r에 해당하는 a, b, r을 구하도록 한다. 이때 iteration을 이용하는 4가지 방법 Gradient Descent, Gauss-Newton, Levenberg Method, Levenverg-Marquardt Method을 살펴보겠다.

실습 코딩

선형 회귀분석 코드 실습에 사용한 방법

경사 하강법(Gradient Descent)
최소제곱법(Least Square)

비선형 회귀분석(circle fitting) 코드 실습에 사용한 방법

경사 하강법(Gradient Descent)
가우스-뉴턴 방법(Gauss-Newton Method)
레벤버그 방법(Levenberg Method)
레벤버그-마쿼르트 방법(Levenberg-Marquardt Method)

1. 필요한 라이브러리 임포트

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from abc import abstractmethod, ABCMeta

numpy와 matplotlib.pyplot은 기본으로 임포트해주었다. 이따 살펴보겠지만, matplotlib.patches는 원 그리기를 위해, abc.abstractmethod와 abc.ABCMeta는 Iteration을 이용하는 방법들에 대하여 필수적으로 들어가야 할 함수를 지정하기 위해 임포트하였다.

2. Iteration 클래스

class Iteration(metaclass=ABCMeta):
    def __init__(self, x, y, a, b, r=1):
        self.x = x
        self.y = y
        self.a = a
        self.b = b
        self.r = r
    @abstractmethod
    def update(self):
        pass
    def plot_mse(self, str):
        plt.title(str)
        plt.xlabel('iteration')
        plt.ylabel('MSE')
        x = np.arange(1, len(self.mses)+1)
        y = np.array(self.mses)
        plt.plot(x, y)
        plt.show()

선형 회귀분석 코드 실습에 사용한 방법

경사 하강법(Gradient Descent)
최소제곱법(Least Square)

비선형 회귀분석(circle fitting) 코드 실습에 사용한 방법

경사 하강법(Gradient Descent)
가우스-뉴턴 방법(Gauss-Newton Method)
레벤버그 방법(Levenberg Method)
레벤버그-마쿼르트 방법(Levenberg-Marquardt Method)

빨간색으로 표시된 Iteration(반복)을 이용하는 방법들에 관한 클래스는 모두 Iteration 클래스를 상속하도록 하였다. Iteration 클래스에는 __init__메서드, update 메서드, plot_mse 메서드가 있다.

__init__ 메서드 : Iteration을 이용하는 방법으로 회귀분석을 하려면 x, y 변수에 들어가는 데이터와 a, b, r 등 각 회귀계수의 초깃값이 처음에 주어져야 한다. 선형 회귀분석의 경우 a, b 두 개의 회귀계수가 필요하고, Circle Fitting의 경우 a, b, r 세 개의 회귀계수가 필요하다. 따라서 a, b, r 세 개의 인스턴스 변수를 지정하되, 선형 회귀분석의 경우 r 값이 입력되지 않아도 메소드가 돌아가도록 r에는 default value를 넣어두었다.
update 메서드 : 모든 Iteration을 이용하는 방법에서 회귀계수의 값을 업데이트하는 과정이 필요하기 때문에 update method를 Abstract Method로 지정하였다.
plot_mse 메서드 : Iteration을 이용하는 방법의 경우 회귀계수의 값을 업데이트할 때마다 오차의 크기가 달라지므로, 오차가 감소하는 모습을 그래프를 통해 한눈에 확인할 수 있도록 Iteration(반복 횟수)-MSE(평균 제곱 오차) 그래프를 그려주는 메서드를 만들었다. 그래프의 제목은 str 인자로 받아 그래프 상단에 출력되도록 하였다.

3. LinearRegression 클래스

class LinearRegression():
    def __init__(self, x, y):
        self.x = x
        self.y = y
        
    def plot(self, str):
        plt.title(str)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.plot(self.x, self.y, 'o', color='red')
        plt.plot(self.x, self.x*self.a + self.b)
        plt.text(self.x.min(), self.a*self.x.min()+self.b, 'y={0}x+{1}'.format(round(self.a, 2), round(self.b, 2)))
        plt.text(self.x.max(), self.a*self.x.max()+self.b, 'MSE : {0}'.format(round(self.mse(), 2)))
        plt.show()
    
    def mse(self):
        val = np.sqrt(np.sum(np.square(self.x * self.a - self.y + self.b)))/2
        return val.sum()/len(self.x)

LinearRegression 클래스에는 선형 회귀를 하는 클래스에 공통으로 필요한 메서드들을 정의하여 상속할 수 있도록 하였다. __init__ 메서드, plot 메서드, mse 메서드가 들어간다.

__init__ 메서드 : 선형 회귀분석을 하려면 공통으로 필요한 것은 x, y 변수에 넣을 수 있는 초기 데이터 쌍들이다. 회귀계수 a, b는 Iteration을 이용한 방법의 경우에만 필요한 것이다. 따라서 LinearRegression 클래스에는 공통으로 필요한 x, y 값만 받아들일 수 있도록 했다.
plot 메서드 : 선형 회귀분석의 결과 만들어진 회귀선과 초기 데이터쌍들의 분포를 시각화하여 확인할 수 있도록 그래프를 그리는 메서드를 정의했다. plt.text를 이용하여 그래프에 회귀선에 해당하는 함수식이 나타나도록 했다.
mse 메서드 : MSE(평균제곱오차)을 구하는 공식은 모든 선형 회귀분석 메서드에 대하여 공통으로 적용할 수 있으므로 LinearRegression 클래스에 정의했다.

4. GradientDescent_line 클래스

(선형 회귀분석 ; 경사하강법 / Linear Regression ; Gradient Descent)

class GradientDescent_line(Iteration, LinearRegression):
    def __repr__(self):
        return "Gradient Descent"
    def setdata(self, step_size):
        self.step_size = step_size

    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())

            update_a = np.sum(-self.x * (self.y - self.x * self.a - self.b))
            update_b = np.sum(-1 * (self.y - self.x * self.a - self.b))
            self.a -= update_a*self.step_size
            self.b -= update_b*self.step_size

Iteration(반복)을 이용하여 선형 회귀분석을 하도록 하는 클래스이므로 Iteration과 LinearRegression 클래스를 상속받았다. 이때 회귀계수 a, b의 초깃값을 받아들이려면 Iteration의 __init__클래스를 이용해야 하므로 Iteration 클래스를 왼쪽에 적었다. __repr__, setdata, update 메서드를 정의했다.

__repr__메서드 : 클래스의 방법 이름을 출력하도록 하였다. 이는 plot 메서드를 이용하여 그래프를 그릴 때 제목으로 각 클래스의 __repr__ 메서드의 출력값을 넣기 위해 지정하도록 했다.
setdata 메서드 : Gradient Descent 방법으로 선형 회귀분석을 할 때 필요한 Step Size 값을 지정하는 메서드다.
update 메서드 : 회귀계수 a, b 값을 업데이트하는 메서드다. 자세한 내용은 위 이론정리 참고.

5. LS_Analytic, LS_Algebric 클래스

(최소제곱법 / Least Square)

class LS_Analytic(LinearRegression):
    def __repr__(self):
        return "Analytic Least Square"
    def solve(self):
        x_square = np.sum(self.x**2)
        x = np.sum(self.x)
        xy = np.sum(self.x*self.y)
        y = np.sum(self.y)

        A = np.array([[x_square, x], 
                      [x, 1]])
        B = np.array([xy, y])

        self.a, self.b = np.linalg.inv(A).dot(B)


class LS_Algebric(LinearRegression):
    def __repr__(self):
        return "Algebric Least Square"
    def solve(self):
        A = np.array([[i, 1] for i in self.x])
        B = np.array([i for i in self.y])
        self.a, self.b = np.dot(np.linalg.pinv(A), B)

LS_Analytic은 최소 제곱법의 해석적 방법, LS_Algebric은 최소 제곱법의 대수적 방법을 이용하는 클래스다. Iteration(반복)을 이용하지 않고 선형 회귀분석을 하도록 하는 클래스이므로 LinearRegression 클래스만을 상속받았다. __repr__, solve 메서드를 정의했다.

__repr__메서드 : 클래스의 방법 이름을 출력하도록 하였다. 이는 plot 메서드를 이용하여 그래프를 그릴 때 제목으로 각 클래스의 __repr__ 메서드의 출력값을 넣기 위해 지정하도록 했다.
solve 메서드 : 최종 회귀계수 a, b 값을 구하는 메서드다. 자세한 내용은 위 이론정리 참고.

6. 선형 회귀 실행

line_x = np.array([11, 12, 13, 14, 15, 16, 17, 18, 19, 20])
line_y = np.array([29, 33, 51, 40, 49, 50, 69, 70, 64, 89])

gdplot = GradientDescent_line(line_x, line_y, 8, 20)
gdplot.setdata(0.0008)
gdplot.update(50000)
gdplot.plot(gdplot)
gdplot.plot_mse(gdplot)

LSplot1 = LS_Analytic(line_x, line_y)
LSplot1.solve()
LSplot1.plot(LSplot1)

LSplot2 = LS_Algebric(line_x, line_y)
LSplot2.solve()
LSplot2.plot(LSplot2)

각 방법에 해당하는 클래스에 대하여 오브젝트를 생성한 후 메서드를 호출하여 선형 회귀를 실행한다.

GradientDescent_line

.setdata({Step Size}) : Step Size 값을 인자로 받아 오브젝트에 인스턴스 변수로 저장한다.
.update({Iteration}) : 실행할 Iteration 횟수를 인자로 받아 회귀계수를 업데이트할 때 그만큼 반복문이 돌도록 한다.
.plot({Object}) : 회귀선 및 입력 데이터쌍을 그래프로 시각화한다. 앞서 생성한 오브젝트를 인자로 받아, 해당 클래스에 __repr__ 메서드의 출력값으로 선언해 둔 문자열이 그래프의 제목이 되도록 한다.
.plot_mse({Object}) : Iteration(반복 횟수)-MSE(평균 제곱 오차)의 관계를 그래프로 시각화한다. 앞서 생성한 오브젝트를 인자로 받아, 해당 클래스에 __repr__ 메서드의 출력값으로 선언해 둔 문자열이 그래프의 제목이 되도록 한다.

LS_Analytic, LS_Algebric

.solve() : 회귀계수 a, b의 값을 구하여 저장하도록 한다.
.plot({Object}) : 회귀선 및 입력 데이터쌍을 그래프로 시각화한다. 앞서 생성한 오브젝트를 인자로 받아, 해당 클래스에 __repr__ 메서드의 출력값으로 선언해 둔 문자열이 그래프의 제목이 되도록 한다.

7. CircleRegression 클래스

class CircleRegression():
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def plot(self, str):
        plt.title(str)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.plot(self.x, self.y, 'o', color='red')
        c = patches.Circle((self.a, self.b), self.r, fc='w', ec='b')
        plt.gca().add_patch(c)
        plt.text(self.a-self.r, self.b, 'center : ({0}, {1})   radius : {2}   MSE={3}'.format(round(self.a, 2), round(self.b, 2), round(self.r, 2), round(self.mse(), 2)))
        plt.show()

    def mse(self):
        val = np.sqrt(np.sum(np.square(np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)**2) - self.r)))/2
        return val/len(self.x)

CircleRegression 클래스는 비선형 회귀 중에서도 Circle Fitting을 하는 클래스에 공통으로 필요한 메서드들을 정의하여 상속할 수 있도록 하였다. __init__ 메서드, plot 메서드, mse 메서드가 들어간다.

__init__ 메서드 : Circle Fitting을 하려면 공통으로 필요한 것은 x, y 변수에 넣을 수 있는 초기 데이터 쌍들이다. 회귀계수 a, b, r은 Iteration을 이용한 방법의 경우에만 필요한 것이다. 따라서 CircleRegression 클래스에는 공통으로 필요한 x, y 값만 받아들일 수 있도록 했다.
plot 메서드 : Circle Fitting의 결과 만들어진 회귀선과 초기 데이터쌍들의 분포를 시각화하여 확인할 수 있도록 그래프를 그리는 메서드를 정의했다. plt.text를 이용하여 회귀의 결과 나온 원의 중심 (a, b) 및 반지름 r의 값이 나타나도록 했다.
mse 메서드 : MSE(평균제곱오차)을 구하는 공식은 모든 Circle Fitting 메서드에 대하여 공통으로 적용할 수 있으므로 CircleRegression 클래스에 정의했다.

8. GradientDescent_circle 클래스

(비선형 회귀분석 Circle Fitting ; 경사하강법 / Nonlinear Regression Circle Fitting ; Gradient Descent)

class GradientDescent_circle(Iteration, CircleRegression):
    def __repr__(self):
        return "Gradient Descent"
    def setdata(self, step_size):
        self.step_size = step_size
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            grad = np.dot(J.T, F)
            update_val = grad*self.step_size
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]

Iteration(반복)을 이용하여 Circle Fitting을 하도록 하는 클래스이므로 Iteration과 CircleRegression 클래스를 상속받았다. 이때 회귀계수 a, b, r의 초깃값을 받아들이려면 Iteration의 __init__클래스를 이용해야 하므로 Iteration 클래스를 왼쪽에 적었다. __repr__, setdata, update 메서드를 정의했다.

__repr__메서드 : 클래스의 방법 이름을 출력하도록 하였다. 이는 plot 메서드를 이용하여 그래프를 그릴 때 제목으로 각 클래스의 __repr__ 메서드의 출력값을 넣기 위해 지정하도록 했다.
setdata 메서드 : Gradient Descent 방법으로 Circle Fitting을 할 때 필요한 Step Size 값을 지정하는 메서드다.
update 메서드 : 회귀계수 a, b, r 값을 업데이트하는 메서드다. 자세한 내용은 위 이론정리 참고.

9. GaussNewton 클래스

(선형 회귀분석 Circle Fitting ; 가우스-뉴턴 방법 / Nonlinear Regression Circle Fitting ; Gauss-Newton Method)

class GaussNewton(Iteration, CircleRegression):
    def __repr__(self):
        return "Gauss-Newton"
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            update_val = np.dot(np.linalg.pinv(J), F)
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]

__repr__메서드 : 클래스의 방법 이름을 출력하도록 하였다. 이는 plot 메서드를 이용하여 그래프를 그릴 때 제목으로 각 클래스의 __repr__ 메서드의 출력값을 넣기 위해 지정하도록 했다.
update 메서드 : 회귀계수 a, b, r 값을 업데이트하는 메서드다. 자세한 내용은 위 이론정리 참고.

10. Levenberg 클래스

(선형 회귀분석 Circle Fitting ; 레벤버그 방법 / Nonlinear Regression Circle Fitting ; Levenberg Method)

class Levenberg(Iteration, CircleRegression):
    def __repr__(self):
        return "Levenberg"
    def setdata(self, damping_factor):
        self.damping_factor = damping_factor
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            curve = np.linalg.inv(np.dot(J.T, J)+np.identity(n=3)*self.damping_factor)
            grad = np.dot(J.T, F)
            update_val = np.dot(curve, grad)
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]

__repr__메서드 : 클래스의 방법 이름을 출력하도록 하였다. 이는 plot 메서드를 이용하여 그래프를 그릴 때 제목으로 각 클래스의 __repr__ 메서드의 출력값을 넣기 위해 지정하도록 했다.
setdata 메서드 : Levenberg Method으로 Circle Fitting을 할 때 필요한 Damping Factor 값을 지정하는 메서드다.
update 메서드 : 회귀계수 a, b, r 값을 업데이트하는 메서드다. 자세한 내용은 위 이론정리 참고.

11. Levenberg_Marquardt 클래스

(선형 회귀분석 Circle Fitting ; 레벤버그-마쿼르트 방법 / Nonlinear Regression Circle Fitting ; Levenberg-Marquardt Method)

class Levenberg_Marquardt(Iteration, CircleRegression):
    def __repr__(self):
        return "Levenberg-Marquardt"
    def setdata(self, damping_factor):
        self.damping_factor = damping_factor
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            curve = np.linalg.inv(np.dot(J.T, J)+np.diag(np.dot(J.T, J))*self.damping_factor)
            grad = np.dot(J.T, F)
            update_val = np.dot(curve, grad)
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]

__repr__메서드 : 클래스의 방법 이름을 출력하도록 하였다. 이는 plot 메서드를 이용하여 그래프를 그릴 때 제목으로 각 클래스의 __repr__ 메서드의 출력값을 넣기 위해 지정하도록 했다.
setdata 메서드 : Levenberg-Marquardt Method으로 Circle Fitting을 할 때 필요한 Damping Factor 값을 지정하는 메서드다.
update 메서드 : 회귀계수 a, b, r 값을 업데이트하는 메서드다. 자세한 내용은 위 이론정리 참고.

12. 비선형 회귀 Circle Fitting 실행

circle_x = np.array([1.2, 3.1, 4.5, 2, 1.4, 0, -1.5, -1.9, -1.5])
circle_y = np.array([2.1, 2.5, 5, 8.3, 8.6, 8.2, 7.3, 5.4, 3.2])

GDcircle = GradientDescent_circle(circle_x, circle_y, 3, 3, 2)
GDcircle.setdata(0.01)
GDcircle.update(100)
GDcircle.plot(GDcircle)
GDcircle.plot_mse(GDcircle)

GNcircle = GaussNewton(circle_x, circle_y, 3, 3, 2)
GNcircle.update(5)
GNcircle.plot(GNcircle)
GNcircle.plot_mse(GNcircle)

Lcircle = Levenberg(circle_x, circle_y, 3, 3, 2)
Lcircle.setdata(0.1)
Lcircle.update(5)
Lcircle.plot(Lcircle)
Lcircle.plot_mse(Lcircle)

LMcircle = Levenberg_Marquardt(circle_x, circle_y, 3, 3, 2)
LMcircle.setdata(0.1)
LMcircle.update(5)
LMcircle.plot(LMcircle)
LMcircle.plot_mse(LMcircle)

각 방법에 해당하는 클래스에 대하여 오브젝트를 생성한 후 메서드를 호출하여 비선형 회귀 Circle Fitting을 실행한다.

GradientDescent_circle, GaussNewton, Levenberg, Levenberg_Marquardt

.setdata({Step Size}) : Step Size 값을 인자로 받아 오브젝트에 인스턴스 변수로 저장한다.
.update({Iteration}) : 실행할 Iteration 횟수를 인자로 받아 회귀계수를 업데이트할 때 그만큼 반복문이 돌도록 한다.
.plot({Object}) : 회귀 결과 생성된 원 및 입력 데이터쌍을 그래프로 시각화한다. 앞서 생성한 오브젝트를 인자로 받아, 해당 클래스에 __repr__ 메서드의 출력값으로 선언해 둔 문자열이 그래프의 제목이 되도록 한다.
.plot_mse({Object}) : Iteration(반복 횟수)-MSE(평균 제곱 오차)의 관계를 그래프로 시각화한다. 앞서 생성한 오브젝트를 인자로 받아, 해당 클래스에 __repr__ 메서드의 출력값으로 선언해 둔 문자열이 그래프의 제목이 되도록 한다.

코드 실행 결과

선형 회귀 (Linear Regression) 결과 - 순서대로

개인적으로 Algebric Least Square결과와 Gradient Descent 엄청 돌린 후 결과는 괜찮은데 Analytic Least Square 결과는 왜 저렇게 나오는지 모르겠다. 아시는 분 있으면 알려주세요 ...^^;

비선형 회귀 (Nonlinear Regression) Circle Fitting 결과 - 순서대로

Gradient Descent보다는 Gauss-Newton Method, Levenberg Method, Levenberg-Marquardt Method를 이용할 때 훨씬 더 빠르게 수렴하는 것을 확인할 수 있다. 더 좋은 비교를 위해 초깃값을 다르게 설정해보았다.

참값과 더 멀게 : a=10, b=10, c=5

전체 소스코드

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from abc import abstractmethod, ABCMeta

class Iteration(metaclass=ABCMeta):
    def __init__(self, x, y, a, b, r=1):
        self.x = x
        self.y = y
        self.a = a
        self.b = b
        self.r = r
    @abstractmethod
    def update(self):
        pass
    def plot_mse(self, str):
        plt.title(str)
        plt.xlabel('iteration')
        plt.ylabel('MSE')
        x = np.arange(1, len(self.mses)+1)
        y = np.array(self.mses)
        plt.plot(x, y)
        plt.show()
    


class LinearRegression():
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def plot(self, str):
        plt.title(str)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.plot(self.x, self.y, 'o', color='red')
        plt.plot(self.x, self.x*self.a + self.b)
        plt.text(self.x.min(), self.a*self.x.min()+self.b, 'y={0}x+{1}'.format(round(self.a, 2), round(self.b, 2)))
        plt.text(self.x.max(), self.a*self.x.max()+self.b, 'MSE : {0}'.format(round(self.mse(), 2)))
        plt.show()
    
    def mse(self):
        val = np.sqrt(np.sum(np.square(self.x * self.a - self.y + self.b)))/2
        return val.sum()/len(self.x)



class GradientDescent_line(Iteration, LinearRegression):
    def __repr__(self):
        return "Gradient Descent"
    def setdata(self, step_size):
        self.step_size = step_size

    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())

            update_a = np.sum(-self.x * (self.y - self.x * self.a - self.b))
            update_b = np.sum(-1 * (self.y - self.x * self.a - self.b))
            self.a -= update_a*self.step_size
            self.b -= update_b*self.step_size


class LS_Analytic(LinearRegression):
    def __repr__(self):
        return "Analytic Least Square"
    def solve(self):
        x_square = np.sum(self.x**2)
        x = np.sum(self.x)
        xy = np.sum(self.x*self.y)
        y = np.sum(self.y)

        A = np.array([[x_square, x], 
                      [x, 1]])
        B = np.array([xy, y])

        self.a, self.b = np.linalg.inv(A).dot(B)


class LS_Algebric(LinearRegression):
    def __repr__(self):
        return "Algebric Least Square"
    def solve(self):
        A = np.array([[i, 1] for i in self.x])
        B = np.array([i for i in self.y])
        self.a, self.b = np.dot(np.linalg.pinv(A), B)



#Circle Fitting



class CircleRegression():
    def __init__(self, x, y):
        self.x = x
        self.y = y
        
    def plot(self, str):
        plt.title(str)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.plot(self.x, self.y, 'o', color='red')
        c = patches.Circle((self.a, self.b), self.r, fc='w', ec='b')
        plt.gca().add_patch(c)
        plt.text(self.a-self.r, self.b, 'center : ({0}, {1})   radius : {2}   MSE={3}'.format(round(self.a, 2), round(self.b, 2), round(self.r, 2), round(self.mse(), 2)))
        plt.show()

    def mse(self):
        val = np.sqrt(np.sum(np.square(np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)**2) - self.r)))/2
        return val/len(self.x)


class GradientDescent_circle(Iteration, CircleRegression):
    def __repr__(self):
        return "Gradient Descent"
    def setdata(self, step_size):
        self.step_size = step_size
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            grad = np.dot(J.T, F)
            update_val = grad*self.step_size
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]


class GaussNewton(Iteration, CircleRegression):
    def __repr__(self):
        return "Gauss-Newton"
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            update_val = np.dot(np.linalg.pinv(J), F)
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]


class Levenberg(Iteration, CircleRegression):
    def __repr__(self):
        return "Levenberg"
    def setdata(self, damping_factor):
        self.damping_factor = damping_factor
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            curve = np.linalg.inv(np.dot(J.T, J)+np.identity(n=3)*self.damping_factor)
            grad = np.dot(J.T, F)
            update_val = np.dot(curve, grad)
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]

class Levenberg_Marquardt(Iteration, CircleRegression):
    def __repr__(self):
        return "Levenberg-Marquardt"
    def setdata(self, damping_factor):
        self.damping_factor = damping_factor
    def update(self, iter):
        self.mses = []
        for i in range(iter):
            self.mses.append(self.mse())
            J_col0 = (self.a - self.x)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col1 = (self.b - self.y)/np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b))
            J_col2 = -np.ones_like(self.x)
            J = np.column_stack((J_col0, J_col1, J_col2))
            F = np.sqrt(np.square(self.x-self.a) + np.square(self.y-self.b)) - self.r
            F.resize(len(F), 1)
            curve = np.linalg.inv(np.dot(J.T, J)+np.diag(np.dot(J.T, J))*self.damping_factor)
            grad = np.dot(J.T, F)
            update_val = np.dot(curve, grad)
            self.a -= update_val[0][0]
            self.b -= update_val[1][0]
            self.r -= update_val[2][0]




circle_x = np.array([1.2, 3.1, 4.5, 2, 1.4, 0, -1.5, -1.9, -1.5])
circle_y = np.array([2.1, 2.5, 5, 8.3, 8.6, 8.2, 7.3, 5.4, 3.2])

GDcircle = GradientDescent_circle(circle_x, circle_y, 3, 3, 2)
GDcircle.setdata(0.01)
GDcircle.update(100)
GDcircle.plot(GDcircle)
GDcircle.plot_mse(GDcircle)

GNcircle = GaussNewton(circle_x, circle_y, 3, 3, 2)
GNcircle.update(5)
GNcircle.plot(GNcircle)
GNcircle.plot_mse(GNcircle)

Lcircle = Levenberg(circle_x, circle_y, 3, 3, 2)
Lcircle.setdata(0.1)
Lcircle.update(5)
Lcircle.plot(Lcircle)
Lcircle.plot_mse(Lcircle)

LMcircle = Levenberg_Marquardt(circle_x, circle_y, 3, 3, 2)
LMcircle.setdata(0.1)
LMcircle.update(5)
LMcircle.plot(LMcircle)
LMcircle.plot_mse(LMcircle)



line_x = np.array([11, 12, 13, 14, 15, 16, 17, 18, 19, 20])
line_y = np.array([29, 33, 51, 40, 49, 50, 69, 70, 64, 89])

gdplot = GradientDescent_line(line_x, line_y, 8, 20)
gdplot.setdata(0.0008)
gdplot.update(50000)
gdplot.plot(gdplot)
gdplot.plot_mse(gdplot)

LSplot1 = LS_Analytic(line_x, line_y)
LSplot1.solve()
LSplot1.plot(LSplot1)

LSplot2 = LS_Algebric(line_x, line_y)
LSplot2.solve()
LSplot2.plot(LSplot2)

저작자표시

현재글[머신러닝] 파이썬 단순 선형 회귀분석 & 비선형 회귀분석 예제 소스코드

빅테크, HTML문서, 백엔드, 장고, 코딩, 웹개발, 프론트엔드, HTML, NLP, 워드임베딩, github, IT, 파이썬, 웹, 프로그래밍, 네이버, Python, git, 카카오, 모듈,

gyuuuna space