이번 포스팅에서는 Python matplotlib 모듈을 사용해서 그래프를 그렸을 때 범례를 추가(adding a legend on the Axes)하는 4가지 방법을 소개하겠습니다. 그리고 범례의 위치, 범례 글자 크기, 범례 배경색, 범례 테두리 색, 범례 그림자 등을 설정하는 방법을 소개하겠습니다. 

 

(1) 범례 추가하기 (adding a legend on the Axes)

    : legend(), legend(handles, labels), legend(handles=handles), legend(labels)

(2) 범례 스타일 설정하기
    : location, fontsize, facecolor, edgecolor, shadow

 

 

(1) 범례 추가하기 (adding a legend on the Axes)

 

matplotlib 으로 그래프를 그린 후 범례를 추가하는 방법에는 (1-1) legend(), (1-2) line.set_label(), (1-3) legend(handles, labels), (1-4) legend(handles=handles) 의 4가지 방법이 있습니다. 이중에서 그래프를 그릴 때 label 을 입력해주고, 이어서 ax.legend() 나 plt.legend() 를 실행해서 자동으로 범례를 만드는 방법이 가장 쉽습니다. 차례대로 4가지 방법을 소개하겠습니다.  결과는 모두 동일하므로 범례가 추가된 그래프는 (1-1)번에만 소개합니다. 

 

먼저, 아래의 2개 방법은 자동으로 범례를 탐지해서 범례를 생성해주는 예제입니다. 

 

(1-1) ax.plot(label='text...') --> ax.legend() 또는 plt.legend()

 

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1)
x = np.arange(100)
y1 = np.random.normal(0, 1, 100)
y2 = np.random.normal(0, 1, 100).cumsum()

## 1. Automatic detection of elements to be shown in the legend
## 1-1. label
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

## labels: A list of labels to show next to the artists. 
ax.plot(x, y1, 'b-', label='Company 1')
ax.plot(x, y2, 'r--', label='Company 2')
ax.legend()
plt.show()

matplotlib legend

 

 

(1-2) line.set_label('text...') --> ax.legend()

 

## 1-2. line.set_label()

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

line, = ax.plot(x, y1, 'b-')
line.set_label('Company 1')

line, = ax.plot(x, y2, 'r--')
line.set_label('Company 2')

ax.legend()
plt.show()

 

 

위의 두 방법은 범례를 자동 탐지해서 생성해줬다면, 아래의 방법부터는 명시적으로 범례의 handles (선 모양과 색깔 등) 와 labels (범례 이름)를 지정해줍니다.

 

(1-3) ax.legend(handles, labels)

 

## 2. Explicitly listing the artists and labels in the legend
## to pass an iterable of legend artists followed by an iterable of legend labels respectively
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

line1, = ax.plot(x, y1, 'b-')
line2, = ax.plot(x, y2, 'r--')

ax.legend([line1, line2], ['Company 1', 'Company 2'])
plt.show()

 

 

 

(1-4) ax.legend(handles=handles)

 

## 3. Explicitly listing the artists in the legend
## the labels are taken from the artists' label properties.
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

line1, = ax.plot(x, y1, 'b-', label='Company 1')
line2, = ax.plot(x, y2, 'r--', label='Company 2')

## handles: A list of Artists (lines, patches) to be added to the legend. 
ax.legend(handles=[line1, line2])
plt.show()

 

 

 

(2) 범례 스타일 설정하기
    : location, fontsize, facecolor, edgecolor, shadow

 

범례의 위치(location) 을 지정하려면 loc 매개변수에 location strings 인 'best', 'upper right', 'upper left', 'lower left', 'lower right', 'right', 'center left', 'center right', 'lower center', 'upper center', 'center' 중에서 하나를 선택해서 입력해주면 됩니다. 기본설정값은 loc='best' 로서 컴퓨터가 알아서 가장 적절한 여백을 가지는 위치에 범례를 생성해주는데요, 나름 똑똑하게 잘 위치를 선택해서 범례를 추가해줍니다. 

 

그밖에 범례의 폰트 크기(fontsize), 범례의 배경 색깔(facecolor, 기본설정값은 'white'), 범례의 테두리선 색깔(edgecolor), 범례 상자의 그림자 여부(shadow, 기본설정값은 'False') 등도 사용자가 직접 지정해줄 수 있습니다. (굳이 범례의 스타일을 꾸미는데 시간과 노력을 많이 쏟을 이유는 없을 것 같기는 합니다만....^^;)

 

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

ax.plot(x, y1, 'b-', label='Company 1')
ax.plot(x, y2, 'r--', label='Company 2')
ax.legend(
    ## location strings
    ## : 'best', 'upper right', 'upper left', 'lower left', 'lower right'
    ##   'right', 'center left', 'center right', 'lower center', 'upper center', 'center'
    loc='center left', # default: loc='best'
    ## fontsize
    ## # 'xx-small','x-small','small','medium','x-large','xx-large'
    fontsize='large', 
    facecolor='yellow', # background color
    edgecolor='black',
    shadow=True) # shadow behind the legend if True

ax.set_title('Adding a Legend with color and shadow', fontsize=18)
plt.show()

matplotlib adding a legend with color and shadow at center left location

 

 

[ Reference ]

* matplotlib.pyplot.legend() : https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.legend.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요. :-)

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 Python matplotlib 모듈로 그래프를 그릴 때,  

 

(1) 눈금 설정하기 : set_xticks(), set_yticks()

(2) 눈금 이름 설정하기 : set_xticklabels(), set_yticklabels()

(3) 축 이름 설정하기 : set_xlabel(), set_ylabel()

(4) 제목 설정하기 : set_title()

 

하는 방법을 소개하겠습니다. 

 

 

먼저, numpy 모듈을 사용해서 표준정규분포 X~N(0, 1) 로부터 난수 100개를 생성해서 matplotlib의 기본 설정으로 선 그래프를 그려보겠습니다. matplotlib의 Figure 객체를 만들고, fig.add_subplot(1, 1, 1) 로 하위 플롯을 하나 만든 다음에 선 그래프를 그렸습니다. 

 

import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

np.random.seed(1)
x=np.arange(100)
y=np.random.normal(0, 1, 100)

## plot with default setting
ax.plot(x, y)
plt.show()

 

 

(1) 눈금 설정하기 : set_xticks(), set_yticks()

 

X축의 눈금을 [0, 50, 100] 으로 설정하고, Y축의 눈금은 [-2, 0, 2] 로 설정해보겠습니다. 

 

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

ax.plot(x, y)

## setting xticks, yticks
ax.set_xticks([0, 50, 100])
ax.set_yticks([-2, 0, 2])

plt.show()

matplotlib.pyplot : ax.set_xticks(), ax.set_yticks()

 

 

 

(2) 눈금 이름 설정하기 : set_xticklabels(), set_yticklabels()

 

X축과 Y축에 눈금이름(xticklabel, yticklabel)을 설정해줄 수도 있습니다. 이때 set_xticks(), set_yticks() 메소드로 눈금의 위치를 먼저 설정해주고, 이어서 같은 개수의 눈금 이름을 가지고 set_xticklabels(), set_yticklabels() 메소드로 눈금 이름을 설정해주면 됩니다. 

 

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

ax.plot(x, y)

## setting xticks, yticks with labels
ax.set_xticks([0, 50, 100])
ax.set_yticks([-2, 0, 2])
ax.set_xticklabels(['start', 'middel', 'end'], fontsize=12)
ax.set_yticklabels(['low', 'zero', 'high'], fontsize=12)

plt.show()

 

 

 

(3) 축 이름 설정하기 : set_xlabel(), set_ylabel()

 

set_xlabel(), set_ylabel() 메소드로 축 이름(xlabel, ylabel)을 설정해 줄 때 fontsize 매개변수를 사용해서 글자 크기를 조정할 수 있습니다. 

 

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

ax.plot(x, y)

## setting xticks, yticks with labels
ax.set_xticks([0, 50, 100])
ax.set_yticks([-2, 0, 2])
ax.set_xticklabels(['start', 'middel', 'end'], fontsize=12)
ax.set_yticklabels(['low', 'zero', 'high'], fontsize=12)

## setting xlabel, ylabel
ax.set_xlabel('Steps', fontsize=16)
ax.set_ylabel('Value', fontsize=16)

plt.show()

 

 

 

(4) 제목 설정하기 : set_title()

 

제목의 텍스트, 폰트 크기(fontsize), 위치(loc={'center', 'left', 'right'}) 를 설정해 줄 수 있습니다. 

 

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1)

ax.plot(x, y)

## setting xticks, yticks with labels
ax.set_xticks([0, 50, 100])
ax.set_yticks([-2, 0, 2])
ax.set_xticklabels(['start', 'middel', 'end'], fontsize=12)
ax.set_yticklabels(['low', 'zero', 'high'], fontsize=12)

## xlabel, ylabel
ax.set_xlabel('Steps', fontsize=16)
ax.set_ylabel('Value', fontsize=16)

# title
ax.set_title('Plo with ticks and labels', 
             fontsize=20, 
             loc='left') # 'center', 'right'

plt.show()

matplotlib: setting title, xlabel, xticks, xticklabels

 

 

 

[ Reference ]

* ax.set_xticks(): https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xticks.html

* ax.set_xticklabels(): https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xticklabels.html

* ax.set_xlabel(): https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xlabel.html

* ax.set_title(): https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.axes.Axes.set_title.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요! :-)

 

728x90
반응형
Posted by Rfriend
,

여러개의 그래프를 위/아래 또는 왼쪽/오른쪽으로 붙여서 시각화를 한 후에 비교를 해보고 싶을 때가 있습니다. 그리고 측정 단위가 서로 같을 경우에는 x축, y축을 여러개의 그래프와 공유하면서 동일한 scale의 눈금으로 시각화 함으로써 그래프 간 비교를 더 쉽고 정확하게 할 수 있습니다. 

 

이번 포스팅에서는

  (1) 여러개의 하위 플롯 그리기: plt.subplots(nrows=2, ncols=2)

  (2) X축과 Y축을 공유하기
        : plt.subplots(sharex=True. sharey=True)

  (3) 여러개의 하위 플롯들 간의 간격을 조절해서 붙이기

        : (3-1) plt.subplots_adjust(wspace=0, hspace=0)

        : (3-2) plt.tight_layout(h_pad=-1, w_pad=-1)

하는 방법을 소개하겠습니다. 

 

 

 

(1) 여러개의 하위 플롯 그릭: plt.subplots(nrows=2, ncols=2)

 

먼저, plt.subplots(nrows=2, ncols=2) 로 2개 행, 2개 열의 layout 을 가지는 하위 플롯을 그려보겠습니다. 평균=0, 표준편차=1을 가지는 표준정규분포로 부터 100개의 난수를 생성해서 파란색(color='b')의 약간 투명한(alpha=0.5) 히스토그램을 그려보겠습니다.

 

for loop 을 돌면서 2 x 2 의 레이아웃에 하위 플롯이 그려질때는 왼쪽 상단이 1번, 오른쪽 상단이 2번, 왼쪽 하단이 3번, 오른쪽 하단이 4번의 순서로 하위 그래프가 그려집니다.  

 

이렇게 하위 플롯이 그려질 때 자동으로 하위 플롯 간 여백(padding)을 두고 그래프가 그려지며, X축과 Y축은 공유되지 않고 각각 축을 가지고 그려집니다.  

 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

## subplots with 2 by 2
## (1) there are paddings between subplots
## (2) x and y axes are not shared
fig, axes = plt.subplots(
    nrows=2, 
    ncols=2, 
    figsize=(10, 10))

for i in range(2):
    for j in range(2):
        axes[i, j].hist(
            np.random.normal(loc=0, scale=1, size=100), 
            color='b', 
            alpha=0.5)

subplots with padding and x, y axes

 

 

(2) X축과 Y축을 공유하기
      : plt.subplots(sharex=True. sharey=True)

(3) 여러개의 하위 플롯들 간의 간격을 조절해서 붙이기

      : (3-1) plt.subplots_adjust(wspace=0, hspace=0)

 

여러개의 하위 플롯 간에 X축과 Y축을 공유하려면 plt.subplots()의 옵션 중에서 sharex=True, sharey=True를 설정해주면 됩니다. 

 

그리고 하위 플롯들 간의 간격을 없애서 그래프를 서로 붙여서 그리고 싶다면 두가지 방법이 있습니다. 먼저, plt.subplots_adjust(wspace=0, hspace=0) 에서 wspace 는 폭의 간격(width space), hspace 는 높이의 간격(height space) 을 설정하는데 사용합니다. 

 

fig, axes = plt.subplots(
    nrows=2, ncols=2, 
    sharex=True, # sharing properties among x axes
    sharey=True, # sharing properties among y axes 
    figsize=(10, 10))

for i in range(2):
    for j in range(2):
        axes[i, j].hist(
            np.random.normal(loc=0, scale=1, size=100), 
            color='b', 
            alpha=0.5)
        
## adjust the subplot layout 
plt.subplots_adjust(
    wspace=0, # the width of the padding between subplots 
    hspace=0) # the height of the padding between subplots

matplotlib subplots adjust

 

 

아래 예제에서는 plt.subplots(sharex=False, sharey=True) 로 해서 X축은 공유하지 않고 Y축만 공유하도록 했습니다. 

그리고 plt.subplots_adjust(wspace=0, hspace=0.2) 로 해서 높이의 간격(height space)에만 0.2 만큼의 간격을 부여해주었습니다. 

 

fig, axes = plt.subplots(
    nrows=2, ncols=2, 
    sharex=False, # sharing properties among x axes
    sharey=True, # sharing properties among y axes 
    figsize=(10, 10))

for i in range(2):
    for j in range(2):
        axes[i, j].hist(
            np.random.normal(loc=0, scale=1, size=100), 
            color='b', 
            alpha=0.5)
        
## adjust the subplot layout 
plt.subplots_adjust(
    wspace=0, # the width of the padding between subplots 
    hspace=0.2) # the height of the padding between subplots

 

 

 

(3) 여러개의 하위 플롯들 간의 간격을 조절해서 붙이기

      : (3-2) plt.tight_layout(h_pad=-1, w_pad=-1)

 

하위 플롯 간 간격을 조절하는 두번째 방법으로는 plt.tight_layout(h_pad, w_pad) 을 사용하는 방법입니다. plt.tight_layout(h_pad=-1, w_pad=-1)로 설정해서 위의 2번에서 했던 것처럼 4개의 하위 플롯 간에 간격이 없이 모두 붙여서 그려보겠습니다. (참고: h_pad=0, w_pad=0 으로 설정하면 하위 플롯간에 약간의 간격이 있습니다.)

 

fig, axes = plt.subplots(
    nrows=2, ncols=2, 
    sharex=True, # sharing properties among x axes
    sharey=True, # sharing properties among y axes 
    figsize=(10, 10))

for i in range(2):
    for j in range(2):
        axes[i, j].hist(
            np.random.normal(loc=0, scale=1, size=100), 
            color='b', 
            alpha=0.5)

## adjusting the padding between and around subplots
plt.tight_layout(
    h_pad=-1, # padding height between edges of adjacent subplots
    w_pad=-1) # padding width between edges of adjacent subplots

matplotlib subplots tight_layout

 

 

X축은 공유하지 않고 Y축만 공유하며, plt.tight_layout(h_pad=3, w_pad=0) 으로 설정해서 높이 간격을 벌려보겠습니다. 그리고 하위 플롯이 그려지는 순서대로 'blue', 'red', 'yellow', 'black' 의 색깔을 입혀보겠습니다. 

 

fig, axes = plt.subplots(
    nrows=2, ncols=2, 
    sharex=False, # sharing properties among x axes
    sharey=True, # sharing properties among y axes 
    figsize=(10, 10))

color = ['blue', 'red', 'yellow', 'black']
k=0
for i in range(2):
    for j in range(2):
        axes[i, j].hist(
            np.random.normal(loc=0, scale=1, size=100), 
            color=color[k], 
            alpha=0.5)
        k +=1

## adjusting the padding between and around subplots
plt.tight_layout(
    h_pad=3, # padding height between edges of adjacent subplots
    w_pad=0) # padding width between edges of adjacent subplots

 

[ Reference ]

- plt.subplots(): https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html
- plt.subplots_adjust(): https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots_adjust.html
- plt.tight_layout(): https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.tight_layout.html

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 Python의 pandas DataFrame 에서 여러개 칼럼의 조건을 일부(any)만 또는 전부(all) 만족하는 행 가져오기하는 방법을 소개하겠습니다. 

 

pandas DataFrame의 any(), all() 메소드를 이용하면 매우 간단하게 조건을 만족하는 행을 색인해서 가져올 수 있습니다. 

 

(1) pandas DataFrame 의 칼럼 중에서 1개 이상의 칼럼이 True 인 행 가져오기: pd.DataFrame.any(axis=1)

(2) pandas DataFrame 의 칼럼 중에서 모든 칼럼이 True 인 행 가져오기: pd.DataFrame.all(axis=1)

 

 

pandas DataFrae any(axis=1), all(axis=1)

 

 

 

먼저, 예제로 사용할 간단한 pandas DataFrame 을 만들어 보겠습니다. 4개의 칼럼을 가지고 있고, 결측값도 하나 넣어습니다. 

 

import pandas as pd
import numpy as np

## creating a sample DataFrame with 4 columns
df = pd.DataFrame({
    'x1': [0, 1, 2, 3, 4], 
    'x2': [-2, -1, 0, 1, 3], 
    'x3': [-4, -3, -2, -1, -4], 
    'x4': [np.nan, 0, 2, 3, -10]
})

print(df)
#    x1  x2  x3    x4
# 0   0  -2  -4   NaN
# 1   1  -1  -3   0.0
# 2   2   0  -2   2.0
# 3   3   1  -1   3.0
# 4   4   3  -4 -10.0

 

 

 

(1) pandas DataFrame 의 칼럼 중에서 1개 이상의 칼럼이 True 인 행 가져오기: pd.DataFrame.any(axis=1)

 

아래처럼 np.abs(df) > 2 를 하면 모든 칼럼의 모든 행에 대해서 절대값(absolute value)이 2보다 큰지 아닌지 여부에 대해 True/False 블리언 값을 반환합니다. 

 

## returns boolean for all columns and rows
np.abs(df) > 2

# 	x1	   x2	   x3	   x4
# 0	False	False	True	False
# 1	False	False	True	False
# 2	False	False	False	False
# 3	True	False	False	True
# 4	True	True	True	True

 

 

이제 칼럼 4개 중에서 절대값(absolute value)이 2보다 큰 값이 단 하나라도 존재하는 행을 가져와 보겠습니다. 이때 '칼럼들에 대해 단 하나라도 존재하면'이라는 조건 판단은 pandas.DataFrame.any(axis=1) 메소드를 사용하면 편리합니다. 

 

any(axis =1) 에서 axis=1 을 설정해주면 칼럼 축을 기준으로 조건 만족여부를 평가합니다. 기본설정값이 axis=0 이므로 반드시 명시적으로 any(axis=1) 처럼 축을 지정해주어야 합니다. 

 

결측값이 있어도 다른 칼럼에서 조건을 만족하면 해당 행을 가져옵니다. 

 

## (1) DataFrame.any(axis=0, bool_only=None, skipna=True, level=None, **kwargs)
## pd.DataFrame.any(): Return whether any element is True, potentially over an axis.

df[(np.abs(df) > 2).any(axis=1)]
#    x1	x2	x3	x4
# 0	 0	-2	-4	NaN
# 1	 1	-1	-3	0.0
# 3	 3	1	-1	3.0
# 4	 4	3	-4	-10.0

 

 

 

pandas.DataFrame.any(axis=1) 메소드를 사용하지 않고, 아래처럼 블리언 값(True=1, False=0)을 칼럼 축으로 더해서(sum(axis=1)), 그 더한 값이 0보다 큰 행을 인덱싱해서 가져오는 방법을 써도 되긴 합니다. 

 

## or equivalantly
df[(np.abs(df) > 2).sum(axis=1) > 0]

#    x1	x2	x3	x4
# 0	 0	-2	-4	NaN
# 1	 1	-1	-3	0.0
# 3	 3	1	-1	3.0
# 4	 4	3	-4	-10.0

 

 

 

 

(2) pandas DataFrame 의 칼럼 중에서 모든 칼럼이 True 인 행 가져오기: pd.DataFrame.all(axis=1)

 

이번에는 pandas.DataFrame.all(axis=1)을 이용해서 DataFrame에 있는 4개의 모든 칼럼이 조건을 만족하는 행만 가져오기를 해보겠습니다. 

 

## DataFrame.all(axis=0, bool_only=None, skipna=True, level=None, **kwargs)
## Return whether all elements are True, potentially over an axis.

df[(np.abs(df) > 2).all(axis=1)]
#    x1	  x2	 x3	  x4
# 4	 4	3	-4	-10.0

 

 

 

아래는 pandas.DataFrame.all() 메소드를 사용하지 않고, 대신에 조건 만족여부에 대한 블리언 값을 칼럼 축으로 전부 더한 후, 이 블리언 합이 칼럼 개수와 동일한 행을 가져와본 것입니다. 위의 pandas.DataFrame.all(axis=1) 보다 코드가 좀더 길고 복잡합니다. 

 

## or, equivalently
df[(np.abs(df) > 2).sum(axis=1) == df.shape[1]]

#     x1	 x2	 x3	 x4
# 4	 4	3	-4	-10.0

 

[ Reference ]

* pandas.DataFrame.any() : https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.any.html  
* pandas.DataFrame.all(): https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.all.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요! :-)

 

728x90
반응형
Posted by Rfriend
,

프로그래밍 코드를 짜다보면 수행 절차나 방법, 사용하는 메소드에 따라서 수행 시간이 차이가 나는 경우가 종종 있습니다. 그리고 성능이 중요해서 여러가지 방법을 테스트해보면서 가장 실행시간이 짧도록 튜닝하면서 최적화하기도 합니다. 

 

이번 포스팅에서는 Python에서 코드를 실행시켰을 때 소요된 시간을 측정하는 2가지 방법을 소개하겠습니다. 

 

(1) datetime.now() 메소드 이용해서 실행 시간 측정하기

(2) %timeit 로 실행 시간 측정하기

 

 

python : measuring the execution time of code snippets

 

먼저, 예제로 사용할 샘플 데이터셋으로서, 1억개의 값을 가지는 xarr, yarr 의 두개의 배열(array)를 만들어 보겠습니다. 그리고 배열 내 각 1억개의 값 별로 True/False 의 조건값을 가지는 cond 라는 배열도 난수를 생성시켜서 만들어보겠습니다. 

 

import numpy as np

## generating sample array with 100 million values
xarr = np.arange(100000000)
yarr = np.zeros(100000000)
cond = np.where(np.random.randn(100000000)>0, True, False)


cond[:10]
# [Out] array([False,  True,  True, False, False,  True,  True,  True,  
#              True, True])

 

 

위에서 만든 1억개의 원소를 가지는 배열을 가지고 조건값으로 True/False 블리언 값 여부에 따라서 True 조건값 이면 xarr 배열 내 값을 가지고, False 조건값이면 yarr 배열 내 값을 가지는 새로운 배열을 만들어보겠습니다. 이때 (1) List Comprehension 방법과, (2) NumPy의 Vectorized Operations 방법 간 수행 시간을 측정해서 어떤 방법이 더 빠른지 성능을 비교해보겠습니다.

(물론, Vectorized Operations이 for loop 순환문을 사용하는 List Comprehension보다 훨~씬 빠릅니다! 눈으로 직접 확인해 보겠습니다. )

 

## Let's compare the elapsed time between 2 methods 
## (list comprehension vs. vectorized operations)

## (1) List Comprehension
new_arr = [(x if c else y) for (x, y, c) in zip(xarr, yarr, cond)]

## (2) Vectorized Operations in NumPy 
new_arr = np.where(cond, xarr, yarr)

 

 

 

(1) datetime.now() 메소드 이용해서 실행 시간 측정하기

 

datetime 모듈은 날짜, 시간, 시간대(time zone) 등을 다루는데 사용하는 모듈입니다 datetime.now() 메소드는 현재의 로컬 날짜와 시간을 반환합니다. 실행 시간을 측정할 코드 바로 앞에 start_time = datetime.now() 로 시작 날짜/시간을 측정해놓고, 실행할 코드가 끝난 다음 줄에 time_elapsed = datetime.now() - start_time 으로 '끝난 날짜/시간'에서 '시작 날짜/시간'을 빼주면 '코드 실행 시간'을 계산할 수 있습니다. 

 

아래 결과를 비교해보면 알 수 있는 것처럼, for loop 순환문을 사용하는 List Comprehension 방법보다 NumPy의 Vectorized Operation이 약 38배 빠른 것으로 나오네요. 

 

## (1) -- measuring the elapsed time using datetime

## (a) List Comprehension
from datetime import datetime
start_time = datetime.now() 
list_comp_for_loop = [(x if c else y) for (x, y, c) in zip(xarr, yarr, cond)]
time_elapsed = datetime.now() - start_time 

print('Time elapsed (hh:mm:ss.ms) {}'.format(time_elapsed))
# Time elapsed (hh:mm:ss.ms) 0:00:17.753036

np.array(list_comp_for_loop)[:10]
# array([0., 1., 2., 0., 0., 5., 6., 7., 8., 9.])



## (b) Vectorized Operations in NumPy 
start_time = datetime.now() 
np_where_vectorization = np.where(cond, xarr, yarr)
time_elapsed = datetime.now() - start_time 

print('Time elapsed (hh:mm:ss.ms) {}'.format(time_elapsed))
# Time elapsed (hh:mm:ss.ms) 0:00:00.462215

np_where_vectorization[:10]
# array([0., 1., 2., 0., 0., 5., 6., 7., 8., 9.])

 

 

(2) %timeit 로 실행 시간 측정하기

 

다음으로 Python timeit 모듈을 사용해서 짧은 코드의 실행 시간을 측정해보겠습니다. timeit 모듈은 터미널의 command line 과 Python IDE 에서 호출 가능한 형태의 코드 둘 다 사용이 가능합니다. 

 

아래에는 Jupyter Notebook에서 %timeit [small code snippets] 로 코드 수행 시간을 측정해본 예인데요, 여러번 수행을 해서 평균 수행 시간과 표준편차를 보여주는 특징이 있습니다. 

 

## (2) measuring the elapsed time using timeit

## (a) List Comprehension
import timeit

%timeit list_comp_for_loop = [(x if c else y) for (x, y, c) in zip(xarr, yarr, cond)]
# 17.1 s ± 238 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


## (b) Vectorized Operations in NumPy 
%timeit np_where_vectorization = np.where(cond, xarr, yarr)
# 468 ms ± 8.75 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

 

 

[Reference]

* Python datetime: https://docs.python.org/3/library/datetime.html

* Python timeit: "measuring the execution time of small code snippets"
   : https://docs.python.org/3/library/timeit.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요! 

 

728x90
반응형
Posted by Rfriend
,

Word counts 할 때 많이 사용하는 코드인데요, 이번 포스팅에서는 

  (1) 리스트에서 원소별 개수를 세서 {Key:Value} 쌍의 Dictionary를 만들고 

  (2) 원소별 개수를 세어놓은 Dictionary에서 개수 상위 n 개의 {Key:Value} 쌍을 가져오기

하는 방법을 소개하겠습니다. 

 

 

counts dictionary, getting top n

 

(1) 리스트에서 원소별 개수를 세서 {Key:Value} 쌍의 Dictionary를 만들기

 

먼저, 예제로 사용할 간단한 리스트를 만들어보겠습니다. 

 

## creating sample lists
my_list = ['a', 'f', 'a', 'b', 'a', 'a', 'c', 'b', 
           'c', 'e', 'a', 'c', 'b', 'f', 'c']
           
print(my_list)
# ['a', 'f', 'a', 'b', 'a', 'a', 'c', 'b', 'c', 'e', 'a', 'c', 'b', 'f', 'c']

 

 

다음으로, 원소별 개수를 세서 저장할 비어있는 Dictionary 인 counts={} 를 만들어놓고, for loop 순환문으로 리스트의 원소를 하나씩 순서대로 가져다가 Dictionary counts 의 Key 값에 해당 원소가 들어있으면 +1을 하고, Key 값에 해당 원소가 안들어 있으면 해당 원소를 Key 값으로 등록하고 1 을 값으로 입력해 줍니다.  

 

def get_counts(seq): 
    counts = {}
    for x in seq:
        if x in counts:
            counts[x] += 1
        else:
            counts[x] = 1
    return counts
    
 
counts = get_counts(my_list)


print(counts)
# {'a': 5, 'f': 2, 'b': 3, 'c': 4, 'e': 1}


## access value by key
counts['a']
# 5

 

 

 

(2) 원소별 개수를 세어놓은 Dictionary에서 개수 상위 n 개의 {Key:Value} 쌍을 가져오기

 

Dictionary를 정렬하는 방법에 따라서 두 가지 방법이 있습니다.

 

(a) sorted() 메소드를 이용해서 key=lambda x: x[1] 로 해서 정렬 기준을 Dictionary의 Value 로 하여 내림차순으로 정렬(reverse=True) 하고, 상위 n 개까지만 슬라이싱해서 가져오는 방법입니다. 

 

## way 1
## reference: https://rfriend.tistory.com/473
def top_n(count_dict, n=3):
    return sorted(count_dict.items(), reverse=True, key=lambda x: x[1])[:n]
    

## getting top 2
top_n(counts, n=2)
# [('a', 5), ('c', 4)]

 

 

 

(b) 아래는 dict.items() 로 (Key, Value) 쌍을 for loop 문을 돌리면서 (Value, Key) 로 순서를 바꿔서 리스트 [] 로 만들고 (list comprehension), 이 리스트에 대해서 sort(reverse=True) 로 Value 를 기준으로 내림차순 정렬한 후에, 상위 n 개까지만 슬라이싱해서 가져오는 방법입니다. 

 

## way2
## reference: https://rfriend.tistory.com/281
def top_n2(count_dict, n=3):
    val_key = [(v, k) for k, v in count_dict.items()]
    val_key.sort(reverse=True)
    return val_key[:n]
    
## getting top 2
top_n2(counts, n=2)
# [(5, 'a'), (4, 'c')]

 

 

[Reference]

* Dictionary 정렬: https://rfriend.tistory.com/473

* List 정렬: https://rfriend.tistory.com/281

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요. 

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 시계열 데이터에서

  (1) TimeStamp 행별로 칼럼별 비율을 구하고 

  (2) 시도표 (time series plot) 를 그리기

하는 방법을 소개하겠습니다. 

 

먼저, 예제로 사용할 간단한 pandas DataFrame을 만들어보겠습니다. index 로 2000년 ~ 2021년까지의 년도를 사용하고, 성별로 'M', 'F'의 두 개의 칼럼에 포아송분포로 부터 난수를 발생시켜서 만든 도수(frequency)를 가지는 DataFrame 입니다. 

 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


## creating a sample DataFrame
ts = np.arange(2000, 2022) # from year 2000 to 2021

print(ts)
# [2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
#  2014 2015 2016 2017 2018 2019 2020 2021]

np.random.seed(1) # for reproducibility
M = np.arange(len(ts)) + np.random.poisson(lam=10, size=len(ts))
F = np.arange(len(ts))[::-1] + np.random.poisson(lam=2, size=len(ts))

df = pd.DataFrame({'M': M, 'F': F}, index=ts)

df.head()
#       M	F
# 2000	9	21
# 2001	7	24
# 2002	9	20
# 2003	12	19
# 2004	13	18

 

 

(1) TimeStamp 행별로 칼럼별 비율을 구하기

 

먼저, pandas DataFrame 에서 합(sum)을 구할 때 각 TimeStamp 별로 칼럼 축(axis = 1) 으로 합을 구해보겠습니다. 

 

## calculating summation by rows
df.sum(axis=1).head()

# 2000    30
# 2001    31
# 2002    29
# 2003    31
# 2004    31
# dtype: int64

 

 

참고로, index 축으로 칼럼별 합을 구할 때는 df.sum(axis=0) 을 해주면 됩니다. sum(axis=0) 이 기본설정값이므로 df.sum() 하면 동일한 결과가 나옵니다. 

 

## summation by index axis
df.sum(axis=0) # default setting

# M    426
# F    274
# dtype: int64

 

 

pandas DataFrame에서 div() 메소드를 사용하면 각 원소를 특정 값으로 나눌 수 있습니다. 가령, 위의 예제 df DataFrame의 각 원소를 10으로 나눈다고 하면 아래처럼 df.div(10) 이라고 해주면 됩니다. (나누어주는 값 '10' 이 broadcasting 되어서 각 원소를 나누어주었음.)

 

## pd.DataFrame.div()
## : Get Floating division of dataframe and other, 
##   element-wise (binary operator truediv).
df.div(10).head()

#         M	F
# 2000	0.9	2.1
# 2001	0.7	2.4
# 2002	0.9	2.0
# 2003	1.2	1.9
# 2004	1.3	1.8

 

 

이제 df DataFrame의 각 원소를 각 원소가 속한 TimeStamp별로 칼럼 축(axis=1)으로 합한 값(df.sum(axis=1))으로 나누어주면 우리가 구하고자 하는 각 TimeStamp별 칼럼별 비율을 구할 수 있습니다. 

 

df.div(df.sum(axis=1), axis=0).head()

#         M	        F
# 2000	0.300000	0.700000
# 2001	0.225806	0.774194
# 2002	0.310345	0.689655
# 2003	0.387097	0.612903
# 2004	0.419355	0.580645

 

 

 

(2) 시도표 (time series plot) 를 그리기

 

pandas DataFrame 의 plot() 메소드를 사용하면 편리하게 시계열 도표를 그릴 수 있습니다. 이때 성별을 나타내는 칼럼 'M', 'F' 별로 선의 모양(line type)과 색깔(color) 을 style={'M': 'b--', 'F': 'r-'} 매개변수를 사용해서 다르게 해서 그려보겠습니다. ('M' 은 파란색 점선, 'F' 는 빨간색 실선)

 

df.div(df.sum(1), axis=0).plot(
    style={'M': 'b--', 'F': 'r-'}, 
    figsize=(12, 8), 
    title='Proportion Trend by Gender')

plt.show()

proportion trend plot by gender

 

[ Reference ]

* pandas.DataFrame.div(): https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.div.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

이전 포스팅에서는 

 (i) 정상확률과정(stationary process)의 정의 (https://rfriend.tistory.com/691)

 (ii) 통계적 가설 검증을 통한 시계열 정상성(stationarity test) 여부 확인 (https://rfriend.tistory.com/694)

하는 방법을 소개하였습니다. 

 

ARIMA 모형과 같은 통계적 시계열 예측 모델의 경우 시계열데이터의 정상성 가정을 충족시켜야 합니다. 따라서 만약 시계열 데이터가 비정상 확률 과정 (non-stationary process) 이라면, 먼저 시계열 데이터 변환을 통해서 정상성(stationarity)을 충족시켜주어야 ARIMA 모형을 적합할 수 있습니다

 

이번 포스팅에서는 Python을 사용하여 

(1) 분산이 고정적이지 않은 경우 분산 안정화 변환 (variance stabilizing transformation, VST)

(2) 추세가 있는 경우 차분을 통한 추세 제거 (de-trend by differencing)

(3) 계절성이 있는 경우 계절 차분을 통한 계절성 제거 (de-seasonality by seaanl differencing)

하는 방법을 소개하겠습니다. 

 

 

[ 비정상확률과정을 정상확률과정으로 변환하기 (Transforming non-stationary to stationary process) ]

 

 

먼저 예제로 사용할 약 판매량 (drug sales) 시계열 데이터를 가져와서 pandas DataFrame으로 만들고, 시계열 그래프를 그려보겠습니다. 

 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 

## getting drug sales dataset 
file_path = 'https://raw.githubusercontent.com/selva86/datasets/master/a10.csv' 
df = pd.read_csv(file_path, 
                 parse_dates=['date'], 
                 index_col='date') 

df.head(12)

#               value
# date	
# 	            3.526591
# 1991-08-01	3.180891
# 1991-09-01	3.252221
# 1991-10-01	3.611003
# 1991-11-01	3.565869
# 1991-12-01	4.306371
# 1992-01-01	5.088335
# 1992-02-01	2.814520
# 1992-03-01	2.985811
# 1992-04-01	3.204780
# 1992-05-01	3.127578
# 1992-06-01	3.270523

## time series plot 
df.plot(figsize=[12, 8]) 
plt.title('Non-Stationary Process: Increasing Variance + Trend + Seasonality', 
          fontsize=16) 
plt.show()

non-stationary process

 

위의 시계열 그래프에서 볼 수 있는 것처럼, (a) 분산이 시간의 흐름에 따라 증가 하고 (분산이 고정이 아님), (b) 추세(trend)가 있으며, (c) 1년 주기의 계절성(seasonality)이 있으므로, 비정상확률과정(non-stationary process)입니다. 

 

KPSS 검정을 통해서 확인해봐도 p-value가 0.01 이므로 유의수준 5% 하에서 귀무가설 (H0: 정상 시계열이다)을 기각하고, 대립가설(H1: 정상 시계열이 아니다)을 채택합니다. 

 

## UDF for KPSS test 
from statsmodels.tsa.stattools import kpss 
import pandas as pd 

def kpss_test(timeseries): 
    print("Results of KPSS Test:") 
    kpsstest = kpss(timeseries, regression="c", nlags="auto") 
    kpss_output = pd.Series(
        kpsstest[0:3], index=["Test Statistic", "p-value", "Lags Used"] ) 
    
    for key, value in kpsstest[3].items(): 
        kpss_output["Critical Value (%s)" % key] = value 
    print(kpss_output)
    
    
## 귀무가설 (H0): 정상 시계열이다
## 대립가설 (H1): 정상 시계열이 아니다 <-- p-value 0.01

kpss_test(df)

# Results of KPSS Test:
# Test Statistic           2.013126
# p-value                  0.010000
# Lags Used                9.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

(1) 분산이 고정적이지 않은 경우 분산 안정화 변환 (variance stabilizing transformation, VST)

 

분산이 고정적이지 않은 경우 멱 변환(Power Transformation)을 통해서 분산을 안정화(variance stabilization) 시켜줍니다. 분산이 고정적이지 않고 추세가 있는 경우 분산 안정화를 추세 제거보다 먼저 해줍니다. 왜냐하면 추세를 제거하기 위해 차분(differencing)을 해줄 때 음수(-)가 생길 수 있기 때문입니다. 

 

power transformation

 

 

원래의 시계열 데이터의 분산 형태에 따라서 적합한 멱 변환(power transformation)을 선택해서 정상확률과정으로 변환해줄 수 있습니다. 아래의 예제 시도표를 참고하세요. 

variance stabilizing transformation (power transfortion)

 

이번 포스팅에서 사용하는 예제는 시간이 흐릴수록 분산이 점점 커지는 형태를 띠고 있으므로 로그 변환(log transformation) 이나 제곱근 변환 (root transformation) 을 해주면 정상 시계열로 변환이 되겠네요. 아래 코드에서는 자연로그를 취해서 로그 변환을 해주었습니다. 

 

## Variance Stabilizing Transformation (VST) by Taking Logarithm
df_vst = np.log(df.value)

df_vst.head()

# date
# 1991-07-01    1.260332
# 1991-08-01    1.157161
# 1991-09-01    1.179338
# 1991-10-01    1.283986
# 1991-11-01    1.271408
# Name: value, dtype: float64


## plotting
df_vst.plot(figsize=(12, 8))
plt.title("Variance Stabilizing Transformation by taking Logarithm", 
          fontsize=16)
plt.show()

variance stabilizing transformation (VST)

 

 

위의 시도표를 보면 시간이 경과해도 분산이 안정화되었음을 알 수 있습니다.  KPSS 검정을 한번 더 해주면 아직 추세(trend)와 계절성(seasonality)가 남아있으므로 여전히 비정상확률과정을 따른다고 나옵니다. 

 

## 귀무가설 (H0): 정상 시계열이다  
## 대립가설 (H1): 정상 시계열이 아니다  <-- p-value 0.01 

kpss_test(df_vst)

# Results of KPSS Test:
# Test Statistic           2.118189
# p-value                  0.010000
# Lags Used                9.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

(2) 추세가 있는 경우 차분을 통한 추세 제거 (de-trend by differencing)

 

차분(differencing)은 현재의 시계열 값에서 시차 t 만큼의 이전 값을 빼주는 것입니다. 

 

1차 차분 = Delta1_Z(t) = Z(t) - Z(t-1) 

2차 차분 = Delta2_Z(t) = Z(t) - Z(t-1) - (Z(t-1) - Z(t-2)) =  Z(t) - 2Z(t-1) + Z(t-2)

 

Python의 diff() 메소드를 사용해서 차분을 해줄 수 있습니다. 이때 차분의 차수 만큼 결측값이 생기는 데요, dropna() 메소드를 사용해서 결측값은 제거해주었습니다. 

 

## De-trend by Differencing
df_vst_diff1 = df_vst.diff(1).dropna()

df_vst_diff1.plot(figsize=(12, 8))
plt.title("De-trend by 1st order Differencing", fontsize=16)
plt.show()

de-trending by 1st order differencing

 

위의 시도표를 보면 1차 차분(1st order differencing)을 통해서 이제 추세(trend)도 제거되었음을 알 수 있습니다. 하지만 아직 계절성(seasonality)이 남아있어서 정상성 조건은 만족하지 않겠네요. 그런데 아래에 KPSS 검정을 해보니 p-value가 0.10 으로서 유의수준 5% 하에서 정상성을 만족한다고 나왔네요. ^^;

 

## 귀무가설 (H0): 정상 시계열이다  <-- p-value 0.10
## 대립가설 (H1): 정상 시계열이 아니다 

kpss_test(df_vst_diff1)

# Results of KPSS Test:
# Test Statistic            0.121364
# p-value                   0.100000
# Lags Used                37.000000
# Critical Value (10%)      0.347000
# Critical Value (5%)       0.463000
# Critical Value (2.5%)     0.574000
# Critical Value (1%)       0.739000
# dtype: float64

 

 

 

(3) 계절성이 있는 경우 계절 차분을 통한 계절성 제거 (de-seasonality by seaanl differencing)

 

아직 남아있는 계절성(seasonality)을 계절 차분(seasonal differencing)을 사용해서 제거해보겠습니다. 1년 12개월 주기의 계절성을 띠고 있으므로 diff(12) 함수로 계절 차분을 실시하고, 12개의 결측값이 생기는데요 dropna() 로 결측값은 제거해주었습니다. 

 

## Stationary Process: De-seasonality by Seasonal Differencing
df_vst_diff1_diff12 = df_vst_diff1.diff(12).dropna()

## plotting
df_vst_diff1_diff12.plot(figsize=(12, 8))
plt.title("De-seasonality by Seasonal Differencing", 
          fontsize=16)
plt.show()

de-seasonality by seasonal differencing

 

위의 시도표를 보면 이제 계절성도 제거가 되어서 정상 시계열처럼 보이네요. 아래에 KPSS 검정을 해보니 p-value 가 0.10 으로서, 유의수준 5% 하에서 귀무가설(H0: 정상 시계열이다)을 채택할 수 있겠네요.

 

## 귀무가설 (H0): 정상 시계열이다  <-- p-value 0.10
## 대립가설 (H1): 정상 시계열이 아니다 

kpss_test(df_vst_diff1_diff12)

# Results of KPSS Test:
# Test Statistic           0.08535
# p-value                  0.10000
# Lags Used                8.00000
# Critical Value (10%)     0.34700
# Critical Value (5%)      0.46300
# Critical Value (2.5%)    0.57400
# Critical Value (1%)      0.73900
# dtype: float64

 

이제 비정상 시계열(non-stationary process)이었던 원래 데이터를 (1) log transformation을 통한 분산 안정화, (2) 차분(differencing)을 통한 추세 제거, (3) 계절 차분(seasonal differencing)을 통한 계절성 제거를 모두 마쳐서 정상 시계열(stationary process) 로 변환을 마쳤으므로, ARIMA 통계 모형을 적합할 수 있게 되었습니다. 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 백색잡음과정(White noise process), 확률보행과정(Random walk process), 정상확률과정(Stationary process)에 대해서 소개하였습니다. (https://rfriend.tistory.com/691

 

지난번 포스팅에서 특히 ARIMA와 같은 시계열 통계 분석 모형이 정상확률과정(Stationary Process)을 가정한다고 했습니다. 따라서 시계열 통계 모형을 적합하기 전에 분석 대상이 되는 시계열이 정상성 가정(Stationarity assumption)을 만족하는지 확인을 해야 합니다. 

 

[ 정상확률과정 (stationary process) 정의 ]

1. 평균이 일정하다. 
2. 분산이 존재하며, 상수이다. 
3. 두 시점 사이의 자기공분산은 시차(時差, time lag)에만 의존한다. 

 

정상 시계열 (stationary time series) 여부를 확인하는 방법에는  3가지가 있습니다. 

  [1] 시계열 그래프 (time series plot)

  [2] 통계적 가설 검정 (statistical hypothesis test)

  [3] 자기상관함수(ACF), 편자기상관함수(PACF)

 

정상성 확인하는 방법 (how to check the stationarity)

 

이번 포스팅에서는  이중에서 통계적 가설 검정 (Statistical Hypothesis Test) 을 이용해 정상성(stationarity) 여부를 확인하는 방법을 소개하겠습니다. 

 

- (a) Augmented Dickey-Fuller (“ADF”) test

- (b) Kwiatkowski-Phillips-Schmidt-Shin (“KPSS”) test

 

 

ADF 검정은 시계열에 단위근(unit root)이 존재하는지의 여부를 검정함으로써 정상 시계열인지 여부를 판단합니다. 단위근이 존재하면 정상 시계열이 아닙니다. 

 

단위근(unit root)이란 확률론의 데이터 검정에서 쓰이는 개념입니다. 주로 ‘단위근 검정’의 형식으로 등장합니다. 일반적으로 시계열 데이터는 시간에 따라 일정한 규칙을 가짐을 가정합니다. 따라서 매우 복잡한 형태를 갖는 시계열 데이터라도 다음과 같은 식으로 어떻게든 단순화시킬 수 있을 것이라 생각해볼 수 있습니다.

즉, ‘t시점의 확률변수는 t-1, t-2 시점의 확률변수와 관계를 가지면서 거기에 에러가 포함된 것’이라는 의미입니다. 여기서 편의를 위해 y0=0이라 가정합니다.  이제 아래의 방정식을 볼까요.

여기서 m=1이 위 식의 근이 된다면 이때의 시계열 과정을 단위근을 가진다고 합니다. 단위근 모형은 주로 복잡한 시계열 데이터를 단순하게나마 계산하려 할 때 사용됩니다.

 

 

Python으로 가상의 시계열 데이터셋(1개의 정상 시계열, 3개의 비정상 시계열)을 만들어서 위의 ADF test, KPSS test 를 각각 해보겠습니다. 위의 정상확률과정의 정의에 따라서, 추세(trend)가 있거나, 분산(variance)이 시간의 흐름에 따라 변하거나(증가 또는 감소), 계절성(seasonality)을 가지고 있는 시계열은 정상성(stationarity) 가정을 만족하지 않습니다. 

 

- (1) 정상 시계열 데이터 (stationary time series)

- (2) 추세를 가진 비정상 시계열 데이터 (non-stationary time series with trend)

- (3) 분산이 변하는 비정상 시계열 데이터 (non-stationary time series with changing variance)

- (4) 계절성을 가지는 비정상 시계열 데이터 (non-stationary time series with seasonality)

 

 

이제 (1)~(4)번 데이터별로 (a) ADF test, (b) KPSS test 로 정상성 여부를 차례대로 가설 검정해보겠습니다. 

 

 

(1) 정상 시계열 데이터 (stationary time series)

 

먼저, 자기회귀과정(auto-regressive process)을 따르는 AR(1) 과정의 시계열 데이터를 임의로 만들어보겠습니다. 

 

import numpy as np
import matplotlib.pyplot as plt

## Stationary Process 
## exmaple: AR(1) process 
np.random.seed(1) # for reproducibility
z_0 = 0 
rho = 0.6 # correlation b/w z(t) and z(t-1) 
z_all = [z_0] 

for i in range(200): 
    z_t = rho*z_all[i] + np.random.normal(0, 0.1, 1) 
    z_all.append(z_t) 
    
## plotting 
plt.rcParams['figure.figsize'] = (10, 6)
plt.plot(z_all)
plt.title("stationary time series : AR(1) process", fontsize=16)

## adding horizonal line at mean position 0 
plt.axhline(0, 0, 200, color='red', linestyle='--', linewidth=2) 
plt.show()

stationary process time series, 정상 시계열 

 

 

 

(1-a) ADF 검정 (Augmented Dickey-Fuller test)

 

Python의 statsmodels 모듈에 있는 adfuller 메소드를 사용해서 ADF 검정을 위한 사용자 정의함수를 정의해보겠습니다. 

 

#! pip install statsmodels

## UDF for ADF test
from statsmodels.tsa.stattools import adfuller
import pandas as pd

def adf_test(timeseries):
    print("Results of Dickey-Fuller Test:")
    dftest = adfuller(timeseries, autolag="AIC")
    dfoutput = pd.Series(
        dftest[0:4],
        index=[
            "Test Statistic",
            "p-value",
            "#Lags Used",
            "Number of Observations Used",
        ],
    )
    for key, value in dftest[4].items():
        dfoutput["Critical Value (%s)" % key] = value
    print(dfoutput)

 

 

이제 위에서 만든 ADF 검정 사용자정의함수 adf_test() 를 사용해서 정상시계열 z_all 에 대해서 ADF 검정을 해보겠습니다. p-value 가 8.74e-11 이므로 유의수준 5% 하에서 귀무가설 (H0: 단위근(unit root)이 존재한다. 즉, 정상 시계열이 아니다)을 기각하고 대립가설 (H1: 단위근이 없다. 즉, 정상 시계열이다.) 을 채택합니다. (맞음 ^_^)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root.(not stationary.)
# -- Alternate Hypothesis: The series has no unit root.(stationary.)

adf_test(z_all)  # p-value 8.740428e-11 => stationary

# Results of Dickey-Fuller Test:
# Test Statistic                -7.375580e+00
# p-value                        8.740428e-11
# #Lags Used                     1.000000e+00
# Number of Observations Used    1.990000e+02
# Critical Value (1%)           -3.463645e+00
# Critical Value (5%)           -2.876176e+00
# Critical Value (10%)          -2.574572e+00
# dtype: float64

 

 

이번에는 Python의 statsmodels 모듈을 사용해서 KPSS 검정 (Kwiatowski-Phillips-Schmidt-Shin Test) 을 위한 사용자 정의함수를 정의해보겠습니다. 

 

## UDF for KPSS test
from statsmodels.tsa.stattools import kpss
import pandas as pd

def kpss_test(timeseries):
    print("Results of KPSS Test:")
    kpsstest = kpss(timeseries, regression="c", nlags="auto")
    kpss_output = pd.Series(
        kpsstest[0:3], index=["Test Statistic", "p-value", "Lags Used"]
    )
    for key, value in kpsstest[3].items():
        kpss_output["Critical Value (%s)" % key] = value
    print(kpss_output)

 

 

이제 정상 시계열 z_all 에 대해서 위에서 정의한 kpss_test() 함수를 사용해서 정상성 여부를 확인해보겠습니다. p-value 가 0.065 로서 유의수준 10% 하에서 귀무가설 (H0: 정상 시계열이다) 를 채택하고, 대립가설 (H1: 정상 시계열이 아니다) 를 기각합니다. (유의수준 5% 하에서는 대립가설 채택). (맞음 ^_^)

* 이때, KPSS 검정은 ADF 검정과는 귀무가설과 대립가설이 정반대이므로 해석에 조심하시기 바랍니다.  

 

kpss_test(z_all) # p-value 0.065035 => stationary at significance level 10%

# Results of KPSS Test:
# Test Statistic           0.428118
# p-value                  0.065035
# Lags Used                5.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

(2) 추세를 가진 비정상 시계열 데이터 (non-stationary time series with trend)

 

다음으로, 추세(trend)를 가지는 가상의 비정상(non-stationary) 시계열 데이터를 만들어보겠습니다. 평균이 일정하지 않으므로 비정상 시계열이 되겠습니다. 

 

## time series with trend
np.random.seed(1) # for reproducibility
ts_trend = 0.01*np.arange(200) + np.random.normal(0, 0.2, 200)

## plotting
plt.plot(ts_trend)
plt.title("time series with trend", fontsize=16)

plt.show()

time series with trend (non-stationary process)

 

 

(2-a) ADF 검정 (ADF test)

 

위의 추세를 가지는 비정상 시계열 데이터에 대해 ADF 검정 (ADF test)를 실시해보면, p-value가 0.96 으로서 유의수준 5% 하에서 귀무가설 (H0: 시계열이 단위근을 가진다. 즉, 정상 시계열이 아니다) 을 채택하고, 귀무가설 (H1: 시계열이 단위근을 가지지 않는다. 즉, 정상 시계열이다.) 을 기각합니다. (맞음 ^_^)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root.(not stationary.)
# -- Alternate Hypothesis: The series has no unit root.(stationary.)

adf_test(ts_trend) # p-value 0.96 => non-stationary

# Results of Dickey-Fuller Test:
# Test Statistic                   0.125812
# p-value                          0.967780
# #Lags Used                      10.000000
# Number of Observations Used    189.000000
# Critical Value (1%)             -3.465431
# Critical Value (5%)             -2.876957
# Critical Value (10%)            -2.574988
# dtype: float64

 

 

 

(2-b) KPSS 검정 (KPSS test)

 

추세(trend)를 가지는 시계열에 대해 KPSS 검정(KPSS test)을 실시해보면, p-value 가 0.01 로서 귀무가설 (H0: 정상 시계열이다) 를 기각하고, 대립가설 (H1: 정상 시계열이 아니다) 를 채택합니다. (맞음 ^_^)

 

## KPSS test for checking the stationarity of a time series. 
## The null and alternate hypothesis for the KPSS test are opposite that of the ADF test.
# -- Null Hypothesis: The process is trend stationary.
# -- Alternate Hypothesis: The series has a unit root (series is not stationary).

kpss_test(ts_trend) # p-valie 0.01 => non-stationary

# Results of KPSS Test:
# Test Statistic           2.082141
# p-value                  0.010000
# Lags Used                9.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

(3) 분산이 변하는 비정상 시계열 데이터 (non-stationary time series with changing variance)

 

다음으로, 추세는 없지만 시간이 흐름에 따라서 분산이 점점 커지는 가상의 비정상 시계열을 만들어보겠습니다. 분산이 일정하지 않기 때문에 정상 시계열이 아닙니다. 

 

## time series with increasing variance
np.random.seed(1) # for reproducibility
ts_variance = []

for i in range(200):
    ts_new = np.random.normal(0, 0.001*i, 200).astype(np.float32)[0]
    ts_variance.append(ts_new)

## plotting    
plt.plot(ts_variance)
plt.title("time series with increasing variance", fontsize=16)

plt.show()

time series with increasing variance (non-stationary process)

 

 

(3-a) ADF 검정 (ADF test)

 

위의 시간이 흐름에 따라 분산이 커지는 비정상 시계열에 대해 ADF 검정(ADF test)을 실시하면, p-value가 5.07e-19 로서 유의수준 5% 하에서 귀무가설(H0: 단위근을 가진다. 즉, 정상 시계열이 아니다)를 기각하고, 대립가설(H1: 단위근을 가지지 않는다. 즉, 정상 시계열이다)를 채택합니다. (땡~! 틀림 -_-;;;)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root. (not stationary)
# -- Alternate Hypothesis: The series has no unit root. (stationary)

adf_test(ts_variance) # p-vaue 5.07e-19 => stationary (Opps, wrong result -_-;;;)

# Results of Dickey-Fuller Test:
# Test Statistic                -1.063582e+01
# p-value                        5.075820e-19
# #Lags Used                     0.000000e+00
# Number of Observations Used    1.990000e+02
# Critical Value (1%)           -3.463645e+00
# Critical Value (5%)           -2.876176e+00
# Critical Value (10%)          -2.574572e+00
# dtype: float64

 

 

(3-b) KPSS 검정 (KPSS test)

 

위의 시간이 흐름에 따라 분산이 커지는 비정상 시계열에 대해 KPSS 검정(KPSS test)을 실시하면, p-value가 0.035 로서 유의수준 5% 하에서 귀무가설(H0: 정상 시계열이다)를 기각하고, 대립가설(H1: 정상 시계열이 아니다)를 채택합니다. (딩동댕! 맞음 ^_^)  (ADF test 는 틀렸고, KPSS test 는 맞았어요.)

 

## KPSS test for checking the stationarity of a time series. 
## The null and alternate hypothesis for the KPSS test are opposite that of the ADF test.
# -- Null Hypothesis: The process is trend stationary.
# -- Alternate Hypothesis: The series has a unit root (series is not stationary).

kpss_test(ts_variance) # p-value 0.035 => not stationary

# Results of KPSS Test:
# Test Statistic           0.52605
# p-value                  0.03580
# Lags Used                3.00000
# Critical Value (10%)     0.34700
# Critical Value (5%)      0.46300
# Critical Value (2.5%)    0.57400
# Critical Value (1%)      0.73900
# dtype: float64

 

 

 

(4) 계절성을 가지는 비정상 시계열 데이터 (non-stationary time series with seasonality)

 

마지막으로, 코사인 주기의 계절성(seasonality)을 가지는 가상의 비정상(non-stationary) 시계열 데이터를 만들어보겠습니다. 

 

## time series with seasonality
## generating x from 0 to 4pi and y using numpy
np.random.seed(1) # for reproducibility
x = np.arange(0, 50, 0.2) # start, stop, step 
ts_seasonal = np.cos(x) + np.random.normal(0, 0.2, 250)

## ploting 
plt.plot(ts_seasonal) 
plt.title("time series with seasonality", fontsize=16)

plt.show()

time series with seasonality (non-stationary process)

 

 

(4-a) ADF 검정 (ADF test)

 

위의 계절성(seasonality)을 가지는 비정상 시계열에 대해서, ADF 검정(ADF test)을 실시하면, p-value가 3.14e-16 로서 유의수준 5% 하에서 귀무가설(H0: 단위근을 가진다. 즉, 정상 시계열이 아니다)를 기각하고, 대립가설(H1: 단위근을 가지지 않는다. 즉, 정상 시계열이다)를 채택합니다. (땡~! 틀림 -_-;;;)

 

## checking the stationarity of a series using the ADF(Augmented Dickey-Fuller) Test
# -- Null Hypothesis: The series has a unit root.(not stationary.)
# -- Alternate Hypothesis: The series has no unit root.(stationary.)

adf_test(ts_seasonal) # p-value 3.142783e-16 => stationary (Wrong result. -_-;;;)

# Results of Dickey-Fuller Test:
# Test Statistic                -9.516720e+00
# p-value                        3.142783e-16
# #Lags Used                     1.600000e+01
# Number of Observations Used    2.330000e+02
# Critical Value (1%)           -3.458731e+00
# Critical Value (5%)           -2.874026e+00
# Critical Value (10%)          -2.573424e+00
# dtype: float64

 

 

 

(4-b) KPSS 검정 (KPSS test)

 

위의 계절성(seasonality)을 가지는 비정상 시계열에 대해서, KPSS 검정(KPSS test)을 실시하면, p-value가 0.1 로서 유의수준 10% 하에서 귀무가설(H0: 정상 시계열이다)를 기각하고, 대립가설(H1: 정상 시계열이 아니다)를 채택합니다. (딩동댕! 맞음 ^_^)  (ADF test 는 틀렸고, KPSS test 는 맞았어요. 유의수준 5% 하에서는 둘 다 틀림.)

 

## KPSS test for checking the stationarity of a time series. 
## The null and alternate hypothesis for the KPSS test are opposite that of the ADF test.
# -- Null Hypothesis: The process is trend stationary.
# -- Alternate Hypothesis: The series has a unit root (series is not stationary).

kpss_test(ts_seasonal) # p-value 0.1 => non-stationary (at 10% significance level)

# Results of KPSS Test:
# Test Statistic           0.016014
# p-value                  0.100000
# Lags Used                9.000000
# Critical Value (10%)     0.347000
# Critical Value (5%)      0.463000
# Critical Value (2.5%)    0.574000
# Critical Value (1%)      0.739000
# dtype: float64

 

 

 

이상의 정상성 여부에 대한 통계적 가설 검정 결과를 보면,

(1) 정상 시계열에 대해 ADF test, KPSS test 모두 정상 시계열로 가설 검정 (모두 맞음 ^_^)

(2) 추세가 있는 시계열에 대해 ADF test, KPSS test 가 모두 정상 시계열이 아니라고 정확하게 가설 검정 (모두 맞음 ^_^)

(3) 분산이 변하는 시계열에 대해 ADF test 는 정상 시계열로 가설 검정 (틀렸음! -_-;;;), KPSS test 는 정상 시계열이 아니라고 가설 검정 (맞음 ^_^)

(4) 계절성이 있는 시계열에 대해 ADF test 는 정상 시계열로 가설 검정 (틀렸음! -_-;;;), KPSS test 는 정상 시계열이 아니라고 가설 검정 (맞음 ^_^)

 

ADF test 는 추세가 있는 비정상 시계열에 대해서는 정상 시계열이 아님을 잘 검정하지만, 분산이 변하거나 계절성이 있는 시계열에 대해서는 정상성 여부를 제대로 검정해내지 못하고 있습니다.

 

반면에 KPSS test 는 위의 4가지의 모든 경우에 정상성 여부를 잘 검정해내고 있습니다. 

 

통계적 가설 검정 외에 시계열 도표 (time series plot)을 꼭 그려보고 눈으로도 반드시 확인해보는 습관을 들이시기 바랍니다. 

 

 

[ Reference ]

- ADF test: https://en.wikipedia.org/wiki/Augmented_Dickey%E2%80%93Fuller_test 

- KPSS test: https://en.wikipedia.org/wiki/KPSS_test

- ADF test and KPSS test in Python: https://www.statsmodels.org/dev/examples/notebooks/generated/stationarity_detrending_adf_kpss.html

- 단위근 (unit root) : https://www.scienceall.com/%EB%8B%A8%EC%9C%84%EA%B7%BCunit-root-%E5%96%AE%E4%BD%8D%E6%A0%B9/

 

다음번 포스팅에서는 정상성 가정을 만족하지 않는 시계열 데이터에 대해 데이터 변환을 통해 정상화 시키는 방법을 소개하겠습니다. 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 Python을 사용해서 웹사이트에서 압축파일을 다운로드해서 압축을 해제하고 데이터셋을 합치는 방법을 소개하겠습니다. 

 

세부 절차 및 이용한 Python 모듈과 메소드는 아래와 같습니다. 

 

(1) os 모듈로 다운로드한 파일을 저장할 디렉토리가 없을 경우 새로운 디렉토리 생성하기

(2) urllib.request.urlopen() 메소드로 웹사이트를 열기 

(3) tarfile.open().extractall() 메소드로 압축 파일을 열고, 모든 멤버들을 압축해제하기

(4) pandas.read_csv() 메소드로 파일을 읽어서 DataFrame으로 만들기

(5) pandas.concat() 메소드로 모든 DataFrame을 하나의 DataFrame으로 합치기

(6) pandas.to_csv() 메소드로 합쳐진 csv 파일을 내보내기

 

 

먼저, 위의 6개 절차를 download_and_merge_csv() 라는 이름의 사용자 정의함수로 정의해보겠습니다.  

 

import os
import glob
import pandas as pd
import tarfile
import urllib.request

## downloads a zipped tar file (.tar.gz) that contains several CSV files, 
## from a public website. 
def download_and_merge_csv(url: str, down_dir: str, output_csv: str):
    """
    - url: url address from which you want to download a compressed file
    - down_dir: directory to which you want to download a compressed file
    - output_csv: a file name of a exported DataFrame using pd.to_csv() method
    """
    
    # if down_dir does not exists, then create a new directory
    down_dir = 'downloaded_data'
    if os.path.isdir(down_dir):
        pass
    else:
        os.mkdir(down_dir)
        
    # Open for reading with gzip compression.
    # Extract all members from the archive to the current working directory or directory path. 
    with urllib.request.urlopen(url) as res:
        tarfile.open(fileobj=res, mode="r|gz").extractall(down_dir)
    
    # concatenate all extracted csv files
    df = pd.concat(
        [pd.read_csv(csv_file, header=None) 
         for csv_file in glob.glob(os.path.join(down_dir, '*.csv'))])
    
    # export a DataFrame to a csv file
    df.to_csv(output_csv, index=False, header=False)

 

참고로, tarfile.open(fileobj, mode="r") 에서 4개의 mode 를 지원합니다. 

tarfile(mode) 옵션
-. mode="r": 존재하는 데이터 보관소로부터 읽기 (read)
-. mode="a": 존재하는 파일에 데이터를 덧붙이기 (append)
-. mode="w": 존재하는 파일을 덮어쓰기해서 새로운 파일 만들기 (write, create a new file overwriting an existing one)
-. mode="x": 기존 파일이 존재하지 않을 경우에만 새로운 파일을 만들기 (create a new file only if it does not already exist)

* for more information on tarfile module: https://docs.python.org/3/library/tarfile.html

 

 

현재 Jupyter Notebook 커널의 디렉토리에는 아래처럼  아직 다운로드한 파일이 없습니다. 

 

jovyan@kubecon-tutorial-0:~$ pwd
/home/jovyan
jovyan@kubecon-tutorial-0:~$ 
jovyan@kubecon-tutorial-0:~$ ls
data  down_merge_csv.ipynb  kale.log  lost+found
jovyan@kubecon-tutorial-0:~$ 
jovyan@kubecon-tutorial-0:~$

 

 

 

 

이제 위에서 정의한 download_and_merge_csv() 를 사용해서 

  (a) url='https://storage.googleapis.com/ml-pipeline-playground/iris-csv-files.tar.gz' 로 웹사이트로 부터 압축파일을 열고 모든 파일들을 해제해서 

  (b) down_dir='downloaded_data' 의 디렉토리에 다운로드하고, 

  (c) output_csv='iris_merged_data.csv'  라는 이름의 csv 파일로 모든 파일을 합쳐서 내보내기

를 해보겠습니다. 

 

download_and_merge_csv(
    url='https://storage.googleapis.com/ml-pipeline-playground/iris-csv-files.tar.gz', 
    down_dir='downloaded_data', 
    output_csv='iris_merged_data.csv')

 

 

아래의 화면캡쳐처럼 'iris_merged_data.csv' 라는 이름의 csv 파일이 새로 생겼습니다. 그리고 'downloaded_data' 라는 폴더도 새로 생겼습니다. 

 

 

 

 

터미널에서 새로 생긴 'downloaded_data' 로 디렉토리를 이동한 다음에, 파일 리스트를 확인해보니 'iris-1.csv', 'iris-2.csv', 'iris-3.csv' 의 3개 파일이 들어있네요. head 로 상위의 10 개 행을 읽어보니 iris 데이터셋이군요. 

 

jovyan@kubecon-tutorial-0:~$ ls
data  downloaded_data  down_merge_csv.ipynb  iris_merged_data.csv  kale.log  lost+found
jovyan@kubecon-tutorial-0:~$ 
jovyan@kubecon-tutorial-0:~$ 
jovyan@kubecon-tutorial-0:~$ cd downloaded_data/
jovyan@kubecon-tutorial-0:~/downloaded_data$ ls
iris-1.csv  iris-2.csv  iris-3.csv
jovyan@kubecon-tutorial-0:~/downloaded_data$ 
jovyan@kubecon-tutorial-0:~/downloaded_data$ 
jovyan@kubecon-tutorial-0:~/downloaded_data$ head iris-1.csv
5.1,3.5,1.4,0.2,setosa
4.9,3.0,1.4,0.2,setosa
4.7,3.2,1.3,0.2,setosa
4.6,3.1,1.5,0.2,setosa
5.0,3.6,1.4,0.2,setosa
5.4,3.9,1.7,0.4,setosa
4.6,3.4,1.4,0.3,setosa
5.0,3.4,1.5,0.2,setosa
4.4,2.9,1.4,0.2,setosa
4.9,3.1,1.5,0.1,setosa
jovyan@kubecon-tutorial-0:~/downloaded_data$

 

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요.  :-)

 

728x90
반응형
Posted by Rfriend
,