"Drift"와 "Shift"는 시계열 데이터 분석에서 중요한 두 가지 개념이며, 시간에 따른 데이터의 다양한 변화 또는 변동을 나타냅니다. 

 

이번 포스팅에서는 

 

1. 시계열 데이터의 Drift 란 무엇인가? 

2. 시계열 데이터의 Shift 란 무엇인가? 

3. 시계열 데이터의 Drift와 Shift 분석은 어떻게 하나?

4. Python 을 이용한 시계열 데이터 Drift, Shift 생성/ 시각화 및 탐지 예시

 

에 대해서 소개하겠습니다. 

 

 

Time Series - Drift, Shift

 


1. 시계열 데이터에서의 Drift 란 무엇인가?  

 

- 정의: Drift 는 시계열 데이터가 시간에 따라 느리고 점진적이며 종종 선형적으로 변화는 것(a slow, gradual, and often linear change in the time series data over time)을 의미합니다. 이는 한 방향으로 지속적인 변화 경향, 즉 상승 또는 하락하는 경향입니다. 

- 예시: 재무 데이터에서, 주식 가치가 몇 년 동안 지속적으로 상승하거나 하락하는 경우, 이는 드리프트의 예입니다. 

- 중요성: 드리프트를 이해하고 식별하는 것은 정확한 예측 및 장기적인 추세를 반영하는 모델을 만드는 데 중요합니다. 드리프트를 무시하면 시간이 지남에 따라 점점 더 부정확해지는 모델이 될 수 있습니다. 




2. 시계열 데이터에서의 Shift 란 무엇인가? 

 

- 정의: Shfit, 종종 "레벨 시프트(level shift)"로 언급되며, 시계열의 평균 또는 분산에서 갑작스러운 변화 (a sudden change in the mean or variance of a time series)를 의미합니다. 드리프트와 달리, 시프트는 갑작스럽게 발생하며 언제든지 일어날 수 있습니다.  

- 유형: 평균 시프트 (Mean Shift)와 분산 시프트 (Variance Shift) 로 나누어볼 수 있습니다. 

  -. 평균 시프트 (Mean Shift): 시리즈의 평균 수준이 갑자기 변경되는 경우입니다. 예를 들어, 성공적인 마케팅 캠페인으로 인해 판매량이 갑자기 증가하는 경우가 이에 해당합니다.  

  -. 분산 시프트 (Variance Shift): 시리즈의 변동성이 변경되는 경우입니다, 예를 들어 주식 가격의 변동성 증가와 같은 경우입니다. 

- 중요성: 시프트를 감지하는 것은 시계열 데이터에서 갑작스러운 변화를 이해하고 대응하는 데 필수적입니다. 시프트는 중대한 사건, 시스템의 구조적 변화 또는 데이터에 영향을 미치는 외부 요인의 변화를 나타낼 수 있습니다.

 

 

3. 시계열 데이터의 Drift와 Shift 분석은 어떻게 하나? 

 

시계열 데이터에서 드리프트와 시프트를 분석하는 것은 시간이 지남에 따라 데이터에서 점진적인 (Drift) 및 급격한 (Shift) 변화를 식별하고 정량화하는 것을 포함합니다. 이 목적을 위해 다양한 방법과 알고리즘이 사용되며, 각각은 데이터 유형과 분석 요구 사항에 적합합니다. 여기에 주요한 몇 가지를 소개합니다. 


3-1. 통계적 공정 관리 방법 (SPC, Statistical Process Control) 

- 제어 차트 (Control Charts): 평균 또는 분산의 변화를 탐지하기 위해 CUSUM (Cumulative Sum Control Chart) 및 EWMA(Exponentially Weighted Moving Average)와 같은 차트가 사용됩니다. 품질 관리에서 널리 사용됩니다. 

- Shewhart 제어 차트 (Shewhart Control Charts): 큰 시프트는 감지할 수 있지만 작은 시프트나 드리프트에는 덜 민감합니다. 


3-2. 시계열 분석 방법 (Time Series Analysis Methods)

- 이동 평균(Moving Average): 단기 변동을 평활화하고 장기 추세나 주기(long-term trends or cycles)를 강조하는 데 도움이 됩니다. 

- 시계열 분해 방법 (Time Series Decomposition Methods): STL (Seasonal and Trend decomposition using Loess)과 같은 기술은 시계열을 추세(Trend), 계절성(Seasonality) 및 잔차(Residual) 요소로 분해하여 드리프트와 시프트를 식별하는 데 도움이 됩니다. 


3-3. 변화점 감지 알고리즘 (Change Point Detection Algorithms)

- Bai-Perron 검정: 선형 모델에서 여러 중단점을 감지하는 데 사용됩니다.

- 변화점 분석 (Changepoint Analysis)PELT (Pruned Exact Linear Time), 이진 분할 (Binary Segmentation), 세그먼트 이웃 (Segment Neighborhoods)과 같은 방법은 대규모 데이터셋에서 다중 변화점을 감지하도록 설계되었습니다.

- CUSUM (Cumulative Sum Control Chart) 및 EWMA (Exponentially Weighted Moving Average): 평균 또는 분산의 변화를 식별하는 변화점 감지에 적응할 수 있습니다.


3-4. 기계 학습 접근 방식 

- 지도 학습 (Supervised Learning): 변화점 전후를 나타내는 레이블이 있는 데이터가 있을 경우, 변화점을 감지하기 위해 기계 학습 모델을 훈련할 수 있습니다. 

- 비지도 학습 (Unsupervised Learning): 클러스터링(예: k-means clustering, DBSCAN)과 같은 기술은 데이터 분포의 변화를 감지하여 시프트를 식별하는 데 사용될 수 있습니다. 


3-5. 신호 처리 기술 (Signal Processing Techniques)

- 푸리에 변환 (Fourier Transform): 주파수 영역의 시프트를 식별하는 데 유용합니다. 

- 웨이블릿 변환 (Wavelet Transform): 비정상 시계열에서 급격한 및 점진적인 변화를 모두 감지하는 데 효과적입니다. 


3-6. 베이지안 방법 (Bayesia Methods) 

- 베이지안 실시간 변화점 감지 (Bayesian Online Change Point Detection): 실시간 데이터 스트림에서 변화점을 감지하는 데 유용한 확률적 접근 방식입니다. 


3-7. 계량경제학에서의 구조적 분석 검정 

- Dickey-Fuller 검정: 시계열 샘플에서 단위근을 테스트하는 데 사용되며, 이는 종종 드리프트의 징후입니다. 

- Chow 검정: 특정 시점에서 구조적 중단이 있는지 여부를 결정하기 위해 설계되었습니다. 


3-8. 회귀 기반 방법 (Regression-based Methods)

- 선형 회귀 (Linear Regression): 기울기 계수 (slope coefficient)를 평가하여 드리프트를 탐지할 수 있습니다. 

- 분할 회귀 (Segmented Regression): 알려지지 않은 시점에서 회귀 모델의 변화를 탐지하는 데 유용합니다. 


적용 방법은 데이터의 특성 (노이즈 수준, 샘플 크기, 계절성 존재 여부 등)과 분석 요구 사항 (실시간 감지, 작은 변화에 대한 민감도 등)에 따라 다릅니다. 실제로 이러한 방법들의 조합이 종종 사용되어 결과를 검증하고 시계열 데이터에 대한 종합적인 이해를 얻습니다. 


요약하자면, Drift와 Shift 분석은 시계열 데이터에서의 변화가 점진적인지 또는 갑작스러운지를 이해하는 데 도움이 됩니다. 이러한 이해는 시계열 데이터를 기반으로 한 정확한 모델링, 예측 및 의사결정에 있어 중요합니다. 



4. Python 을 이용한 시계열 데이터 Drift, Shift 생성/ 시각화 및 탐지 예시 

 

4-1. Python을 이용한 시계열 데이터 Drift, Mean Shift, Variance Shift 생성/ 시각화

 

## Creating a time series with Drift, Mean Shift, Variance Shift
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Setting a seed for reproducibility
np.random.seed(1004)

# (1) Creating a time series with drift
time = np.arange(100)
drift = 0.5 * time  # Linear drift, Trend
noise = np.random.normal(0, 0.5, 100)
time_series_drift = drift + noise


# (2) Creating a time series with Mean Shift
shift_point = 50
time_series_mean_shift = np.random.normal(0, 0.5, 100)
time_series_mean_shift[shift_point:] += 10  # Adding a shift at the midpoint


# (3) Creating a time series with Variance Shift
first_variance = 0.5
second_variance = 5
time_series_first = np.random.normal(0, first_variance, 50)
time_series_second = np.random.normal(0, second_variance, 50)
time_series_var_shift = np.concatenate((time_series_first, time_series_second), axis=0)


# Plotting the time series
plt.figure(figsize=(8, 12))

plt.subplot(3, 1, 1)
plt.plot(time, time_series_drift)
plt.title("Time Series with Drift")

plt.subplot(3, 1, 2)
plt.plot(time, time_series_mean_shift)
plt.title("Time Series with Mean Shift")

plt.subplot(3, 1, 3)
plt.plot(time, time_series_var_shift)
plt.title("Time Series with Variance Shift")

plt.tight_layout()
plt.show()

 

Time Series - Drift, Mean Shift, Variance Shift

 

 

 

4-2. Python 으로 회귀모형 기반 Drift 분석

 

## Analysis of Time Series Drift using Linear Regression Model
from sklearn.linear_model import LinearRegression

# Reshape data for sklearn
X = time.reshape(-1, 1)
y = time_series_drift

# Linear regression
model = LinearRegression()
model.fit(X, y)

# Slope coefficients
print('---' * 10)
print(f'Slope Coefficient: {model.coef_[0]:.3f}')
print('---' * 10)

# Predicted trend line
trend_line = model.predict(X)

# Plotting
plt.figure(figsize=(6, 4))
plt.plot(time, time_series_drift, label='Time Series')
plt.plot(time, trend_line, label='Trend Line', color='red')
plt.title("Drift Detection")
plt.legend()
plt.show()

 

Detecting Time Series Drift using Linear Regression Method

 

 

 

4-3. Python으로 시계열 Shift 분석

 

- Mean Shift 탐지

 

# Detect Mean Shift
def detect_mean_shift(series, window=5):
    for i in range(window, len(series) - window):
        before = series[i-window:i]
        after = series[i:i+window]
        if np.abs(np.mean(after) - np.mean(before)) > 9:  # Threshold for shift
            return i
    return None

mean_shift_point_detected = detect_mean_shift(time_series_mean_shift)

print("----" * 10)
print(f"Mean Shift detected at point: {mean_shift_point_detected}")
print("----" * 10)

# ----------------------------------------
# Shift detected at point: 50
# ----------------------------------------

 

 

 

- PELT (Pruned Exact Linear Time), Binary Segment 알고리즘을 이용한 변곡점 탐지

 

##-- install "ruptures" module at terminal 
!pip install ruptures


## -- change point detection using PELT (Pruned Exact Linear Time) Algorithm
# Specify the PELT model and fit it to the data
model = "l2"  # Change to "rbf" for the radial basis function cost model
algo = rpt.Pelt(model=model).fit(time_series_mean_shift)

# Retrieve the change points
result = algo.predict(pen=10)  # Adjust the penalty value as needed

# Print the detected change points
print("Change points:", result)
# Change points: [50, 100]


## -- change point detection using Binary Segmentation Algorithm
# Specify the Binary Segment model and fit it to the data
model = "l2"  # Change to "rbf" for the radial basis function cost model
algo = rpt.Binseg(model=model).fit(time_series_mean_shift)

# Retrieve the change points
result = algo.predict(pen=10)  # Adjust the penalty value as needed

# Print the detected change points
print("Change points:", result)
# Change points: [50, 100]

 

 

 

- Variance Shift 탐지 

 

# Detect Variance Shift
def detect_var_shift(series, window=5):
    for i in range(window, len(series) - window):
        before = series[i-window:i]
        after = series[i:i+window]
        if np.abs(np.var(after) - np.var(before)) > 30:  # Threshold for shift
            return i + int(window/2) + 1
    return None

var_shift_point_detected = detect_var_shift(time_series_var_shift)

print("----" * 10)
print(f"Variance Shift detected at point: {var_shift_point_detected}")
print("----" * 10)

# ----------------------------------------
# Variance Shift detected at point: 50
# ----------------------------------------

 

 

 

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

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

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 베타분포 (Beta Distribution)란 무엇이고, 왜 중요한가에 대해 다뤄보겠습니다. 

 

 

1. 베타 분포 (Beta Distribution)는 무엇인가?

 

베타 분포는 구간 [0, 1] 에서 정의된 연속 확률 분포입니다. 일반적으로 Beta(α,β)로 표기되며, 여기서 α와 β는 양의 형태 매개변수입니다. 베타 분포의 확률 밀도 함수(PDF)는 다음과 같습니다:

베타 분포 확률밀도함수



여기서 B(α,β)는 베타 함수로, 확률 밀도 함수 아래의 면적이 1이 되도록 하는 정규화 상수 역할을 합니다.


베타 분포는 형태 매개변수의 값에 따라 다양한 모양을 가지며 다음과 같은 특징을 갖습니다:

  • α=β=1인 경우, 베타 분포는 구간 [0,1]에서 균일 분포입니다.
  • α와 β가 증가함에 따라 분포는 평균 주변에 더 피크한 형태를 가집니다.

 

 


2. 베이지안 통계(Bayesian Statistics)에서 베타분포(Beta Distribution)는 왜 중요한가? 

 


베타 분포는 베이지안 통계에서 여러 이유로 중요합니다:


(1) 켤레 사전 분포 (Conjugate Prior)


베타 분포는 이항, 음이항 이항, 베르누이 및 기하 분포에 대한 켤레 사전 분포입니다(The beta distribution is the conjugate prior for the binomial, negative binomial, Bernoulli, and geometric distributions). 이는 베이지안 분석에서 베타 분포를 사전 분포로 사용하고 관측된 데이터로 업데이트하면 사후 분포도 베타 분포가 되는 특성을 의미합니다. 이러한 특성은 베이지안 추론에서 업데이팅 과정을 단순화합니다.


(2) 사전 신념 표현 (Expressing Prior Beliefs)


베타 분포는 이진 결과(성공/실패, 앞면/뒷면 등)에 대한 성공 확률에 대한 사전 신념을 모델링하는 데 자주 사용됩니다. 형태 매개변수 α와 β(shape parameters α and β)에 적절한 값을 선택함으로써 연구자는 한 방향으로의 강한 신념부터 더 불확실하거나 흩어진 신념까지 다양한 사전 신념을 표현할 수 있습니다.


(3) 베이지안 업데이팅 (Bayesian Updating)


베타 분포는 베이지안 분석에서 사전 지식을 분석에 통합하는 데 도움이 됩니다. 새로운 데이터를 관측한 후에는 사전 분포를 사후 분포로 업데이트하며, 베타 분포의 수학적 속성은 이 업데이트 과정을 분석적으로 다루기 쉽게 만듭니다.

 


(4) 베타-이항 모델 (Beta-Binomial Model)

 

베타 분포는 이항 분포와 밀접한 관련이 있습니다. 베이지안 분석에서 사전으로 사용되어진다면 베르누이 분포와 함께 사용되는 경우가 많으며, 결과적인 모델은 베타-이항 모델 (beta-binomial model)이라고 불립니다. 이 모델은 베이지안 맥락에서 이항 확률 변수의 분포를 모델링하는 데 널리 사용됩니다. 

요약하면 베타 분포는 베이지안 통계에서 사전 신념을 모델링하고 이를 관측 데이터로 업데이트하는 데 유연하고 분석적으로 편리한 방법을 제공합니다. 그 켤레성(conjugate nature)은 베이지안 추론에서 수학적 계산을 간단하게 만듭니다.

 


3. Python을 이용해서 베타 분포의 파라미터별 그래프 그려보기

X 축 값의 구간은 [0, 1] 사이의 실수로 하고, 형태 매개변수 α와 β(shape parameters α and β) 를 변화시켜가면서 베타 분포를 그래프로 그려보겠습니다. 가령, 아래 그래프에서 보는 바와 같이,  α와 β의 값이 각각 1인 경우 확률값이 모두 0으로 동일합니다. (즉, 사전 정보가 없다는 의미입니다) 

α와 β의 값에 따라 매우 다양한 형태의 확률 분포를 표현할 수 있어서 유용합니다. 

 

# Plot the beta distribution for each set of alpha and beta values
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta

# Define a range of x values
x = np.linspace(0, 1, 1000)

# Define different alpha and beta values
d = [.5, 1, 2, 4]

fig = plt.figure(figsize=(16, 12))

for i in np.arange(16):
    a = d[i % 4]
    b = d[i // 4]
    ax = fig.add_subplot(4, 4, 1 + i)
    ax.set_title('Beta Distribution(a=' + str(a) + ', b=' + str(b) + ')')
    ax.set_xlabel('X')
    ax.set_ylabel('P(X)')
    ax.grid()
    ax.plot(x, beta(a, b).pdf(x), 'b--')

plt.subplots_adjust(top=.95, bottom=.05, hspace=.3, wspace=.3)
plt.show()

 

베타 분포 (Beta Distribution for different alpha and beta values)

 

 

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

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

 

 

728x90
반응형
Posted by Rfriend
,

확률과 우도(likelihood)는 통계 및 확률 이론에서 사용되는 관련된 개념으로, 각각 독특한 의미와 적용을 갖고 있습니다.

 

 

1. 확률 (Probability)

 

- 정의: 확률은 특정 이벤트가 발생할 가능성을 측정하며, 0에서 1 사이의 숫자로 나타납니다. 0은 불가능함을 나타내고, 1은 확실함을 나타내며, 0과 1 사이의 값은 다른 정도의 가능성을 나타냅니다.

 

- 해석: 예를 들어, 어떤 이벤트의 확률이 0.8이면 해당 이벤트가 발생할 확률이 80%인 것을 의미합니다.

 

 

 

2. 우도 (likelihood):

 

- 정의: 반면에 우도는 관찰된 데이터가 특정 가설이나 모델을 지원하는 정도를 측정합니다. 우도는 확률이 아니며 0에서 1의 범위에 제한되지 않습니다.

 

- 해석: 우도 함수는 특정 모델이나 가설이 관찰된 데이터를 얼마나 잘 설명하는지를 나타냅니다. 우도가 높을수록 데이터가 해당 가설과 일치하는 정도가 높습니다. 

 

likelihood mathematical expression

 

 

 

3. 주요 차이점

 

(1) 중점 사항 (Focus)

 

- 확률 (Probability): 알려진 확률에 기반하여 미래 이벤트의 가능성을 예측하는 데 중점을 둡니다.

- 우도 (Likelihood): 관찰된 데이터와 특정 가설 또는 모델 간의 일치를 평가하는 데 중점을 둡니다.

 

 

(2) 척도 (Scale)

 

- 확률 (Probability): 0에서 1까지의 범위를 가지며 이벤트 발생의 가능성을 나타냅니다.

- 우도 (Likelihood) : 특정 가설을 지원하는 정도를 나타내며 특정 범위에 제한되지 않습니다.

 

 

(3) 해석 (Interpretation)

 

- 확률 (Probability): 미래 이벤트의 확실성 또는 가능성을 측정합니다.

- 우도 (Likelihood) : 주어진 가설이 관찰된 데이터를 얼마나 잘 설명하는지를 나타냅니다.

 

 

(4) 응용 (Applications)

 

- 확률 (Probability): 예측 모델링, 의사 결정, 확률 게임 등에서 널리 사용됩니다.

- 우도 (Likelihood) : 통계적 추론, 가설 검정, 모델 적합 등에서 널리 사용됩니다.

 

 

 

4. 예시

 

동전 던지기를 고려해봅시다. 앞면이 나올 확률은 0.5입니다(공정한 동전을 가정). 이제 10번의 동전 던지기에서 8번이 앞면이 나왔다면, 해당 가설에 따라 동전이 공정한지 아니면 앞면이 더 많이 나오는 편향된 동전인지를 설명하는데 우도를 계산할 수 있습니다.

 

요약하면, 확률은 미래 이벤트의 가능성을 예측하는 데 사용되며, 우도는 관찰된 데이터와 특정 가설 간의 적합성을 평가하는 데 사용됩니다. 확률은 0에서 1 사이의 값으로 표현되며, 우도는 가설에 대한 지지 정도를 측정하며 특정 범위에 제한되지 않습니다.

 

 

 

5. 최대우도 추정법 (MLE, Maximum Likelihood Estimation)

 관측 결과를 바탕으로 관측되었을 가능성이 가장 큰 모수값으로 모수를 추정하는 방법을 최대가능도법(Maximum Likelihood Method) 이라고 합니다. 이를 수리적으로 정리해보면 다음과 같습니다. 

 

Maximum Likelihood Estimation, MLE, 최대우도추정법

 


실제로는 우도 함수를 계산을 간소화하기 위해 로그-우도( Log-Likelihood )로 변환되는 경우가 많습니다. 로그는 단조 증가 함수이므로 로그-우도는 동일한 최대값을 갖습니다.

 

우도는 통계적 추론에서 중요한 역할을 합니다. 다양한 모델이나 가설을 비교할 때 연구자들은 종종 그들의 우도를 비교하며, 높은 우도를 갖는 모델이 관측된 데이터에 대한 더 나은 설명으로 간주됩니다. 

 

 

 

6. 정규분포로 부터의 확률표본으로 평균과 분산의 최대가능도추정량 구하기

 

수식 입력기로 하나씩 수식 입력하는게 힘들어서 손으로 작성했습니다. ^^; 

이게 대학원 중간고사 볼 때 주관식으로 다 외워서 썼던 것인데, 다시 보니 새롭네요. ㅋㅋ

 

 

 

 

7. Python을 이용한 정규분포로 부터의 확률표본으로 평균과 분산의 최대가능도추정량 구하기

 

Python의 scipy 모듈을 사용해서 최대가능도추정법(MLE) 으로 모수 평균과 분산의 추정량을 구해보겠습니다. 

 

모집단의 평균이 4.5, 분산이 2.0인 정규분포로부터 표본을 1,000개 추출하여, 이 확률표본으로 부터 모집단의 분포의 모수 평균(mu), 분산(sigma^2)을 최대가능도추정(MLE) 기법을 사용해 추정해보겠습니다. 최대가능도추정량이 평균 4.48, 분산 2.0 으로서 모집단의 모수와 매우 근접하게 잘 추정되었음을 알 수 있습니다. 

 

import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize

# Generate synthetic data from a normal distribution
np.random.seed(1004)
true_mean = 4.5
true_std = 2.0
sample_size = 1000
data = np.random.normal(true_mean, true_std, sample_size)

# Define the negative log-likelihood function for a normal distribution
def neg_log_likelihood(params, data):
    mean, std = params
    log_likelihood = -np.sum(norm.logpdf(data, mean, std))
    return log_likelihood

# Initial guess for mean and std
initial_params = [1, 1]

# Use MLE to estimate mean and std
result = minimize(neg_log_likelihood, initial_params, args=(data,), method='L-BFGS-B')
estimated_mean, estimated_std = result.x

# Display the results
print(f"True Mean: {true_mean}, True Standard Deviation: {true_std}")
print(f"Estimated Mean: {estimated_mean}, Estimated Standard Deviation: {estimated_std}")

# True Mean: 4.5, True Standard Deviation: 2.0
# Estimated Mean: 4.485386426327861, Estimated Standard Deviation: 2.0015927811977487

 

 

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

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

 

728x90
반응형
Posted by Rfriend
,

상관관계와 인과관계는 통계 및 연구에서 두 변수 간의 다른 유형의 관계를 설명하는 두 가지 중요한 개념입니다. 


1. 상관관계 (Correlation)

 

- 정의: 상관관계는 두 변수가 함께 변하는 정도를 설명하는 통계적 측정입니다. 다시 말하면 한 변수의 변화가 다른 변수의 변화와 어떤 정도로 연관되어 있는지를 양적으로 표현합니다. 


- 측정지표: 상관계수(correlation coefficients)는 두 변수 간의 연관 정도를 측정하는 지표입니다. 상관을 측정하는 여러 방법 중에서 가장 일반적인 것은 피어슨 상관계수입니다. 피어슨 상관계수는 -1에서 1까지의 숫자로, 변수 간의 관계의 강도와 방향을 나타냅니다.

 

- 상관계수 해석:
  - 양의 상관관계 (+1에 가까움): 한 변수가 증가하면 다른 변수도 증가하는 경향이 있습니다. 
  - 음의 상관관계 (-1에 가까움): 한 변수가 증가하면 다른 변수는 감소하는 경향이 있습니다.
  - 제로 상관관계 (0에 가까움): 변수 간에 선형 관계가 없습니다.

 

- 예시: 아이스크림 판매와 익사 사고 수 사이에 양의 상관관계가 있다고 가정해 봅시다. 둘 다 여름 동안 증가하는 경향이 있습니다. 그러나 이는 아이스크림을 사 먹는 것이 익사를 유발한다거나 그 반대로 되는 것은 아니라는 의미입니다. 단지 상관관계일 뿐입니다.

 



2. 인과관계 (Causation)

 

- 정의: 인과관계는 두 변수 간의 원인과 결과 관계를 나타냅니다. 인과관계에서는 하나의 변수의 변화가 다른 변수에 직접적으로 영향을 주어 원인과 결과의 연결이 있습니다. 인과관계는 한 변수의 변화가 다른 변수에 변화를 일으키는 상황을 의미하며, 인과관계를 확립하려면 단순히 상관관계를 관찰하는 것 이상의 노력이 필요합니다. 이는 논리적인 연결을 입증하고 다른 가능한 설명을 배제하는 것을 포함합니다. 


- 인과관계 판단 기준: 원인을 추론하기 위해 종종 고려되는 여러 기준이 있습니다. 이에는 다음이 포함됩니다:

 

  (1) 시간적 우선순위 (Temporal precedence)

       : 원인은 효과보다 앞서야 함
  (2) 원인과 효과의 공변성 (Covariation of cause and effect)

       : 원인의 변화가 효과의 변화와 관련이 있어야 함
  (3) 대안적 설명의 제거 (Elimination of alternative explanations)

       : 다른 요인들이 원인으로 제외됨

 

- 예시: 흡연과 폐암은 인과 관계가 있습니다. 여러 연구에서 흡연이 폐암 발생 위험을 증가시킨다는 것이 확인되었습니다. 이 관계는 단순한 통계적 연관이 아니라 흡연이 폐암 발생 가능성에 영향을 미치는 알려진 메커니즘이 있습니다.

 

인과관계의 종류

 

3. 주요 차이점

 

(1) 관계의 성격 (Nature of Relationship)
- 상관관계: 두 변수 간의 통계적 관련성이나 관계를 나타냅니다.
- 인과관계: 변화한 변수가 직접적으로 다른 변수에 영향을 미치는 인과 관계를 나타냅니다.


(2) 방향 (Direction)
- 상관관계: 양의 경우 (두 변수가 함께 증가하거나 감소) 또는 음의 경우 (한 변수가 증가하면 다른 변수가 감소)가 될 수 있습니다.
- 인과관계: 원인 변수의 변화가 결과 변수에 변화를 일으키는 특정한 방향성이 있습니다.


(3) 시사점 (Implications)
- 상관관계: 인과관계를 의미하지 않습니다. 상관관계만으로는 인과관계의 증거를 제공하지 않습니다.
- 인과관계: 직접적이고 종종 예측적인 관계를 의미합니다. 인과관계를 확립하기 위해서는 일반적으로 보다 엄격한 실험 설계와 증거가 필요합니다.


(4) 예시
- 상관관계: 아이스크림 판매와 익사 사고 수는 서로 인과관계 없이 상관관계가 있을 수 있습니다.
- 인과관계: 흡연이 폐암을 유발하는 것은 인과관계의 한 예입니다.

통계적 관계를 해석할 때 주의가 필요하며, 상관관계만을 근거로 인과관계를 가정하는 것은 피해야 합니다. 보다 심층적인 연구, 실험, 대안적 설명의 고려가 종종 필요합니다.



4. Python 을 이용한 상관관계(Correlation) 분석 예시

 

import pandas as pd

# Example data as Pandas DataFrame
df = pd.DataFrame({
    'data1': [1, 2, 3, 4, 5], 
    'data2': [2, 2.4, 3, 3.2, 3.5], 
    'data3': [4, 3, 3, 2, 1.5]})

# Calculate the correlation matrix
correlation_matrix = df.corr()

# Print the correlation matrix
print("Correlation Matrix:")
print(correlation_matrix)

# Correlation Matrix:
#           data1     data2     data3
# data1  1.000000  0.985104 -0.973329
# data2  0.985104  1.000000 -0.933598
# data3 -0.973329 -0.933598  1.000000

 

 

5. Python 을 이용한 인과관계(Causation) 분석 예시

 

Python의 statsmodels 모듈을 사용해서 선형회귀모형(linear regression)을 적합해보겠습니다. 

 

선형회귀모형을 독립변수와 종속변수 간에 선형적인 관계를 가정하며, 인과관계에세는 원인이 결과보다 시간면에서 앞서서 일어나야 하므로 아래 예시에서는 설명변수 X가 목표변수 Y보다 먼저 발생한 데이터라고 가정하겠습니다.  만약 관측되지 않거나 통제되지 않은 교란 변수 (Confounding Variable) 가 있다면, 회귀 모델에서 추정된 계수는 편향될 수 있습니다. 모든 관련 변수를 통제하는 것은 인과 추론에 있어서 필수적이며, 아래 예시에 교란 변수 Z를 추가하여 선형회귀모형을 적합하였습니다. 

 

## install statsmodels in terminal first if you don't have it
! pip install numpy pandas statsmodels



## linear regression
import numpy as np
import pandas as pd
import statsmodels.api as sm

# Generate synthetic data for demonstration
np.random.seed(42)
X = np.random.rand(100, 1) * 10  # Independent variable
epsilon = np.random.randn(100, 1) * 2  # Error term
Y = 2 * X + 1 + epsilon  # Dependent variable with a causal relationship

# Create a DataFrame
df = pd.DataFrame({'X': X.flatten(), 'Y': Y.flatten()})

# Add a confounding variable
Z = np.random.randn(100, 1)  # Confounding variable
df['Z'] = Z.flatten()

# Fit a linear regression model
X_with_intercept = sm.add_constant(df[['X', 'Z']])
model = sm.OLS(df['Y'], X_with_intercept).fit()

# Display regression results
print(model.summary())

#                             OLS Regression Results                            
# ==============================================================================
# Dep. Variable:                      Y   R-squared:                       0.908
# Model:                            OLS   Adj. R-squared:                  0.906
# Method:                 Least Squares   F-statistic:                     479.8
# Date:                Sun, 10 Dec 2023   Prob (F-statistic):           5.01e-51
# Time:                        06:30:59   Log-Likelihood:                -200.44
# No. Observations:                 100   AIC:                             406.9
# Df Residuals:                      97   BIC:                             414.7
# Df Model:                           2                                         
# Covariance Type:            nonrobust                                         
# ==============================================================================
#                  coef    std err          t      P>|t|      [0.025      0.975]
# ------------------------------------------------------------------------------
# const          1.4314      0.342      4.182      0.000       0.752       2.111
# X              1.9081      0.062     30.976      0.000       1.786       2.030
# Z             -0.0368      0.167     -0.220      0.826      -0.368       0.295
# ==============================================================================
# Omnibus:                        0.975   Durbin-Watson:                   2.295
# Prob(Omnibus):                  0.614   Jarque-Bera (JB):                0.878
# Skew:                           0.227   Prob(JB):                        0.645
# Kurtosis:                       2.928   Cond. No.                         10.7
# ==============================================================================

# Notes:
# [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

 

 

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

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

 

728x90
반응형
Posted by Rfriend
,

베이지안 통계는 확률이나 정보의 새로운 증거에 따라 확률을 갱신하는 통계학의 한 분야입니다. 이는 18세기 수학자 토마스 베이즈(Thomas Bayes)의 이름에서 따왔으며 통계적 추론에 대한 확률적인 프레임워크를 제공합니다.

 

베이지안 통계의 주요 개념을 하나씩 소개하고, 예제를 하나 풀어보겠습니다. 

 

 

1. 베이즈 정리 (Bayes' Theorem) 및 베이지안 통계 주요 개념

 

베이지안 통계의 핵심은 베이즈 정리로, 사건의 확률을 사건과 관련된 조건에 대한 이전 지식을 기반으로 설명합니다. 수식은 다음과 같습니다. 여기서 는 사후 확률(조건 B가 주어졌을 때 A의 확률), 는 우도(조건 A가 주어졌을 때 B의 확률), 는 사전 확률, 는 B의 주변 확률입니다. 

 

 

베이즈 정리 (Bayes' Theorem)

 

 

 

- 사전 확률 (Prior Probability): 사전 확률은 새로운 증거를 고려하기 전의 사건에 대한 초기 믿음 또는 확률을 나타냅니다. 이는 이전 지식이나 경험에 기초합니다.


- 우도(Likelihood): 우도는 특정 가설이나 모델로 관측된 데이터를 설명할 수 있는 정도를 나타냅니다. 특정 가설이나 모델이 주어졌을 때 관측된 데이터의 확률을 나타냅니다.


- 사후 확률(Posterior Probability): 사후 확률은 새로운 증거를 고려한 후의 가설이나 사건의 갱신된 확률입니다. 베이즈 정리를 사용하여 계산됩니다.


- 사후 추론(Posterior Inference): 베이지안 통계에서는 매개변수의 점 추정 대신에 사후 분포에 중점을 둡니다. 이 분포는 가능한 값의 범위와 해당 확률을 제공합니다. 


- 베이지안 갱신(Bayesian Updating): 새로운 데이터가 생기면 베이지안 통계는 믿음과 확률을 계속해서 갱신할 수 있습니다. 한 번의 분석에서 얻은 사후 분포는 다음 분석의 사전 분포가 될 수 있습니다.

 

 

베이지안 통계는 데이터가 제한적인 경우나 사전 지식을 통합해야 하는 경우에 특히 유용합니다. 기계 학습, 데이터 과학, 경제학 및 의학 연구를 포함한 여러 분야에서 응용되며 불확실성을 일관되고 확률적으로 처리하는 유연성과 능력으로 알려져 있습니다. 

 

 

 

2. 빈도론자(Frequentist) vs. 베이지안(Bayesian) 통계 추론*

 

 통계추론은 크게 빈도론자(Frequentist), 베이지안(Bayesian)에 의한 추론으로 구분합니다. 예를 들어서 비교해보겠습니다. ( * 예시 출처: '통계학의 개념과 제문제', 이긍희, 김훈, 김재희, 박진호, 이재용 공저, KNOU출판부) 

 

 "예를 들어, 어느 학교에 학생이 입학하였습니다. 학교는 학업 능력이 낮은 학생은 탈락시키고 학업 능력이 높은 학생은 조기에 상급학교에 진학시킨다고 합니다. 두 학생의 학업능력을 평가하기 위해 여러 번의 시험과 과제물을 평가하였습니다. 

 빈도론자(Frequentist)의 추론은 학생의 학업능력은 고정되어 있다고 가정하고, 일정기간 동안 여러 번의 시험과 과제물 채점 결과로 그 학생의 능력을 평가합니다. 시험 회수를 늘릴수록 두 학생의 평가결과는 공정하게 비교된다고 봅니다. 

 반면에 베이지안(Bayesian)은 학생의 학업능력은 고정되어 있지 않다고 가정하고 시험과 과제물 채점 시마다 그 학생의 능력평가를 변경합니다. 먼저 학생의 능력은 모두 같다고 가정합니다(또는 분석자가 사전정보를 이용한다면 한 학생이 다른 학생보다 우수하다고 가정합니다). 첫 시험을 보았는데 어떤 학생이 시험을 잘 본다면 그 학생의 능력은 다른 학생보다 우수한 것으로 생각합니다. 다음 시험을 본 후 그 학생의 능력에 대한 판단을 수정합니다." 

 

 

베이지안의 추론 방식이 일면 합리적이고, 또 보통 사람이 추론하는 방식과 유사함을 알 수 있습니다. 특히, 관측할 수 있는 데이터의 개수가 제한적이고 적을 때 사전적으로 알고 있는 지식이나 주관적인 믿음을 사전 확률로 이용하고, 새로 관측된 데이터로 사후 확률을 업데이트 하는 베이지안 분석법이 매우 강력합니다. 그리고 사후 확률도 단 하나의 점 추정이 아니라 확률적으로 생각하는 방식도 유효하구요. 

 

빈도론자와 베인지안 간에 서로 논박하는 주요 내용은 아래를 참고하세요. 

 

 빈도론자(Frequentist)는 고정된 모수를 무한히 반복되는 표본에 대한 통계량의 표본분포를 바탕으로 추정하거나 검정합니다. 

 반면 베이지안(Bayesian)은 표본확률에 사전확률을 더하여 추정합니다(위의 Bayes' Theorem). 모수는 임의적이어서 확률분포를 가지고 있으며, 모든 추정과 검정은 주어진 데이터와 모수의 사전확률을 바탕으로 한 사후확률에 기반해서 진행됩니다. 

 빈도론자(Frequentist)는 베이지안의 결과가 지나치게 모수의 사전분포에 의존해서 결과가 일정하지 않고 계산시간이 많이 든다고 비판합니다. 

 반면 베이지안(Bayesian)은 빈도론자가 주어진 정보를 활용하지 않아 올바른 추정에 어려움이 있다고 비판합니다. 

 

 

(* 출처: '통계학의 개념과 제문제', 이긍희, 김훈, 김재희, 박진호, 이재용 공저, KNOU출판부) 

 

 

 

3. 베이지안 통계 추론에 의한 사후 확률 계산 예시

 

의료 암 진단 사례로 베이지안 갱신의 예를 들어보겠습니다. 

 

(문제) 의료 데이터에 의하면 사전에 알려진 어떤 암에 걸릴 확률이 0.1% 라고 합니다.  

암에 걸렸는지를 검사하는 진단에서 실제 암이 걸린 사람은 90%의 확률로 양성이 나오고, 실제로 암에 안 걸린 건강한 사람은 5%의 확률로 양성(오진단)이 나온다고 합니다.  

이럴 때 이 암 진단에서 양성이라고 판정이 나왔을 때, 실제로 암에 걸렸을 확률을 계산하세요. 

 

이해하기 쉽도록 도식화해서 순서대로 풀어보겠습니다. 

 

(1) 사전 확률 (Prior probability) 설정

 

아직은 개별 환자에 대한 정보가 없는 상태에서, 기존의 의료 데이터로 부터 얻은 사전 확률은 암에 걸릴 확률이 0.1%, 암에 걸리지 않을 확률이 99.9%라는 것을 이미 알고 있습니다. 

 

사전확률 (prior probability)

 

 

 

(2) 검사 정밀도에 따른 조건부 확률 (Likelihood)

 

위 문제에 보면, 암에 걸렸다는 조건이 주어졌을 때 검사가 양성일 조건부 확률(P(D|H))이 90%, 암에 걸리지 않고 건강하다는 조건이 주어졌을 때 검사가 양성일 조건부 확률이 5% 라고 합니다. 이를 표로 표현하면 아래와 같습니다. 

 

검사 정밀도에 따른 조건부 확률

 

 

(3) 4개 유형별 각각의 확률 계산

 

위에 주어진 확률들로 부터 아래와 같이 4개의 유형이 존재하게 되며, 각 유형별로 확률을 계산해볼 수 있습니다. 

 

- 암이면서 양성일 확률     = 0.001 * 0.90 = 0.09%

- 암이면서 음성일 확률     = 0.001 * 0.10 = 0.01%

- 건강하면서 양성일 확률 = 0.999 * 0.05 = 4.995%

- 건강하면서 음성일 확률 = 0.999 * 0.95 = 94.91%

 

4개 유형별 각각의 확률

 

 

 

(4) 검사가 양성이므로 음성 영역은 제거

 

우리가 관심있어하는 사람의 진단 결과가 양성이라고 했으므로, 위의 (3)번에서 진단 결과가 음성인 영역은 이 환자에게는 해당사항이 없으므로 제거합니다. 

 

 

 

 

(5) 정규화(normalization)한 후의 사후 확률 (Posterior Probability) 계산

 

진단 결과가 양성인 영역만 남겨둔 상태에서, 암이면서 양성인 확률 (=0.001 * 0.90 = 0.09%) 과 건강하면서 양성인 확률을 (0.999 * 0.05 = 4.995%) 모두 더해서, 이 값으로 암이면서 양성인 확률과 건강하면서 양성인 확률을 나누어 주면 정규화(normalization)이 된 사후 확률을 구할 수 있습니다. 

 

정규화(normalization) 후의 사후 확률(posterior probability)

 

 

위에서 손으로 풀어본 결과 암 진단 결과가 양성으로 나온 사람이 있다고 했을 때, 실제 암에 걸렸을 사후 확률은 0.09 / (0.09 + 4.995) = 1.76% 이네요. (실제 건강한데 진단 결과 양성이 나왔을 사후 확률은 4.995 / (0.09+4.994) = 98.23% (=1-1.76%) 입니다)

 

90%의 정확도로 암을 진단하는 검사의 결과에서 양성이 나온다면 일반인들은 '아, 이제 나 죽나보구나. 흑.. ㅠ.ㅠ' 하고 앞이 깜깜할 거 같습니다. 그런데 베이지안(Bayesian) 통계에 대한 지식이 있는 사람이라면, '흠... 사전확률을 보니 대부분의 사람은 암이 아니고 극소수만 암이기 때문에, 진단 결과가 양성이라는 새로운 정보를 업데이트 해도... 실제로 암일 사후 확률이 그리 높지는 않겠네. 추가 검사를 기다려보도록 하지' 라고 생각할 수 있겠습니다. 베이지안 통계를 알면 쫄지 않아도 됩니다! ㅋㅋ

 

위의 사례를 보면 통계 지식이 부족한 일반인들이 생각하는 결과와 베이지안 통계로 추정한 확률이 차이가 많이 나고, 베이지안 통계가 일반인들의 확률에 대한 상식에 반하는거 아닌가하는 생각마저 들지도 모르겠습니다. ^^; 

 

 

 

4. Python을 이용한 베이지안 갱신 (Bayesian Updating) 계산

 

베이즈 정리를 이용하여 사전확률, 조건부 확률이 주어졌을 때 사후확률을 구하는 사용자 정의 함수를 정의해보겠습니다. 

 

## Bayesian Update using Bayes's Theorem
import numpy as np

def bayesian_update(prior_prob, likelihood, verbose=True):
    # Bayes' theorem: P(H|D) = P(D|H) * P(H) / P(D)
    # Posterior Probability = Likelihood * Prior Probability / Evidence

    # Calculate the unnormalized posterior probability
    posterior_unnormalized = likelihood * prior_prob
    
    # Calculate the evidence (marginal likelihood) using the law of total probability
    evidence = sum(posterior_unnormalized)

    # Calculate the normalized posterior probability
    posterior_prob = posterior_unnormalized / evidence

    if verbose:
        print('Likelihood:', likelihood)
        print('Prior Probability:', prior_prob)
        print('Posterior Unnormalized:', posterior_unnormalized)
        print('Evidence(Marginal Likelihood):', evidence)
        print('Posterior Probability:', posterior_prob)

    return posterior_prob

 

 

어느 특정 암에 걸릴 확률이 0.1% 라고 하는 사전확률이 주어지고,

암에 걸렸을 때 암 검사에서 양성이 나올 조건부 확률이 90.0%,

암에 걸리지 않았을 때 양성이 나올 조건부 확률이 5% 라고 했을 때, 

만약 어떤 사람이 암 검사에서 양성이 나왔다면 실제로 암에 걸렸을 사후 확률을 계산해 보겠습니다. 

 

# Example values
prior_probability = np.array([0.001, 0.999])  # Prior probability for H0 and H1
likelihood_given_h0 = np.array([0.90, 0.05])  # Likelihood of the data given H0
#likelihood_given_h1 = np.array([0.10, 0.95])  # Likelihood of the data given H1

# Update the posterior probability
posterior_probability_h0 = bayesian_update(prior_probability, likelihood_given_h0)

# Print the results
print("-----" * 10)
print("Posterior Probability for H0:", posterior_probability_h0)


# Likelihood: [0.9  0.05]
# Prior Probability: [0.001 0.999]
# Posterior Unnormalized: [0.0009  0.04995]
# Evidence(Marginal Likelihood): 0.05085
# Posterior Probability: [0.01769912 0.98230088]
# --------------------------------------------------
# Posterior Probability for H0: [0.01769912 0.98230088]

 

 

암 검사에서 양성이 나왔을 때 실제로 암에 걸렸을 확률은 1.76% 네요.  

 

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

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

 

728x90
반응형
Posted by Rfriend
,

중심극한정리(Central Limit Theorem, 이하 CLT)는 통계학의 기본 개념으로, 특히 표본 크기가 충분히 큰 경우 모집단의 표본 평균 분포를 설명합니다.

 

1. 중심극한정리(Central Limit Theorem, CLT) 개념

 

(1) 무작위 추출 (Random Sampling): 중심극한정리는 모양이 어떤 분포든지 상관없이 주어진 모집단에서 고정된 크기의 무작위 표본을 추출하고 각 표본의 평균을 계산한다고 가정합니다.

 

(2) 표본 평균 분포 (Distribution of Sample Mean): CLT는 원래 모집단 분포의 모양과 상관없이 표본 평균의 분포가 샘플 크기가 증가함에 따라 정규 분포를 근사화한다고 말합니다. 

 

(3) 크기가 큰 표본 (Large Size of Sample): 정규 분포 근사는 표본 크기가 충분히 큰 경우에 특히 잘 적용되며 일반적으로 n ≥ 30입니다. 그러나 경우에 따라서는 적은 표본 크기에서도 CLT 근사가 가능할 수 있습니다. 특히 기본 모집단 분포가 극도로 편향되지 않은 경우입니다. 

 

 

 

이는 표본 평균의 분포(X bar)가 평균이 모집단 평균(mu)과 같고 표준 편차(sigma/sqrt(n))가 모집단 표준 편차를 표본 크기의 제곱근으로 나눈 값인 정규 분포에 근사적으로 따른다는 것을 의미합니다.

 

중심극한의 정리의 중요성은 그 널리 활용성에 있습니다. 중심극한의 정리를 통계학자와 연구자들은 원래 모집단 분포의 모양을 모르더라도 표본 평균의 분포에 대해 특정 가정을 할 수 있게 됩니다. 이는 통계적 추론(statistical inference), 가설 검정(hypothesis testing) 및 신뢰 구간을 추정(estimate of confidence interval)하는 데 중요하며, 많은 통계적 방법이 정규성 가정(hypothesis of normal distribution)에 의존하기 때문에 실무에서 이러한 통계적 기법의 기초 역할을 합니다.

 

 

 

2. Python을 이용한 중심극한정리 시뮬레이션 

 

다음은 Python을 사용해서 균등분포(uniform distribution)를 따르는 모집단(population)으로 부터 각 표본 크기가 30개인 표본을 1,000번 추출하는 시뮬레이션을 하여 각 표본의 평균을 계산하고, 표본 평균의 분포(distribution of sample mean)를 히스토그램으로 시각화해 본 것입니다. 

 

import numpy as np
import matplotlib.pyplot as plt

# Parameters
population_size = 100000  # Size of the population
sample_size = 30        # Size of each sample
num_samples = 1000      # Number of samples to generate

# Generate a non-normally distributed population (e.g., uniform distribution)
population = np.random.uniform(0, 1, population_size)

# Initialize an array to store the sample means
sample_means = []

# Generate samples and calculate means
for _ in range(num_samples):
    sample = np.random.choice(population, size=sample_size, replace=False)
    sample_mean = np.mean(sample)
    sample_means.append(sample_mean)

# Plot the population and the distribution of sample means
plt.figure(figsize=(12, 6))

# Plot the population distribution
plt.subplot(1, 2, 1)
plt.hist(population, bins=30, color='blue', alpha=0.7)
plt.title('Distribution of Uniform Distribution (Population)')

# Plot the distribution of sample means
plt.subplot(1, 2, 2)
plt.hist(sample_means, bins=30, color='green', alpha=0.7)
plt.title('Distribution of Sample Means')

plt.show()

 

central limit themrem simulation

 

왼쪽 히스토그램이 균등분포를 따르는 모집단의 것이고, 오른쪽 히스토그램은 균등분포를 따르는 모집단으로 부터 무작위 추출한 표본의 평균(sample mean)의 분포를 나타내는 것입니다. 오른쪽의 표본 평균의 분포는 평균 mu=(1+0)/2 = 0.5 를 중심으로 좌우 대칭의 종모양의 정규분포를 따름을 알 수 있습니다. 

 

위의 Python codes 에서 모집단(population) 을 이산형, 연속형 확률분포 별로 바꿔서 시뮬레이션을 해보시면 중심극한정리를 눈으로 확인해볼 수 있습니다. 

(가령, 베타분포의 경우 population = np.random.beta(0.8, 0.8, population_size) 처럼 수정하시면 됩니다)

 

 

 

3. 중심극한정리를 이용한 문제풀이 예시

 

중심극한정리를 이용하여 이항분포의 확률값을 정규분포로 근사해서 구할 수 있습니다. 이항분포는 n 이 비교적 작을 때 (통상 n <= 25) 정확한 확률값을 얻는데 유용하지만 n 이 클 때는 계산이 복잡해집니다. 이 경우 정규분포를 이용하여 이항분포의 확률의 근사값을 구할 수 있습니다. 즉, n이 클수록 또한 p값이 0.5에 가까울수록 이항분포는 정규분포화 유사해집니다. 이산형 확률변수 Y가 이항분포 B(n, p)를 따를 때, 이를 표준화(Z)한 후 표준정규분포 Z~N(0, 1)에 근사시킵니다. 

 

연속형 변수를 사용하여 이산형 변수를 근사할 때에는 연속성 수정계수(continuity correction factor)를 사용하여 오차를 줄입니다. 부등호에 따라 0.5를 가감함으로써 이항분포의 구간을 정규분포의 구간으로 수정합니다. 

 

[예시 문제] 한국인 성인 50%가 최소한 한 개의 신용카드를 가지고 있다고 가정하자.만일 30명의 성인을 표본추출할 때 19명에서 20명 사이의 성인이 최소한 한 개의 신용카드를 소지하고 있을 확률을
(1) 이항분포를 사용하여 구하고, 
(2) 정규근사를 사용하여 구하라. 

 

 

(1) 이항분포를 사용한 풀이

 

이산형 확률변수 Y가 최소한 한 개의 신용카드를 가지고 있을 성인의 수라고 하면, Y~B(30, 0.50)이 됩니다. 이에 대한 확률은 다음과 같이 계산합니다. 

 

이항분포를 사용한 풀이

 

 

 

(2) 정규근사를 사용한 풀이

 

Y~B(30, 0.50)일 때 E(Y)=n*p=30*0.50=15, Var(Y)=n*p*(1-p)=30(0.50)(1-0.50)=7.5 이고, 표준화의 변수를 Z라고 할 때, 이를 정규분포로 근사시킨 후 연속성 수정(+-0.5)을 하여 계산하면 다음과 같습니다. 

 

정규근사를 이용한 풀이

 

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

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

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 통계학에서 사용하는 p-값(p-value)이란 무엇이고, 어떻게 해석하는지 소개하겠습니다. 

 

1. 통계학에서 p-value 는 무엇이며, 무엇에 사용하는가? 

2. 독립된 두 표본 간 평균 차이 t-검정에서 p-value 해석 예시

3. Python으로 독립 두 표본 t-검정 (independent two-sample t-test)

4. 빈도론자(Frequentist)와 베이지안(Bayesian) 간 p-value 해석의 차이점은 무엇인가?

 

 

 

1. 통계학에서 p-value 는 무엇이며, 무엇에 사용하는가? 


통계학에서 p-값(확률 값)은 귀무가설(null hypothesis, H0)에 대한 증거를 평가하는 데 도움이 되는 지표입니다. 귀무가설은 차이나 효과가 없다는 것을 나타내는 명제로, 연구자들은 이 귀무가설을 기각할 충분한 증거가 있는지를 판단하기 위해 통계 검정을 사용합니다. 

p-값은 귀무가설이 참일 때 관측 결과나 그보다 더 극단적인 결과가 나올 확률을 나타냅니다. 다시 말해, 이는 귀무가설에 대한 증거의 강도를 측정하는 것입니다.

 


주로 다음과 같은 방식으로 사용합니다. 

(1) 만약 p-값이 작다면(일반적으로 미리 정의된 유의수준, 예를 들어 0.05보다 낮다면), 이는 관측 결과가 무작위로 발생할 가능성이 적다는 것을 나타내며, 이는 귀무가설(null hypothesis, H0)을 기각하고 대립가설(alternative hypothesis, H1)을 채택합니다.

(2) 만약 p-값이 크다면, 이는 관측 결과가 귀무가설과 상당히 일치한다는 것을 나타내며, 귀무가설을 기각할 충분한 증거가 없다는 것을 의미합니다.

 


p-값이 작다고 해서 특정 가설이 참임을 증명하는 것은 아닙니다. 이는 단순히 데이터가 귀무가설과 일치하지 않는다는 것을 나타냅니다. 또한, 유의수준(예: 0.05)의 선택은 어느 정도 임의적이며 연구자의 판단이나 해당 분야의 관례에 따라 달라집니다.

연구자들은 p-값을 신중하게 해석하고 귀무가설 검정에 관한 결정을 내릴 때 다른 관련 정보와 함께 고려해야 합니다. p-값은 효과의 크기나 실제적인 중요성에 대한 정보를 제공하지 않으며, 단지 귀무가설에 대한 증거의 강도(the evidence against a null hypothesis)를 나타냅니다

 



2. 독립된 두 표본 간 평균 차이 t-검정에서 p-value 해석 예시

 

예를 들어 혈압을 낮추는 새로운 약물에 대한 임상시험을 다루는 예를 살펴보겠습니다. 연구자들은 새로운 약물이 플라시보에 비해 혈압을 낮추는 데 효과적인지를 테스트하고자 합니다. 이 맥락에서 p-값의 사용과 해석은 다음과 같을 수 있습니다.

시나리오:

- 귀무가설 (H0): 새로운 약물은 혈압에 아무런 영향을 미치지 않으며, 관측된 차이는 무작위로 발생한 것이다.
- 대립가설 (H1): 새로운 약물은 플라시보에 비해 혈압을 낮추는 데 효과적이다. 즉, 플라시보 그룹보다 신약 그룹의 혈압이 더 낮다.


실험 설계:

- 연구자들은 참가자들을 새로운 약물을 받는 그룹과 플라시보를 받는 그룹으로 무작위로 배정합니다.
- 특정한 치료 기간 이후 각 참가자의 혈압 변화를 측정합니다.


데이터 분석:

- 각 10개의 샘플 데이터를 수집하고 분석한 후, 연구자들은 두 그룹 간의 혈압 변화의 평균을 비교하기 위해 통계적 검정(예: t-검정)을 수행합니다.
- 예를 들어, 검정 결과로 얻은 t-통계량이 3.372, p-값이 0.0017인 경우를 가정해 봅시다.


해석:

- 얻은 p-값(0.0034)이 선택된 유의수준(예: alpah=0.05)보다 작습니다.  
- 빈도주의적 해석: 이는 만약 귀무가설이 사실이라면(즉, 새로운 약물이 효과가 없다면), 실제로 관측된 데이터보다 더 극단적인 데이터를 관측할 확률이 0.17%밖에 되지 않는다는 것을 나타냅니다. 0.0017이 유의수준 0.05 보다 작으므로 연구자들은 귀무가설을 기각하기로 결정할 수 있습니다.
- 결정: 연구자들은 귀무가설을 기각하고 새로운 약물이 플라시보에 비해 혈압을 낮추는 데 효과적임을 시사하는 충분한 증거가 있다고 결론지을 수 있습니다. 

 

아래에는 t-분포표를 수작업으로 해석하는 방법입니다. 

t-분포표의 행은 자유도(degree of freedom) 으로서, 만약 샘플의 관측치 개수가 주어지면 (자유도 =  n-1) 로 계산해줍니다. t-분포표의 열은 유의수준 알파(alpha)이고, 표의 안에 있는 숫자는 가설 검정에 사용하는 t-통계량입니다. (정규분포표에서는 행과 열이 z-통계량이고, 표의 가운데 값이 p-value인데요, t-분포표는 이와 달라서 좀 헷갈리는 면이 있습니다.) 샘플로 부터 자유도와 t-통계량를 계산하면, 아래의 t-분포표를 보고 상단 열에서 p-값을 찾으면 됩니다.  

 

이번 예에서는 가령, 관측치가 각 그룹별로 10개라고 하면, 자유도 = n-1 = 9 가 됩니다. 

그리고 t-통계량을 계산했더니 3.372가 나왔다고 했을 때, 아래 표에서 t-통계량 3.372 는 자유도 9인 행에서 보면 3.250보다는 크고, 3.690보다는 작은 값입니다. 따라서 상단의 알파(alpha) 값은 0.005와 0.0025 사이의 값이라고 추정할 수 있습니다. 그런데 대립가설이 "H1: 새로운 약물은 플라시보에 비해 혈압을 낮추는 데 효과적이다. 즉, 플라시보 그룹보다 신약 그룹의 혈압이 더 낮다." 이므로 단측검정 (alternative = 'greater') 을 사용하므로, p-value=0.005/2=0.0025 보다 작고 ~ p-value=0.0025/2=0.00125 보다는 큰 값으로 볼 수 있습니다. 

 

유의수준 0.05보다는 훨씬 작은 값이므로, 귀무가설을 기각하고 대립가설을 채택합니다. 즉, 신약이 혈압을 낮추는데 효과가 있다고 판단할 수 있습니다. 

 

 

t-분포표 해석하기

 


p-값 자체만으로는 효과의 크기나 임상적 중요성에 대한 정보를 제공하지 않습니다. 연구자들은 통계 검정을 기반으로 의사결정을 내릴 때 효과 크기, 신뢰구간 및 실제적인 중요성도 고려해야 합니다. 또한 베이지안 프레임워크에서는 관측된 데이터를 기반으로 사전 신념을 업데이트하는 것이 포함될 수 있습니다.

 

 

3. Python으로 독립 두 표본 t-검정 (independent two-sample t-test)

 

scipy 모듈의 scipy.stats.stats.ttest_ind() 메소드를 사용해서 독립된 두 표본의 평균에 차이가 있는지를 t-검정해 보겠습니다.

 

- 귀무가설(H0): Group 1과 Group 2의 평균에는 차이가 없다. (mu_g1 = mu_g2)

- 대립가설(H1): Group 1가 Group 2의 평균보다 크다. (mu_g1 > mu_g2)

 

t 통계량이 3.37이고 p-값(p-value)이 0.0034 로서, 유의수준(alpha) 0.05보다 p-값이 더 작으므로 귀무가설을 기각(reject null hypothesis)하고, 대립가설을 채택(accept alternative hypothesis)합니다. 

 

## independent two-sample t-test
import numpy as np
import scipy.stats as stats

# Example data for two independent samples
group1 = [25, 30, 32, 28, 34, 26, 30, 29, 31, 27]
group2 = [21, 24, 28, 25, 30, 22, 26, 23, 27, 24]

print('Mean of Group 1:', np.mean(group1))
print('Mean of Group 2:', np.mean(group2))

# Perform independent two-sample t-test
t_statistic, p_value = stats.ttest_ind(
    group1, group2, 
    equal_var=True,
    alternative='greater', # alternative hypothesis: mu_1 is greater than mu_2
)

# Print the results
print(f'T-statistic: {t_statistic}')
print(f'P-value: {p_value}')

# Check the significance level (e.g., 0.05)
alpha = 0.05
if p_value < alpha:
    print('Reject the null hypothesis. \
There is a significant difference between the groups.')
else:
    print('Fail to reject the null hypothesis. \
There is no significant difference between the groups.')
    
# Mean of Group 1: 29.2
# Mean of Group 2: 25.0
# T-statistic: 3.3723126837047914
# P-value: 0.0016966232763530164
# Reject the null hypothesis. There is a significant difference between the groups.

 




4. 빈도론자(Frequentist)와 베이지안(Bayesian) 간 p-value 해석의 차이점은 무엇인가? 

 

빈도론자(Frequentist)와 베이지안(Bayesian) 간의 통계학적 접근 방식에서 p-값의 해석은 다릅니다. 각 접근 방식이 p-값을 어떻게 해석하는지 간략하게 살펴보겠습니다.

(1) 빈도주의적 해석 (Frequentist Interpretation):

- 빈도주의 통계학에서 p-값은 귀무가설이 참일 때 관측된 데이터나 그보다 더 극단적인 데이터가 나올 확률로 간주됩니다.
- 의사결정 과정은 일반적으로 특정한 유의수준(alpha), 보통 0.05,을 기준으로 귀무가설을 기각하거나 기각하지 않는 것으로 구성됩니다. 
- 만약 p-값이 선택된 유의수준보다 작거나 같으면 귀무가설이 기각됩니다. p-값이 유의수준보다 크면 귀무가설이 기각되지 않습니다.
- 빈도주의 통계학에서는 확률을 가설에 할당하지 않으며 대신 특정 가설 하에서 데이터의 확률에 중점을 둡니다.


(2) 베이지안 해석 (Bayesian Interpretation):

- 베이지안 통계학에서는 관측된 데이터를 기반으로 가설에 대한 사전 신념(prior beliefs)을 업데이트하는 것에 중점을 둡니다.
- 베이지안 분석은 가설에 대한 사전 신념과 관측된 데이터 모두를 고려하여 가설에 대한 사후 확률 분포(posterior probability distribution for the hypothesis)를 제공합니다.
- 베이지안은 의사결정을 위해 일반적으로 p-값 만을 고려하지 않습니다. 대신 전체 사후 분포를 고려하며 이전 정보를 통합합니다. 
- 베이지안 분석은 데이터가 주어졌을 때 가설이 참일 확률을 직접 계산합니다. 이는 빈도주의적 접근과 달리 가설이 주어진 데이터에 대한 확률만을 고려하는 것이 아닙니다.
- 가설 검정을 위해 때로는 베이지안 분석에서는 가설 간의 가능도 비율(ratio of the likehoods under different hypotheses)인 베이즈 팩터(Bayes factor)를 사용하기도 합니다.

요약하면 핵심적인 차이점은 확률 해석에 있습니다. 빈도주의자들은 확률을 사건의 장기적 빈도로 간주하며, 베이지안은 확률을 신념이나 불확실성의 척도로 해석합니다. 빈도주의자들은 p-값을 가설에 대한 의사결정에 사용하며, 베이지안안은 베이지안 추론을 사용하여 직접적으로 신념을 업데이트하고 표현합니다.

 

 

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

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 독립된 2개 표본의 평균 차이를 검정하는 t-test 에 대해서 알아보았습니다. 
이번 포스팅에서는 Python을 이용한 짝을 이룬 표본의 평균 차이를 검정하는 paired t-test 에 대해서 소개하겠습니다. 
(R을 이용한 짝을 이룬 표본에 대한 평균 차이 검정은 https://rfriend.tistory.com/128 를 참고하세요)

짝을 이룬 t-test(paired t-test)는 짝을 이룬 측정치(paired measurements)의 평균 차이가 있는지 없는지를 검정하는 방법입니다. 데이터 값이 짝을 이루어서 측정되었다는 것에 대해 예를 들어보면, 사람들의 그룹에 대해 신약 효과를 알아보기 위해 (즉, 신약 투약 전과 후의 평균의 차이가 있는지) 신약 투약 전과 후(before-and-after)를 측정한 데이터를 생각해 볼 수 있습니다.  

 

 

[ One sample t-test vs. Independent samples t-test vs. Paired samples t-test ]

paired samples t-test, image source: https://datatab.net/tutorial/paired-t-test

 


짝을 이룬 t-test(paired t-test)는 종속 표본 t-test (the dependent samples t-test), 짝을 이룬 차이 t-test (the paired-difference t-test), 매칭된 짝 t-test (the matched pairs t-test), 반복측정 표본 t-test (the repeated-sample t-teset) 등의 이름으로도 알려져있습니다. 동일한 객체에 대해서 전과 후로 나누어서 반복 측정을 합니다.  

 

아래의 도표에는 짝을 이룬 두 표본의 대응 비교 데이터셋에 대한 모습입니다. (Xi, Yi) 가 동일한 대상에 대해서 before-after 로 반복 측정되어서, 각 동일 객체의 전-후의 차이 (즉, Di = Xi - Yi) 에 대해서 검정을 진행하게 됩니다. 

 

paired t-test



[ 가정사항 (Assumptions) ]

(1) 측정 대상이 독립적(independent)이어야 합니다. 하나의 객체에 대한 측정은 어떤 다른 객체의 측정에 영향을 끼치지 않아야 합니다.  
(2) 각 짝을 이룬 측정치는 동일한 객체로 부터 얻어야 합니다. 예를 들면, 신약의 효과를 알아보기 위해 투약 전-후(before-after) 를 측정할 때 동일한 환자에 대해서 측정해야 합니다. 
(3) 짝을 이뤄 측정된 전-후의 차이 값은 정규분포를 따라야한다는 정규성(normality) 가정이 있습니다. 만약 정규성 가정을 충족시키지 못하면 비모수 검정(nonparametric test) 방법을 사용해야 합니다. 

 


[ (예제) 신약 치료 효과 여부 검정 ]

 

새로운 당뇨병 치료제를 개발한 제약사의 예를 계속 들자면, 치료에 지대한 영향을 주는 외부요인을 통제하기 위해 10명의 당뇨병 환자를 선별하여 1달 동안 '위약(placebo)'을 투여한 기간의 혈당 (Xi)과 동일 환자에게 '신약(new medicine)'을 투여한 1달 기간 동안의 혈당 수치(Yi)를 측정하여 짝을 이루어 혈당 차이를 유의수준 5%에서 비교하는 방법이 짝을 이룬 표본에 대한 검정이 되겠습니다. (palacebo 와 신약 투여 순서는 무작위로 선정. 아래 예는 그냥 예시로 아무 숫자나 입력해본 것임. 혈당 수치 이런거 전 잘 몰라요. ^^;)

 

* 귀무가설 (Null Hypothesis, H0): 신약 투입 효과가 없다 (Mu1 = Mu2, ie. Difference=0)
* 대립가설 (Alternative Hypothesis, H1): 신약 투입 효과가 있다 (Mu1 > Mu2, ie. Difference > 0, right-sided test)

 

 

[ Python scipy 모듈을 이용한 paired t-test 실행 ]

 

scipy 모듈의 scipy.stats.ttest_rel 메소드를 사용해서 쌍을 이룬 t-test 를 실행합니다. "Calculate the t-test on TWO RELATED samples of scores, a and b." 라는 설명처럼 TWO RELATED samples 에서 rel 을 타서 메소드 이름을 지었습니다. (저라면 ttest_paired 라고 메소드 이름 지었을 듯요...) 

 

alternative='two-sided' 가 디폴트 설정인데요, 이번 예제에서는 'H1: 신약이 효과가 있다. (즉, 신약 먹기 전보다 신약 먹은 후에 혈당이 떨어진다)' 는 가설을 검정하기 위한 것이므로 alternative='greater' 로 설정을 해주었습니다. 

 

## -- Paired t-test

import numpy as np
from scipy import stats

## sample data-set (repeated measurements of Before vs. After for the same objects)
bef = np.array([51.4, 52.0, 45.5, 54.5, 52.3, 50.9, 52.7, 50.3, 53.8, 53.1])
aft = np.array([50.1, 51.5, 45.9, 53.1, 51.8, 50.3, 52.0, 49.9, 52.5, 53.0])


## paired t-test using Python scipy module
# H0: New medicine is not effective (i.e., no difference b/w before and after)
# H1: New medicine is effective (i.e., there is difference b/w before and after)
stat, p_val = stats.ttest_rel(bef, aft, alternative='greater')

print('statistic:', stat, '   p-value:', p_val)
# statistic: 3.550688262985491    p-value: 0.003104595950799298

 

분석 결과 p-value 가 0.003 으로서 유의수준 0.05 하에서 귀무가설 (H0: 신약은 효과가 없다. 즉, before와 after의 차이가 없다) 을 기각(reject)하고, 대립가설(H1: 신약은 효과가 있다. 즉, before 보다 after의 혈당 수치가 낮아졌다)을 채택(accept) 합니다. 




[ Reference ]
* The Paired t-test
  : https://www.jmp.com/en_nl/statistics-knowledge-portal/t-test/paired-t-test.html
* scipy.stats.ttest_rel 
  (Calculate the t-test on TWO RELATED samples of scores, a and b.)
  : https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_rel.html

 

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

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

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 Python을 사용해서 두 집단 간 평균이 같은지 아니면 다른지를 검정하는 t-test 를 해보겠습니다. 

 

연속형 확률분포인 t-분포 (Student's t-distribution) 에 대해서는 https://rfriend.tistory.com/110 를 참고하세요. 

R을 사용한 독립된 두 집단간 모평균 차이에 대한 검정은 https://rfriend.tistory.com/127 를 참고하세요. 

 

모집단의 평균과 분산에 대해서는 알지 못하는 경우가 많으므로, 보통은 모집단에서 무작위로 표본을 추출(random sampling)해서 모집단의 평균과 분산을 추정합니다. 표본의 크기가 작은 집단 간 평균의 차이가 있는지를 검정할 때 t-분포에 기반한 t-통계량(t-statistics)을 사용하여 검정을 합니다. 

 

t-검정은 대상 표본 집단이 1개인지 2개인지에 따라서 아래와 같이 구분할 수 있습니다. 

  * One-sample t-test : 모집단의 평균이 귀무가설의 특정 평균 값과 같은지를 검정

  * Two-sample t-test: 두 모집단의 평균이 같다는 귀무가설을 검정

 

One-sample t-test와 Two-sample t-test에서 사용하는 통계량에 대해서는 아래에 정리해보았습니다. 

 

 

여기서부터는 독립된 두 표본 간의 평균 차이에 대한 t-검정 (independent two-sample t-test) 에 대해서만 자세하게 소개하도록 하겠습니다. 

 

(1) Two-sample t-test 의 가설 (Hypothesis)

 

 - 귀무가설 (Null Hypothesis, H0): Mu1 = M2 (두 모집단의 평균이 같다)

 - 대립가설 (Alternative Hypothesis, H1)

    -. 양측검정 대립가설 (two-sided test H1): Mu1 ≠ Mu2 (두 모집단의 평균이 같지 않다)

    -. 우측검정 대립가설 (right-tailed test H1): Mu1 > M2 (모집단1의 평균이 모집단2의 평균보다 크다)

    -. 좌측검정 대립가설 (left-tailed test H1): M1 < M2 (모집단1의 평균이 모집단2의 평균보다 작다) 

 

t-test 를 통해 나온 p-value 가 유의수준보다 작으면 귀모가설을 기각하고 대립가설을 채택(즉, 두 모집단의 평균이 차이가 있다)하게 됩니다. 

 

 

 

(2) Two-sample t-test 의 가정사항 (Assumptions)

 

Two-sample t-test 의 결과가 유효하기 위해서는 아래의 가정사항을 충족시켜야 합니다. 

 

 (a) 한 표본의 관측치는 다른 표본의 관측치와 독립이다. (independent) 

 (b) 데이터는 정규분포를 따른다. (normally distributed)

 (c) 두 집단의 표본은 동일한 분산을 가진다. (the same variance).

       (--> 이 가설을 만족하지 못하면 Welch's t-test 를 실행합니다.)

 (d) 두 집단의 표본은 무작위 표본추출법을 사용합니다. (random sampling)

 

정규성 검정(normality test)을 위해서 Kolmogorov-Smirnov test, Shapiro-Wilk test, Anderson-Darling test 등을 사용합니다. 등분산성 검정(Equal-Variance test) 을 위해서 Bartlett test, Fligner test, Levene test 등을 사용합니다. 

 

 

 

(3) Python을 이용한 Two-sample t-test 실행 

 

(3-1) 샘플 데이터 생성

 

먼저 numpy 모듈을 사용해서 정규분포로 부터 각 관측치 30개를 가지는 표본을 3개 무작위 추출해보겠습니다. 이중 표본집단 2개는 평균과 분산이 동일한 정규분포로 부터 무작위 추출하였으며, 나머지 1개 집단은 평균이 다른 정규분포로 부터 무작위 추출하였습니다. 

 

## generating sample dataset
import numpy as np

np.random.seed(1004) # for reproducibility
x1 = np.random.normal(loc=0, scale=1, size=30) # the same mean
x2 = np.random.normal(loc=0, scale=1, size=30) # the same mean
x3 = np.random.normal(loc=4, scale=1, size=30) # different mean


x1
# array([ 0.59440307,  0.40260871, -0.80516223,  0.1151257 , -0.75306522,
#        -0.7841178 ,  1.46157577,  1.57607553, -0.17131776, -0.91448182,
#         0.86013945,  0.35880192,  1.72965706, -0.49764822,  1.7618699 ,
#         0.16901308, -1.08523701, -0.01065175,  1.11579838, -1.26497153,
#        -1.02072516, -0.71342119,  0.57412224, -0.45455422, -1.15656742,
#         1.29721355, -1.3833716 ,  0.3205909 , -0.59086187, -1.43420648])

x2
# array([ 0.60998011,  0.51266756,  1.9965168 ,  1.42945668,  1.82880165,
#        -1.40997132,  0.49433367,  0.9482873 , -0.35274099, -0.15359935,
#        -1.18356064, -0.75440273, -0.85981073,  1.14256322, -2.21331694,
#         0.90651805,  2.23629   ,  1.00743665,  1.30584548,  0.46669171,
#        -0.49206651, -0.08727244, -0.34919043, -1.11363541, -1.71982966,
#        -0.14033817,  0.90928317, -0.60012686,  1.03906073, -0.03332287])

x3
# array([2.96575604, 4.15929405, 4.33053582, 4.02563551, 3.90786096,
#        3.08148823, 4.3099129 , 2.75788362, 3.66886973, 2.35913334,
#        3.72460166, 3.94510997, 5.50604364, 2.62243844, 2.74438348,
#        4.16120867, 3.57878295, 4.2341905 , 2.79844805, 5.48131392,
#        4.29105321, 4.4022031 , 3.58533963, 5.00502917, 5.45376705,
#        3.92961847, 4.52897801, 1.62104705, 3.24945253, 5.10641762])


## Box plot
import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(10, 8))
sns.boxplot(data=[x1, x2, x3])
plt.xlabel("Group", fontsize=16)
plt.ylabel("Value", fontsize=16)
plt.xticks([0, 1, 2], ["x1", "x2", "x3"], fontsize=14)
plt.show()

 

(3-2) t-test 가설 충족 여부 검정

 

t-검정의 가정사항으로서 정규성 검정(normality test)과 등분산성 검정 (equal variance test) 을 Python의 scipy 모듈을 사용해서 수행해보겠습니다. 

 

* Kolmogorov-Smirnov Test 정규성 검정 

  - (귀무가설, H0): 집단의 데이터가 정규 분포를 따른다. 

  - (대립가설, H1): 집단의 데이터가 정규 분포를 따르지 않는다. 

 

아래에 x1, x2, x3 의 세 집단에 대한 K-S 정규성 검정 결과 p-value 가 모두 유의수준 0.05 보다 크므로 귀무가설을 채택하여 세 집단의 데이터가 정규 분포를 따른다고 볼 수 있습니다. 

 

## (1) Normality test using Kolmogorov-Smirnov Test
import scipy.stats as stats

t_stat_x1, p_val_x1 = stats.kstest(x1, 'norm', args=(x1.mean(), x1.var()**0.5))
t_stat_x2, p_val_x2 = stats.kstest(x2, 'norm', args=(x2.mean(), x2.var()**0.5))
t_stat_x3, p_val_x3 = stats.kstest(x3, 'norm', args=(x3.mean(), x3.var()**0.5))

print('[x1]  t-statistics:', t_stat_x1, '  p-value:', p_val_x1)
print('[x2]  t-statistics:', t_stat_x2, '  p-value:', p_val_x2)
print('[x3]  t-statistics:', t_stat_x3, '  p-value:', p_val_x3)

# [x1]  t-statistics: 0.13719205314969185   p-value: 0.577558008887932
# [x2]  t-statistics: 0.11086245840821829   p-value: 0.8156064477001308
# [x3]  t-statistics: 0.09056001868899977   p-value: 0.9477307432911599

 

 

다음으로 집단 x1과 x2, 집단 x1과 x3에 대한 등분산 가정 검정 결과, p-value 가 모두 유의수준 0.05 보다 크므로 두 집단 간 분산이 같다고 할 수 있습니다. (귀무가설 H0: 두 집단 간 분산이 같다.)

 

## (2) Equal variance test using Bartlett's tes
var_test_stat_x1x2, var_test_p_val_x1x2 = stats.bartlett(x1, x2)
var_test_stat_x1x3, var_test_p_val_x1x3 = stats.bartlett(x1, x3)

print('[x1 vs. x2]', 'statistic:', var_test_stat_x1x2, '  p-value:', var_test_p_val_x1x2)
print('[x1 vs. x3]', 'statistic:', var_test_stat_x1x3, '  p-value:', var_test_p_val_x1x3)

# [x1 vs. x2] statistic: 0.4546474955289549   p-value: 0.5001361557169177
# [x1 vs. x3] statistic: 0.029962346601998174   p-value: 0.8625756934286083

 

처음에 샘플 데이터를 생성할 때 정규분포로 부터 분산을 동일하게 했었으므로 예상한 결과대로 잘 나왔네요. 

 

 

 

(3-3) 독립된 두 표본에 대한 t-test 평균 동질성 여부 검정

 

이제 독립된 두 표본에 대해 t-test 를 실행해서 두 표본의 평균이 같은지 다른지 검정을 해보겠습니다. 

 

 - (귀무가설 H0) Mu1 = Mu2 (두 집단의 평균이 같다)

 - (대립가설 H1) Mu1 ≠ Mu2 (두 집단의 평균이 다르다) 

 

분산은 서로 같으므로 equal_var = True 라고 매개변수를 설정해주었습니다. 

그리고 양측검정(two-sided test) 을 할 것이므로 alternative='two-sided' 를 설정(default)해주면 됩니다. (왜그런지 자꾸 에러가 나서 일단 코멘트 부호 # 로 막아놨어요. scipy 버전 문제인거 같은데요... 흠... 'two-sided'가 default 설정이므로 # 로 막아놔도 문제는 없습니다.)

 

## (3) Identification test using Independent 2 sample t-test

## x1 vs. x2
import scipy.stats as stats
t_stat, p_val = stats.ttest_ind(x1, x2, 
                                #alternative='two-sided', #‘less’, ‘greater’
                                equal_var=True)

print('t-statistic:', t_stat, '   p-value:', p_val)
#t-statistic: -0.737991822907993    p-value: 0.46349499774375136
#==> equal mean


## x1 vs. x3
import scipy.stats as stats
t_stat, p_val = stats.ttest_ind(x1, x3, 
                                #alternative='two-sided', #‘less’, ‘greater’
                                equal_var=True)

print('t-statistic:', t_stat, '   p-value:', p_val)
#t-statistic: -15.34800563666855    p-value: 4.370531118607397e-22
#==> different mean

 

(3-1)에서 샘플 데이터를 만들 때 x1, x2 는 동일한 평균과 분산의 정규분포에서 무작위 추출을 하였으며, x3만 평균이 다른 정규분포에서 무작위 추출을 하였습니다. 

 

위의 (3-3) t-test 결과를 보면 x1, x2 간 t-test 에서는 p-value 가 0.46으로서 유의수준 0.05 하에서 귀무가설(H0)을 채택하여 두 집단 x1, x2 의 평균은 같다고 판단할 수 있습니다. 

 

x1, x3 집단 간 t-test 결과를 보면 p-value 가 4.37e-22 로서 유의수준 0.05 하에서 귀무가설(H0)을 기각(reject)하고 대립가설(H1)을 채택(accept)하여 x1, x3 의 평균이 다르다고 판단할 수 있습니다. 

 

 

[ Reference ]

* Wikipedia Student's t-test: https://en.wikipedia.org/wiki/Student%27s_t-test

* Python scipy.stats.ttest_ind 메소드
: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.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
,