파이썬을 이용하여 허블 - 르메트르의 법칙 그래프 그리기 - curve_fit 을 이용한 추세선 그리기를 이용하여

<허블 - 르메트르의 법칙 (reference : space.fm)>

허블 - 르베트르의 법칙은 지구과학 1 교과서에서 매우 중요하게 다루어지는 내용일 뿐 아니라, 우주 팽창과 관련하여 과학적으로도 매우 중요하고, 천문학사에서도 굉장히 중요한 내용으로 다루어 집니다. 이번 포스팅에서는 허블 - 르메트르의 법칙을 파이썬을 이용해 그래프를 그리는 방법을 알아보고자 합니다. 그동안 파이썬 관련 포스팅 중 가장 어려운 내용이 될 것 같아 조금 두렵기도 하지만, 최대한 쉽게 풀어 설명해 보겠습니다.

 

  1. 자료 수집

  자료는 SIMBAD에서 추출하였습니다. SIMBAD는 천문학 관측 및 연구 결과를 수록한 데이터 베이스라고 생각하면 되는데, 여기에 방대한 천문학적 자료가 수록되어 있습니다. 외부은하의 적색편이나 거리 정보도 수록되어 있을 것으로 생각하고 열심히 찾았는데, 겨우겨우 찾은것이 Tully라는 학자가 1988년에 연구한 데이터 뿐이었습니다. 분명 다른 자료가 있을 것으로 생각되나 이것밖에 못 찾았으니 우선 이걸로 연습을 해 보겠습니다. (다만, 오래전의 관측 결과라 허블 상수가 지금의 연구 결과와 다소 상이합니다.)

 

관심있는 분들은 SIMBAD에서 직접 찾아보셔도 되고, 우선 제가 여기에 실습용 csv 파일을 올리겠습니다. 아래 데이터를 다운받아 사용하시면 됩니다. 아무런 바이러스도 심지 않았으니 안심하셔도 되요.ㅎㅎ

tullydata.csv

0.00MB

SIMBAD에서 어떤 데이터를 다운하느냐에 따라 얻을 수 있는 데이터의 종류가 많이 달라지는데, 저의 경우 은하의 거리, 후퇴속도 뿐 아니라 적경, 적위, 후퇴속도 에러 등등 불필요한 정보도 혹시나 하여 함께 받았습니다. 파일을 열어보면

요런 아이들이 있는데, header 정보의 이름은 작업하기 편하게 제가 인위적으로 조금 바꿨습니다. SIMBAD에서 원본을 다운받으시면 아마 저런 header 이름은 아닐 것입니다. 적당히 편집해야 하므로, 그냥 제가 올려놓은 저 파일을 쓰는게 나을 것입니다.

** 자료 출처(reference : SIMBAD, VII/145/catalog, Nearby Galaxies Catalogue (NBG) (Tully 1988))

반응형

  2. 파이썬 코드 입력

이제 코드를 짜고 그래프를 불러올 차례입니다. 코드는 일반적으로 많이 이용되는 그래프를 그리는 코드와, 라이브러리 pandas를 이용한 csv 파일 불러오는 코드가 대부분인데

 

1. 보통 은하의 후퇴속도와 거리 그래프에서, 관측한 은하의 후퇴속도와 거리는 점으로 찍어 표현하고,

2. 여러개의 점을 추세선을 그려(fitting) 평균한 값을 사용합니다.

3. 그리고, 여기서는 오차 막대(error bar)까지 이용하여, 후퇴속도의 관측 오차까지 표현해 보고자 합니다.

 

1,2번은 scipy라는 라이브러리에 있는 curve_fit이라는 기능을 이용할 것이며, error bar는 matplotlib에 있는 기능을 이용해 볼 것입니다.

 

우선 전체 코드를 보겠습니다.

 

import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

def objective(x, a,b):
    return a*x

data=pd.read_csv("C:\\111\\tullydata.csv")
x=data['distance']
y=data['velocity']
plt.figure(figsize=(9,7))

plt.errorbar(x, y, yerr=data['err'], fmt='.k', capsize=1.5, elinewidth=0.5)

popt,pcov = curve_fit(objective, x, y, method='trf', bounds=([78,-0.000001],[79,0.000001]))

plt.plot(x,objective(x,*popt),color='red',linestyle='-', linewidth=0.5)

c,d=popt
plt.grid(ls='-.')
plt.xlim(0)
plt.ylim(0)
plt.title("Hubble's law",size=18)
plt.xlabel('distance(Mpc)', size=15)
plt.ylabel('velocity(km/s)', size=15)

f=str(round(c,1))

print('V='+f+'*d')
plt.text(22, 1500, 'V='+f+'*d', size=15)
plt.savefig('C:\\111\\hubble.jpg', dpi=300)
plt.show()

이 코드를 돌리면 아래처럼 나옵니다.

  3. 파이썬 코드 해석

1) 중요한 코드에만 색을 칠했습니다. 우선

 

def objective(x, a,b):
    return a*x 

에서는 추세선을 그려 fitting할 방법을 결정하는 겁니다. 허블상수는 1차함수 그래프로 표현되기 때문에, def를 이용하여 추세선이 1차함수로 그려지게 하였습니다. y 절편은 필요가 없어 매게변수를 x,a만 해도 될 것 같은데, 2개로만 하니 오류가 뜨더군요. 대신, fitting할 함수에는 y절편을 표시하지 않아 별 문제가 없습니다.

 

2) 에러바(error bar) 그리기

 

plt.errorbar(x, y, yerr=data['err'], fmt='.k', capsize=1.5, elinewidth=0.5)

 

이번 포스팅에서 첫 번째로 중요한 코드로, 에러바를 그리는 코드입니다. 에러바를 그릴려면, 3개의 매게변수가 필요한데, x, y는 그냥 x, y축에 해당하는 값을, yerr는 y축 방향의 에러 값을 나타냅니다. x축 방향의 에러값은 xerr로 하면 됩니다. 

 

코드를 자세히 보시면, x는 x=data['distance']이고, y는 y=data['velocity'] 입니다. 각각 거리와 후퇴속도를 의미하겠죠?

 

fmt는 그래프를 그릴 형식입니다. 빈 칸으로 두면, 점으로 찍지 않고 선으로 모두 그려버려 그래프가 이상해 집니다. fmt='.k'에서 .은 점으로 찍어 그래프를 그려라 하는 것을 의미하고, k는 점의 색을 의미합니다.

 

capsize는 에러바를 자세히 보면 에러바의 양 끝단에 짧은 가로줄이 있습니다. 이것의 크기를 말합니다. elinewidth는 에러바의 두께입니다. 기본값으로 두어도 좋지만, 그래프 전체 크기에 비해 이쁘게 그래프가 잘 나오지 않을 수 있으니 지정해 주는것이 좋습니다.

 

3) 추세선(fitting line) 기본 값 지정하기

 

popt,pcov = curve_fit(objective, x, y, method='trf', bounds=([78,-0.000001],[79,0.000001]))

 

이번 포스팅에서 두 번째로 중요한 코드입니다. scipy에 있는 curve_fit입니다. 앞서 def에서 정의한 a*x를 위 데이터에 맞춤(fitting)하여, 알맞은 a값을 찾아내는 것입니다. objective는 def에서 정의한 함수이고, 이 함수에 맞추어 x,y값에 적합한 fitting을 하라는 의미입니다. 그래서 curve_fit 안에 기본 매개변수가 3개가 들어간 것입니다.

 

매개변수 이후 나오는 method는 fitting의 방법을 지정해 주는 것인데, 크게 lm, trf, dogbox라는 3가지 방법이 있습니다. 각각에 대한 설명까지 다루는 것은 글의 목적을 한참 벗어나니 생략하고, method를 따로 지정해주지 않으면 기본적으로는 lm으로 피팅합니다. 그런데 뒤에 bounds값이 주어진 경우 기본인 lm으로는 피팅이 되지 않습니다. bounds값이 있는 경우에는 trf나 dogbox로만 피팅이 가능합니다. 그런데 trf로 하든 dogbox로 하든 피팅 결과는 크게 다르지 않습니다.

 

다음으로 bounds를 보겠습니다. bounds에서는 a*x+b라는 피팅 함수에서 a,b의 최소값과 최대값을 지정해 주는 것입니다. 위에서 함수를 a*x라고 하기는 했지만, 엄연히 말하면 def에서 b까지 제시하였고, curve_fit 자체가 기본적으로 b값까지 제시하지 않으면 오류가 나기 때문에, bounds에서도 b값을 넣어 주어야 합니다. 

 

첫 번째 [ ]는 a, b의 최소값이며

두 번째 [ ]는 a, b의 최대값을 의미합니다. 

 

b값에 -0.0000001, 0.0000001 같이 집어넣은건, y절편 값을 최대한 0에 맞추기 위해 구간을 저리 지정하였습니다.

 

 

4) 추세선(fitting line) 그리기

 

plt.plot(x,objective(x,*popt),color='red',linestyle='-', linewidth=0.5)

 

이제 curve_fit을 통해, 저 데이터에 적합한 a, b 값이 계산되었고, 이에 따른 일종의 오차(정확히는 공분산값)까지 계산되었습니다. 이 결과를 바탕으로 그래프를 그리면 됩니다.

 

curve_fit에서 계산되는 데이터는 크게 a,b에 해당되는 값 하나와 오차 값 하나입니다. 그래서 위에서 popt, pcov = 이라고 2개의 변수를 잡아준 건데, 첫 번째 변수인 popt에는 a, b에 해당하는 fittig 결과가, 두 번째 변수인 pcov에는 오차(공분산)값이 지정되게 됩니다.

 

이제 이 값을 plot 하면 됩니다. 일반적인 plot과 동일한데, 아주 약간 다릅니다.

 

일반적으로 plot 다음에는 x값과 y값이 나옵니다. 예를들어 plt.plot(x,y) 이런식으로요. 그런데 이번에는 x값에 fitting한 새로운 y값이 필요한 것입니다. 그 값은 popt에 저장되어 있습니다. 값이 2개로. 그런데 이건 경우에 따라 2개 이상이 될 수 도 있습니다. 예를들어 2차함수라면 최대 3개까지도 가능합니다. 그래서 코드가 저리 짜였는데, 코드를 잘 보면

 

plt. plot 다음으로, 첫 번째는 x값이, 두 번째는 objective(x, *popt)가 들어가 있습니다. 여기서 x값을 a*x라는 objective 함수에 집어넣었을 때 새롭게 계산된 y값이 자리하게 됩니다. 그 뒤로는 단순히 그래프 색이나 라인 스타일, 라인 폭 등을 지정해 주는 값들입니다.

 

 

5) 오차를 최소화 하고 확인하려면

 

실제 fitting을 수행할 때에는 실측 결과와 모델간의 오차가 가장 적을 때 까지 모델을 수정해가며 반복작업을 해야합니다. 여기서야 어차피 맛보기로 하는 것이지만, 좀 더 과학적인 작업을 원하신다면, pcov에 저장되어있는 오차값을 꼼꼼히 확인하면서 bounds의 최소 최대값을 지정해 줘야합니다. (만약 bounds에서 최소, 최대값을 지정해주지 않으면 파이썬이 알아서 계산하긴 합니다만, y절편까지 계산하는 등 불필요한 작업을 하여 마음에 안들 수 있습니다.)

 

pcov에서는 이 경우에 총 4개의 공분산 값을 계산해 주는데,  우리에게 필요한건 x에 대한 y의 변화값에 대한 공분산 값입니다. pcov에 저장되는 공분산값에 대한 글은 이 글의 범위를 너무 벗어나기 때문에, 여기서 자세히 논하지는 않겠습니다. 어쨌든 우리에게 필요한 값의 표준편차를 계산하는 코드는 아래 코드입니다. 

 

perr = np.sqrt(np.diag(pcov))

 

요 코드를 짠 뒤, perr값을 확인하며 표준편차가 가장 작아질 때 까지 bounds의 최소 최대값을 수정해가면서 피팅을 반복하면 될 것 입니다.

 

6) 긴 글을 마치며

  예제 중심으로 파이썬과 지구과학을 설명하다 보니, 코드에 대한 자세한 설명을 하기 어렵습니다. 지구과학과 파이썬을 연계하는것이 목적이니, 코드에 대한 자세한 설명은 각 라이브러리 홈페이지나 다른 인터넷 글을 보시는 수 밖에 없습니다..ㅠㅠ 마지막으로 본 포스팅을 쓰기 위해, 적절한 은하 데이터를 찾는 과정에서 탈락해버린 그래프 하나를 마저 소개합니다. 이 녀석은 버리긴 했지만, 그리긴 했으니 여기다 함께 실어 두겠습니다.