이번 포스팅에서는 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
,

이번 포스팅에서는 Python pandas 의 DataFrame에서 

 

(1) 그룹별로 x 칼럼을 정렬 후 누적 비율을 구한 후 (calculating the cumulative proportion by groups)  

(2) 그룹별로 특정 분위수의 위치 구하기 (getting the indices for a specific quantile p by groups)

 

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

 

그룹별로 연산을 수행하므로 pandas.DataFrame.groupby().apply(UDF) 형식으로 구문을 작성할 거예요. 

 

 

[ pandas DataFrame에서 그룹별로 정렬 후 누적 비율을 구한 후에 --> 그룹별로 특정 분위수 위치 구하기 ]

pandas getting cumulative proportion and quantile

 

 

 

먼저, 예제로 사용하기 위해 그룹('grp') 칼럼별 값('x')을 가지는 간단한 pandas DataFrame을 만들어보겠습니다. 

 

import numpy as np
import pandas as pd

df = pd.DataFrame({
    'grp': ['a', 'a', 'a', 'a', 'a', 'b', 'b', 'b', 'b', 'b'], 
    'x': [3, 1, 2, 7, 4, 4, 2, 5, 9, 7]})
    
print(df)
#   grp  x
# 0   a  3
# 1   a  1
# 2   a  2
# 3   a  7
# 4   a  4
# 5   b  4
# 6   b  2
# 7   b  5
# 8   b  9
# 9   b  7

 

 

 

(1) 그룹별로 x 칼럼을 정렬 후 누적 비율을 구한 후 (calculating the cumulative proportion by groups)  

 

그룹별로 x 칼럼에 대한 누적 비율을 구하기 위해, 먼저 그룹별로 x 칼럼의 비율(proportion)을 계산해서 'prop' 라는 칼럼을 추가해보겠습니다. x_prop() 라는 사용자 정의 함수를 정의한 후, df.groupby('grp').apply(x_prop) 처럼 그룹에 apply() 메소드로 사용자 정의 함수를 적용해서 연산을 했습니다. 

 

## adding the proportion column by group
def x_prop(group):
    group['prop'] = group.x / group.x.sum()
    return group

df = df.groupby('grp').apply(x_prop)


print(df)
#   grp  x      prop
# 0   a  3  0.176471
# 1   a  1  0.058824
# 2   a  2  0.117647
# 3   a  7  0.411765
# 4   a  4  0.235294
# 5   b  4  0.148148
# 6   b  2  0.074074
# 7   b  5  0.185185
# 8   b  9  0.333333
# 9   b  7  0.259259


## checking the sanity
df.groupby('grp').prop.sum()
#      grp
# a    1.0
# b    1.0
# Name: prop, dtype: float64

 

 

앞에서 계산한 그룹별 x 칼럼의 비율 'prop'을 그룹별로 내림차순(descending order)으로 정렬해서 보면 아래와 같습니다. 

 

## sorting in descending order by prop 
df.sort_values(by=['grp', 'prop'], ascending=False)

#     grp	x	prop
# 8	b	9	0.333333
# 9	b	7	0.259259
# 7	b	5	0.185185
# 5	b	4	0.148148
# 6	b	2	0.074074
# 3	a	7	0.411765
# 4	a	4	0.235294
# 0	a	3	0.176471
# 2	a	2	0.117647
# 1	a	1	0.058824

 

 

pandas 의 cumsum() 메소드를 사용해서 그룹별 x칼럼의 비율 'prop'의 누적 합계 (cumulative sum) 인 'cum_prop' 를 그룹별로 계산해보겠습니다. 역시 비율 'prop'에 대해서 누적 합계(cum_sum)를 구하는 사용자 정의 함수 'cum_prop()'를 먼저 정의한 후에, 이를 df.groupby('grp').apply(cum_prop) 처럼 apply() 메소드에 사용자 정의함수를 적용해서 계산했습니다. 

 

## sorting in descending order by prop and calculating the cumulative sum of prop
def cum_prop(group):
    group['cum_prop'] = group.sort_values(
        by='prop', ascending=False).prop.cumsum()
    return group

df = df.groupby('grp').apply(cum_prop)


df.sort_values(by=['grp', 'cum_prop'])

#     grp	x	prop	        cum_prop
# 3	a	7	0.411765	0.411765
# 4	a	4	0.235294	0.647059
# 0	a	3	0.176471	0.823529
# 2	a	2	0.117647	0.941176
# 1	a	1	0.058824	1.000000
# 8	b	9	0.333333	0.333333
# 9	b	7	0.259259	0.592593
# 7	b	5	0.185185	0.777778
# 5	b	4	0.148148	0.925926
# 6	b	2	0.074074	1.000000

 

 

 

위의 예시는 간단한 편이므로 아래처럼 사용자 정의 함수를 정의하는 대신에 apply() 메소드 안에 바로 lambda 로 'prop'에 대해서 내림차순 정렬 후 누적 합계를 구하는 함수를 바로 써줘도 됩니다. 

 

## or, equivalentsly, using lambda function for cumulative proportion
df.groupby('grp').apply(lambda x: x.sort_values(by='prop', ascending=False).prop.cumsum())

# grp   
# a    3    0.411765
#      4    0.647059
#      0    0.823529
#      2    0.941176
#      1    1.000000
# b    8    0.333333
#      9    0.592593
#      7    0.777778
#      5    0.925926
#      6    1.000000
# Name: prop, dtype: float64

 

 

 

(2) 그룹별로 특정 분위수의 위치 구하기 (getting the indices for a specific quantile p by groups)

 

이제 위에서 구한 '그룹별 비율의 누적 합계('cum_prop')'에 대해서 pandas.Series.searchsorted(values, side='left') 메소드를 사용해서 특정 비율이 들어갈 위치를 구해보겠습니다.

 

비율에 대해 내림차순 정렬 후 누적 합계를 구한 값에 대해 특정 값이 들어갈 위치를 구하는 것이므로, 결과적으로 자료 크기 순서에 따른 위치값인 분위수(quantile) 를 구할 수 있게 됩니다. 인덱스가 '0'부터 시작하므로 위치를 구하기 위해서 반환받는 값에 '+1' 을 해주었습니다. 

 

그룹별로 특정 분위수의 위치를 구하고 싶으므로, 분위수를 구하는 사용자 정의 함수인 'quantile_idx()' 함수를 정의한 후에, 이를 df.groupby('grp').apply(quantile_idx, p) 처럼 apply() 메소드에 사용자 정의 함수와 매개변수 p를 입력해서 적용해주면 되겠습니다.

 

그룹별로 분위수 p=0.2, p=0.5, p=0.7 인 곳의 위치를 구해보니 잘 작동하는군요. 

 

## pandas.Series.searchsorted
## Series.searchsorted(value, side='left', sorter=None)[source]
## Find indices where elements should be inserted to maintain order.

def quantile_idx(group, p=0.5):
    group = group.sort_values(by='cum_prop', ascending=True)
    return group.cum_prop.searchsorted(p) + 1
    

## getting the index of quantile p=0.2 by group
df.groupby('grp').apply(quantile_idx, p=0.2)

# grp
# a    1
# b    1
# dtype: int64


## getting the index of quantile p=0.5 by group
df.groupby('grp').apply(quantile_idx, p=0.5)

# grp
# a    2
# b    2
# dtype: int64


## getting the index of quantile p=0.7 by group
df.groupby('grp').apply(quantile_idx, p=0.7)

# grp
# a    3
# b    3
# dtype: int64

 

 

[Reference]

* pandas.Series.searchsorted() method:  https://pandas.pydata.org/docs/reference/api/pandas.Series.searchsorted.html

 

 

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

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

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 

 (1) 그룹별로 x 칼럼을 기준으로 내림차순 정렬 후 (sorting by x in ascending order)

 (2) 그룹별로 y 칼럼의 첫번째 값, 마지막 값을 DataFrame에 칼럼 추가하기

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

 

(방법 1) pandas.DataFrame 의 transform('first', 'last') 메소드를 사용하는 방법

(방법 2) 그룹별 y의 첫번째 값, 마지막 값을 구해 DataFrame을 만들고, merge() 메소드로 합치는 방법

 

pandas DataFrame sort_values() groupby() transform('first') transform('last')

 

 

먼저, 예제로 사용할 간단한 DataFrame을 만들어보겠습니다. 

 

import numpy as np
import pandas as pd

df = pd.DataFrame({
    'grp': ['a', 'a', 'a', 'b', 'b', 'b'], 
    'x': [2, 3, 1, 4, 6, 5], 
    'y': [10, 20, 30, 40, 50, 60]
}) 

df

#   grp	x	y
# 0	a	2	10
# 1	a	3	20
# 2	a	1	30
# 3	b	4	40
# 4	b	6	50
# 5	b	5	60

 

 

 

(방법 1) pandas.DataFrame 의 transform('first', 'last') 메소드를 사용하는 방법

 

그룹별로 'x' 칼럼을 기준으로 내림차순으로 정렬하려면 df.sort_values(by=['grp', 'x']) 메소드를 사용합니다. 

 

## sorting by 'grp' and 'x' in ascnding order
df.sort_values(by=['grp', 'x'])

#   grp	x	y
# 2	a	1	30
# 0	a	2	10
# 1	a	3	20
# 3	b	4	40
# 5	b	5	60
# 4	b	6	50

 

 

그러면, 이제 x를 기준으로 내림차순 정렬한 후에 'grp' 그룹별로 y 칼럼의 첫번째 값('first')과 마지막 값('last')을 가져다가 기존의 df DataFrame에 새로운 칼럼을 추가해 보겠습니다. groupby('grp') 메소드로 'grp' 그룹별 연산을 하게 되고, transform('first')는 첫번째 값을 가져다가 DataFrame에 칼럼을 추가하며, transform('last')는 마지막 값을 가져다가 DataFrame에 칼럼을 추가합니다. 

 

## adding columns of the first and last value of y by group
df['y_first'] = df.sort_values(by=['grp', 'x'])\
    .groupby('grp').y.transform('first')
df['y_last'] = df.sort_values(by=['grp', 'x'])\
    .groupby('grp').y.transform('last')


 df
# 	grp	x	y	y_first	y_last
# 0	a	2	10	30	20
# 1	a	3	20	30	20
# 2	a	1	30	30	20
# 3	b	4	40	40	50
# 4	b	6	50	40	50
# 5	b	5	60	40	50

 

 

 

(방법 2) 그룹별 y의 첫번째 값, 마지막 값을 구해 DataFrame을 만들고, merge() 메소드로 합치는 방법

 

두번째 방법은 그룹별로 x를 기준으로 정렬 후 그룹별로 y 값의 첫번째 값과 마지막 값을 구해서 별도의 DataFrame을 만든 후에, 이를 원래의 DataFrame에 merge() 하는 것입니다. DB의 테이블을 join 하는 것과 유사한 방식이예요. 

 

## creating a sample DataFrame with 2 groups
df = pd.DataFrame({
    'grp': ['a', 'a', 'a', 'b', 'b', 'b'], 
    'x': [2, 3, 1, 4, 6, 5], 
    'y': [10, 20, 30, 40, 50, 60]
}) 


## making a DataFrame with the first and last values of y by groups
y_first = df.sort_values(by='x').groupby('grp').y.first()
y_last = df.sort_values(by='x').groupby('grp').y.last()

df_grp_fst_lst = pd.DataFrame({
    'y_first': y_first, 
    'y_last': y_last
})

df_grp_fst_lst
#     y_first	y_last
# grp
# a	      30	  20
# b	      40	  50

 

 

 

pd.merge(DataFrame1, DataFrame2, how='left', on='key') 방식으로 key를 기준으로 Left Join 하면 되겠네요. 

 

## merging df_grp_fst_lst to df DataFrame by left join on 'grp'
df2 = pd.merge(df, df_grp_fst_lst, how='left', on='grp') 
# or, equivalently: df2= df.merge(df_grp_fst_lst, how='left', on='grp')


df2
# 	grp	x	y	y_first	y_last
# 0	a	2	10	30	20
# 1	a	3	20	30	20
# 2	a	1	30	30	20
# 3	b	4	40	40	50
# 4	b	6	50	40	50
# 5	b	5	60	40	50

 

 

* pandas DataFrame merge(): https://rfriend.tistory.com/258

* pandas DataFrame transform(): https://rfriend.tistory.com/403

 

 

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

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

 

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
,

이번 포스팅에서는 Python의 matplotlib 모듈을 사용해서 

 

  (1) 그래프에 수평선 추가하기 (adding horizontal lines)

  (2) 그래프에 수직선 추가하기 (adding vertical lines) 

 

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

 

예제로 사용할 샘플 데이터셋을 정규분포로부터 난수를 생성해서 100개 샘플을 추출하고, 점 그래프를 그려보겠습니다. 

이 기본 점 그래프에 수평선과 수직선을 차례대로 추가해보겠습니다. 

 

import numpy as np
import matplotlib.pyplot as plt

## generating random numbers
np.random.seed(1004)
x = np.random.normal(0, 1, 100)

## plotting the original data
plt.figure(figsize = (10, 6))
plt.plot(x, linestyle='none', marker='o', color='gray')
plt.show()

 

 

 (1) 그래프에 수평선 추가하기 (adding horizontal lines)

 

(a) plt.axhline(y, xmin, xmax) : 축을 따라서 수평선을 추가, xmin 과 xmax 는 0~1 사이의 값을 가짐

(b) plt.hlines(y, xmin, xmax) : xmin ~ xmax 까지 각 y 값의 수평선을 추가

(c) plt.plot((x1, x2), (y1, y2)) : (x1, x2), (y1, y2) 좌표를 연결하는 선 추가

 

(a) 번의 plt.axhline() 은 y축에서 부터 수평선이 시작하고, xmin~xmax 로 0~1 사이의 비율 값을 가지는 반면에, (b)번의 plt.hlines() 는 xmin 값 위치에서 부터 수평선이 시작하고, xmin~xmax 값으로 좌표값을 받는다는 차이점이 있습니다. 

(c) 번의 plt.plot() 은 단지 수평선, 수직선 뿐만이 아니라 범용적으로 두 좌표를 연결하는 선을 추가할 수 있습니다. 

 

plt.figure(figsize = (10, 6))
plt.plot(x, linestyle='none', marker='o', color='gray')
plt.title("Plot with Horizontal Lines", fontsize=16)

## (1) adding a horizontal line across the axis
## xmin and xmax should be b/w 0 and 1
plt.axhline(y=3, xmin=0, xmax=1, color='blue', linestyle='solid')
plt.axhline(y=2, xmin=0.1, xmax=0.9, color='blue', linestyle='dashed')

## (2) adding a horizontal line at each y from xmin to xmax
plt.hlines(y=0, xmin=0, xmax=50, color='red', linestyle='dotted')

## (3) adding a horizontal line using (x1, x2), (y1, y2) coordinates
plt.plot((50, 100), (-2, -2), color='black', linestyle='dashdot')

plt.show()

 

horizontal lines using matplotlib

 

 

(2) 그래프에 수직선 추가하기 (adding vertical lines) 

 

(a) plt.axvline(x, ymin, ymax) : 축을 따라서 수직선을 추가, ymin 과 ymax 는 0~1 사이의 값을 가짐

(b) plt.vlines(x, ymin, ymax) : ymin ~ ymax 까지 각 x 값의 수평선을 추가

(c) plt.plot((x1, x2), (y1, y2)) : (x1, x2), (y1, y2) 좌표를 연결하는 선 추가

 

(a) 번의 plt.axvline() 은 x축에서 부터 수평선이 시작하고, ymin~ymax 로 0~1 사이의 비율 값을 가지는 반면에, (b)번의 plt.vlines() 는 ymin 값 위치에서 부터 수평선이 시작하고, ymin~ymax 값으로 좌표값을 받는다는 차이점이 있습니다. 

(c) 번의 plt.plot() 은 단지 수평선, 수직선 뿐만이 아니라 범용적으로 두 좌표를 연결하는 선을 추가할 수 있습니다. 

 

plt.figure(figsize = (10, 6))
plt.plot(x, linestyle='none', marker='o', color='gray')
plt.title("Plot with vertical Lines", fontsize=16)

## (1) adding a vertical line across the axis
## ymin and ymax should be b/w 0 and 1
plt.axvline(x=0, ymin=0, ymax=1, color='blue', linestyle='solid')
plt.axvline(x=10, ymin=0.1, ymax=0.9, color='blue', linestyle='dashed')

## (2) adding a vertical line at each y from xmin to xmax
plt.vlines(x=50, ymin=0, ymax=3, color='red', linestyle='dotted')

## (3) adding a vertical line using (x1, x2), (y1, y2) coordinates
plt.plot((100, 100), (0, -3), color='black', linestyle='dashdot')

plt.show()

 

vertical lines using matplotlib

 

 

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

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

 

728x90
반응형
Posted by Rfriend
,