최근에 개인적인 필요에 의해 nonlinear regression 을 접하고 있는데 아래 포스팅을 보면서 개념을 잡아가고 있습니다.

http://darkpgmr.tistory.com/58

nonlinear regression 알고리즘 중 유명한 것이 Levenberg-Marquardt 라고 하는데 Levenberg-Marquardt 알고리즘은 gradient descent 와 Gauss–Newton 방법이 합쳐진 형태라고 하더군요. 그래서 일단 gradient descent 는 아니까 Gauss–Newton 방법이 무엇인지에 대해 공부를 하기 위해 검색을 하던 중 위 블로그를 발견하게 되었고 백문이 불여일타라고 일단 코딩해 봐야 이해가 되니까 위 포스팅에서 맨 아래 있던 예제를 MATLAB 으로 돌려 봤습니다.

예제는 원 방정식에 대해 가우스-뉴턴 방법을 적용해서 파라미터를 찾는 것입니다.

원의 중심과 반지름을 a, b, r 이라고 하고 이를 벡터로 표현하면 x=[a; b; r] 입니다.

원의 방정식(F)는 아시는 바와 같이 아래 수식과 같고~ 아래 수식에서 대문자 X 와 소문자 x 가 다르다는것을 염두에 두시기 바랍니다…. 다른 문자로 쓸걸 그랬나…. ??

위 수식에 대하여 변수 a, b, r 에 대한 자코비안(J)는 아래와 같습니다.

Gauss–Newton 업데이트 수식은 아래와 같습니다. 아래 수식에서 ‘ 는 Transpose 를 의미 하고 inv() 함수는 matrix inverse 함수입니다.

x = x - inv(J(x)' * J(x))*J(x)'*F(x)

MATLAB 에서 inv 함수를 사용하는것 보다는 역슬래쉬(\) 를 사용하는게 더 추천하는 방식이기 때문에 위 수식은 아래와 같이 구현 했습니다.

x = x - (J(x)' * J(x))\J(x)'*F(x)


반지름이 3 이고 중심이 (x,y) = (1,3) 인 원의 방정식을 하나 만들었고 제가 작성한 MATLAB 코드는 아래와 같습니다.

NOP = 100;

center = [1,3];

radius = 3;

style='b-';

THETA=linspace(0,2*pi,NOP);

RHO=ones(1,NOP)*radius;

[X,Y] = pol2cart(THETA,RHO);

X=X+center(1);

Y=Y+center(2);

H=plot(X,Y,style);

grid on

axis square;

axis equal;

다음으로 해당 원의 방정식에 대해 Gauss–Newton 방정식을 적용하는 코드는 아래와 같습니다.

XX=X(:);

YY=Y(:);

x = [0;1;15];

% 원의 방정식에 대한 가우스-뉴턴 방법

for n=1:60

Y=YY(n:(n+10));

X=XX(n:(n+10));

F = sqrt((X-x(1)).^2 + (Y-x(2)).^2)-x(3);

J = [(x(1)-X)./sqrt((X-x(1)).^2 + (Y-x(2)).^2) (x(2)-Y)./sqrt((X-x(1)).^2 + (Y-x(2)).^2) -ones(size(X))];

x = x - (J' * J)\J'*F;

end

x

위 코드에서는 초기값은 0, 1, 15 로 했고 돌려보니 아래와 같이 원의 중심과 반지름을 정확히 추정함을 확인 할 수 있었습니다.

x =

1.0000

3.0000

3.0000

그런데 다들 아시겠지만 Gauss–Newton 방식은 초기값에 따라 발산의 우려가 있다고 합니다. 저도 실험을 해보니 초기값을 다른 값들로 해서 실험을 많이 하다 보니(ex x = [0;8;15];) 로 한 경우에는 추정값이 발산하더군요.



MATLAB 에서 txt 파일을 출력하는 경우 다음 과정에 따라 처리 한다.

 

  • fopen() 함수를 이용하여 파일을 쓰기 모드('w')로 연다.

 

  • fprintf() 함수를 이용하여 파일에 내용을 쓴다.

   

  • fclose() 를 이용하여 file handle 을 닫아 준다.

  



 

다음과 같은 데이터를 이용하여 위 과정 대로 test_file.txt 파일에 저장해 보자.

 

 

저장하는 과정은 다음 코드와 같다.

 

 

 

저장된 파일 test_file.txt 을 열어서 확인해 보면 다음과 같다.

 

7 8 

3 7 

7 1 

1 5 


위 결과를 보면 정상적으로 저장이 안 된 것을 확인 할 수 있다.



 

이는 fprintf 함수에서 %d %d 를 써서 10진수 숫자를 저장할 때

 

각 열 벡터 별로 즉 MATLAB 이 일반적으로 처리하는 행 순서대로 저장하기 때문이다.

 

따라서 matrix 값을 정확히 저장하기 위해서는 transpose 를 취해준 후 저장 해야 한다.

 

fprintf() 문장의 올바른 코드는 다음과 같다.

 

fprintf(file_h,'%d %d \n', x.'); % 10진수 이므로 %d 사용

 


텍스트 파일을 열어서 결과를 확인하면 다음과 같이 정상적으로 저장 된 것을 확인 할 수 있다.

 

7 7 

8 1 

3 1 

7 5  

    


'programming language > MATLAB' 카테고리의 다른 글

MATLAB nargin, nargout  (6) 2011.04.24
MATLAB laplace, inverse laplace transform  (0) 2011.04.24
MATLAB 변수를 저장 하자~ mat file  (0) 2011.04.21
MATLAB NaN  (0) 2011.04.20
MATLAB 파일 출력  (0) 2011.04.19
MATLAB 파일 읽기 importdata()  (0) 2011.04.18
MATLAB binomial r.v. generation  (4) 2011.04.17
MATLAB cell class  (6) 2011.04.11
MATLAB taylor, Maclaurin serise  (0) 2011.04.10

MATLAB 의 rot90() 함수는 matrix 를 시계 반대방향으로 돌리는 기능을 하는 함수입니다.

 

x=magic(4)

x =

16 2 3 13

5 11 10 8

9 7 6 12

4 14 15 1

 

위 x 에 대하여 다음과 같이 하면 시계 방대방향으로 90 도 돌리게 됩니다.

 

x1=rot90(x)

x1 =

13 8 12 1

3 10 6 15

2 11 7 14

16 5 9 4

 

위 x 에 90 도씩 3번 시계 방대방향으로 돌리고 싶다면 다음과 같이 하면 됩니다.

 

x2=rot90(x,3)

x2 =

4 9 5 16

14 7 11 2

15 6 10 3

1 12 8 13

 

 

다음으로 설명한 명령어는 triu(), tril() 함수 입니다.

 

triu() 함수는 upper triangular 부분을 구하는 함수이며 tril() 함수는 lower triangular 부분을 구하는 함수입니다.

 

triu() 함수를 실행하면 다음과 같이 됩니다. diagonal 을 포함하여 upper triangular 부분만 남아 있고 나머지는 다 0 이 되게 됩니다.

 

x_triu=triu(x)

x_triu =

16 2 3 13

0 11 10 8

0 0 6 12

0 0 0 1

 

 

triu() 함수에 diagonal 로 부터의 offset 을 줘서 다음과 같이 diagonal 부분을 제외하고 선택도 가능합니다.

 

x_triu2=triu(x,1)

x_triu2 =

0 2 3 13

0 0 10 8

0 0 0 12

0 0 0 0

 

triu() 함수에 – 의 옵셋을 주게 되면 diagonal 아래 방향으로 옵셋을 주게 됩니다.

 

x_triu3=triu(x,-1)

x_triu3 =

16 2 3 13

5 11 10 8

0 7 6 12

0 0 15 1

 

tril() 함수는 다음과 같이 사용합니다.

 

x_tril= tril(x)

x_tril =

16 0 0 0

5 11 0 0

9 7 6 0

4 14 15 1

 

triu() 함수와 사용방법은 비슷하지만 옵셋을 줄 때는 조심하셔야 합니다.

 

x_tril1= tril(x,1)

x_tril1 =

16 2 0 0

5 11 10 0

9 7 6 12

4 14 15 1

 

triu() 함수에서 + 옵셋을 줬을 때는 triangular 의 크기가 줄어들었는데 tril() 함수에서는 + 옵셋을 줬을 때 triangular 의 크기가 커지는 것을 알 수 있습니다.

 

tril(), triu() 함수의 옵셋은 무조건 diagonal 을 기준으로 + 면 위쪽으로 – 면 아래 쪽으로 triangle 의 크기가 변한다고 알고 계시면 됩니다.

 

tril() 함수에서 -1로 옵셋을 주니깐 삼각형이 줄어들게 됩니다.

 

x_tril2= tril(x,-1)

x_tril2 =

0 0 0 0

5 0 0 0

9 7 0 0

4 14 15 0

 

 

MATLAB 에서 '(쉼표) 는 hermitian transpose 를 의미 합니다.

 

hermitian transpose 란 complex conjugate transpose 를 하는 것입니다.

 

다음과 같은 matrix 에 대해 hermitian transpose 를 취해 보죠

 

x=magic(3)+i*magic(3)

x =

8.0000 + 8.0000i 1.0000 + 1.0000i 6.0000 + 6.0000i

3.0000 + 3.0000i 5.0000 + 5.0000i 7.0000 + 7.0000i

4.0000 + 4.0000i 9.0000 + 9.0000i 2.0000 + 2.0000i

 

위 결과를 보면 imaginary 의 부호가 바뀌었고 transpose 된 것을 확인 할 수 있습니다.

 

x'

ans =

8.0000 - 8.0000i 3.0000 - 3.0000i 4.0000 - 4.0000i

1.0000 - 1.0000i 5.0000 - 5.0000i 9.0000 - 9.0000i

6.0000 - 6.0000i 7.0000 - 7.0000i 2.0000 - 2.0000i

 

 

complex conjugate 는 안 하고 transpose 만 할 때는 .' (점 - 쉼표)를 이용해야 합니다.

 

x.'

ans =

8.0000 + 8.0000i 3.0000 + 3.0000i 4.0000 + 4.0000i

1.0000 + 1.0000i 5.0000 + 5.0000i 9.0000 + 9.0000i

6.0000 + 6.0000i 7.0000 + 7.0000i 2.0000 + 2.0000i

 

 

MATLAB 에서 대소문자 변환 관련한 함수로는 upper(), lower() 함수가 있습니다.

 

upper()함수는 대문자로 변환해 주며

 

X=upper('I am a man')

X =

I AM A MAN

 

lower() 함수는 소문자로 변환 해주는 함수 입니다.

 

lower(X)

ans =

i am a man

  


+ Recent posts