R에서 CDA를 통해 회귀식 구현하기

앞선 글에서 좌표하강알고리즘(Coordinate Descent Algorithm, CDA)에 대해서 알아보았다. 

https://mn0317.tistory.com/10

좌표하강알고리즘(Coordinate Descent Algorithm, CDA)

이번에는 구체적으로 R에서 CDA를 통해 얻은 식과 lm함수를 통해 얻은 회귀계수를 비교해봄으로써, CDA알고리즘을 구현해보도록 하겠다. 

 

0. 부분잔차의 규칙

 앞선 글에서 우리의 부분잔차는 다음과 같이 정의된다. 

다음의 그림은 그저 r번째의 부분잔차와 r-1번째의 부분잔차 사이의 규칙을 나타낸 것이다.

 부분잔차에 해당하는 식을 r번째와 r-1번째를 빼보면 위와 같은 식이 발생한다. 

즉, 한글로 쓰면 다음과 같다. 

부분잔차= 직전부분잔차-직전 갱신된 이전좌표 계수의 적합 + 이전반복에 계산된 해당계수의 적합

이제 부터 코딩으로 직접 실행해보도록 하겠다.

 

1. 데이터 생성 및 lm 함수 사용

코드를 설명하자면 우선, x값으로 -5부터 5까지 300개의 수들을 생성해서 저장하였다. 그 후 2x+3의 일차식의 값을 f로 저장하였다. 그리고 오차를 발생시키기 위해, 300개의 오차들을 발생시켰고 이를 e로 저장하였다. y의 값은 f와 오차의 합으로 저장하였다. 

x와 y의 산점도와 실제 f와 x와의 관계의 선을 나타내면 다음과 같이 나타난다.

 

x와 y의 값의 산포도. 빨간선은 x와 f의 관계이다.

 

이를 lm 함수를 이용해서 나타내면 다음과 같다. 

 

lm함수를 이용해 추정 회귀계수들을 구해보았다.

 

약간의 오차가 있긴 있으나, y절편의 값은 3에, 베타 2의 값은 2에 근접한 것을 알 수 있다.

 

2. CDA 함수 구현

 다음은 CDA 함수를 직접 구현한 내용이다. 

CDA 알고리즘의 최종 함수이다.

코드 한줄 한줄 설명해보도록 하겠다. 

 

먼저 함수의 변수들로 기저함수인 B와 y의 값인 y, stopping rule의 기준을 설정해주는 epsilon, 최대 반복수를 나타내는 max_iter을 설정하였다. 

여기서 J는 베타 값들의 개수이다. 앞선 글에서 살펴보았듯, 함수 f는 기저함수 B와 회귀계수 beta들의 곱의 선형결합이므로, 행렬의 곱을 하기 위해선 회귀계수의 개수는 기저함수행렬 B의 열의 개수와 같아야 한다. 따라서 회귀계수의 총 개수인 J를 기저함수 행렬 B의 열의 개수로 설정해 두었다. 

 

이제 다음으로는 초기값 설정이다. 

베타의 초기값은 영벡터로 설정해 두었다. 그리고 초기 잔차의 값은 (y-적합값)에서 초기 베타들의 값이 0이므로 그냥 y에 해당하는 것을 알 수 있다.

밑의 잔차 이전의 값은 초반에는 값이 갱신이 되지 않았으므로 그냥 잔차제곱의 합의 절반으로 설정해 두었다. 

식으로 살펴보면 위의 식에서 beta j의값들이 0이므로 그냥 y의 값들의 제곱의 합인 것을 알 수 있다.

 

이제 핵심인 for 구문이다. CDA는 결국 반복을 통해 최소해를 찾아가는 과정이므로 for 반복문을 사용하였다. 

첫번째 반복문에서 먼저 반복 횟수인 r은 1부터 설정해둔 값인 max_iter만큼 반복한다. 

그리고 그 안에서 beta 1의 값부터 beta J의 값들이 각각 갱신됨으로 이를 for 구문으로 또 설정하였다. 

이제 그 안에서 앞서 살펴보았던 부분잔차를 partial_residual로 설정해준다. 

그리고 이전 글에서 구했던 beta j의 값을 구하는 식을 설정해둔다.

거기서 구해진 베타 j의 값을 집어넣어 다시 잔차를 갱신한다. 

이를 beta J의 값이 나올때까지 갱신하면 반복 한 번이 진행된 것이다. 이렇게 갱신된 새로운 잔차의 값의 제곱합의 절반으로 RSS_new를 생성한다.

그리고 for 구문이 이 새로운 잔차제곱합 RSS_new와 RSS_old의 값의 차이가 우리가 설정해둔 epslion보다 작을 경우 첫번째 for 구문 반복을 멈추도록 설정해둔다. 

끝으로, 갱신된 RSS_new의 값을 다시 RSS_old로 넣어줌으로써, 다음 반복에서 발생하는 새로운 RSS_new의 값과 비교 가능하도록 만들어준다. 

 

3. lm 함수와 비교

 이제 만든 함수 CDA_lm에 똑같은 값을 넣어 R패키지에 있는 lm함수와 비교해보도록 하겠다. 

일단 기저함수 B를 다음과 같이 설정해준다. 

 

기저함수는 우리가 추정하는 회귀식이 단순선형회귀임으로, B1은 1, B2는 x가 될 것이다. 

그래야 베타1과 베타2와 해당 기저함수들의 각각의 곱의 합으로 함수가 정의될 수 있기 때문이다. 

 

설정한 기저함수 B와 앞서 저장한 y의 값을 CDA_lm함수에 넣어보았다. (여기서 쨰는 오타입니다.)

 설정한 기저함수 행렬 B와 앞서 저장한 y의 값을 집어 넣어 결과를 보니, 회귀계수의 값들이 나오는 것을 볼 수 있다. 

2번째 바퀴수만에 해당 값들이 나온 것을 알 수 있다. 

 

이제 lm함수의 결과값과 비교해보도록 하겠다.

 

cda_lm의 함수값이다.

앞서 설정한 lm함수의 함수값이다.

비교해보면 알겠지만, 두 회귀계수가 동일하게 나오는 것을 알 수 있다. 

CDA회귀선이다.

그림을 그려봐도 아까 lm을 통해 구한 회귀직선과 동일하게 나오는 것을 알 수 있다. 즉, 좌표하강알고리즘이 제대로 구현되었다는 것을 알 수 있다.

 

오늘은 CDA를 R로 직접구현하여 단순선형회귀식의 회귀계수들을 직접 구하는 과정까지 코드를 짜보았다. 다음 글에서는 LASSO 회귀분석을 CDA로 구현해보려 한다. 

 

LASSO 회귀분석에 대한 내용은 밑의 글을 참고하면 좋을 듯하다.

https://mn0317.tistory.com/6

Lasso 회귀분석