[확률] 9.1 확률분포의 추정

JKH·2024년 11월 27일

확률

목록 보기

19/25

9.1 확률분포의 추정

분석하고자 하는 데이터가 어떤 확률변수로부터 실현된 표본이다.

데이터 분석의 첫 번째 가정이다. 이 말은 우리가 정말 관심이 있는 것이 지금 손에 가지고 있는 데이터 즉, 하나의 실현체에 불과한 표본이 아니라 그 뒤에서 이 데이터를 만들어내는 확률변수의 분포라는 뜻이다. 확률론적인 관점에서 볼 때 데이터는 이 확률변수의 분포를 알아내기 위한 일련의 참고 자료일 뿐이다. 따라서 우리는 데이터 즉 표본으로부터 확률변수의 분포를 알아내야 한다.

확률분포의 결정

확률분포를 알아내는 일은 다음처럼 두 작업으로 나뉜다.

확률변수가 베르누이분포, 이항분포, 정규분포 등에서 어떤 확률분포를 따르는지 알아낸다.
데이터로부터 해당 확률분포의 모수의 값을 구한다.

첫 번째 작업 즉, 확률변수가 어떤 확률분포를 따르는가는 데이터가 생성되는 원리를 알거나 데이터의 특성을 알면 추측할 수 있다. 히스토그램을 그려서 확률분포의 모양을 살펴보고 힌트를 얻을 수도 있다.

데이터는 0 또는 1 뿐이다. $\rightarrow$ 베르누이분포
데이터는 카테고리 값이어야 한다. $\rightarrow$ 카테고리분포
데이터는 0과 1 사이의 실수 값이어야 한다. $\rightarrow$ 베타분포
데이터는 항상 0 또는 양수이어야 한다. $\rightarrow$ 로그정규분포, 감마분포, F분포, 카이제곱분포, 지수분포, 하프코시분포 등
데이터가 크기 제한이 없는 실수다. $\rightarrow$ 정규분포 또는 스튜던트 t분포, 코시분포, 라플라스분포 등

이 규칙에는 예외가 있을 수 있다. 예를 들어 항상 양수인 데이터인 경우에도 정규분포로 모형화가 가능하다면 정규분포를 사용할 수 있다. 정규분포와 스튜던트 t분포와 같이 둘 중 어느 것인지 구분하기 힘든 경우에는 뒤에서 설명할 정규성 검정이나 KS검정을 사용한다.

연습 문제 9.1.1

다음 코드로 사이킷런 패키지가 제공하는 보스턴 집값 데이터를 살펴보고 각각의 데이터의 확률분포 특성을 설명하라. 각 데이터에 적합한 확률분포가 존재한다면 어떤 확률분포인지도 설명하라.

from sklearn.datasets import load_boston
boston = load_boston()
dfX = pd.DataFrame(boston.data, columns=boston.feature_names)
dfy = pd.DataFrame(boston.target, columns=["MEDV"])
df = pd.concat([dfX, dfy], axis=1)

✏️

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

data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)
data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])
target = raw_df.values[1::2, 2]

dfX = pd.DataFrame(data) # Create DataFrame for features
dfy = pd.DataFrame(target, columns=["MEDV"]) # Create DataFrame for target
df = pd.concat([dfX, dfy], axis=1) # Concatenate features and target

# Calculate the number of rows and columns for subplots
num_cols = 4  # Number of columns in the subplot grid
num_rows = int(np.ceil(len(df.columns) / num_cols))  # Number of rows

# Create the subplots
fig, axes = plt.subplots(num_rows, num_cols, figsize=(15, 10))
axes = axes.flatten()  # Flatten the axes array for easier iteration

# title
title = ["CRIM","ZN","INDUS","CHAS","NOX","RM","AGE",
         "DIS","RAD","TAX","PTRATIO","B","LSTAT","MEDV"]


# Iterate through columns and plot on subplots
for i, column in enumerate(df.columns):
    sns.histplot(df[column], kde=True, bins=30, ax=axes[i])
    axes[i].set_title(f'Distribution of {title[i]}')

# Adjust layout and show the plot
plt.tight_layout()
plt.show()

📊

데이터의 확률분포 특성 분석 및 적합한 확률분포 제안

CRIM (범죄율)
- 분포: 지수분포 또는 로그 정규분포와 유사.
- 특징: 낮은 값에 데이터가 집중되고, 극단적으로 큰 값이 적은 비대칭 분포.
ZN (대형 주택 비율)
- 분포: 베르누이 분포 (0이 많은 경우) 또는 감마분포 (양의 값).
- 특징: 많은 0 값이 있고, 일부 지역에 양의 값.
INDUS (비소매 업종 비율)
- 분포: 정규분포 또는 베타분포.
- 특징: 비율로 제한된 범위에 있으나, 특정 중심에 데이터 집중.
CHAS (찰스강 인접 여부)
- 분포: 베르누이 분포 (이진값).
- 특징: 0 또는 1로만 이루어진 값.
NOX (일산화질소 농도)
- 분포: 베타분포.
- 특징: 0~1 사이에 존재하며, 농도 값이 특정 구간에 몰림.
RM (평균 방 개수)
- 분포: 정규분포.
- 특징: 평균값 근처에 집중된 연속형 분포.
AGE (1940년 이전 주택 비율)
- 분포: 베타분포.
- 특징: 0~100% 사이 비율 값.
DIS (고용센터 거리)
- 분포: 지수분포 또는 로그 정규분포.
- 특징: 짧은 거리에 집중되고, 긴 꼬리를 가짐.
RAD (고속도로 접근 용이성)
- 분포: 포아송 분포.
- 특징: 이산적인 값으로 이루어짐.
TAX (재산세율)
- 분포: 정규분포 또는 감마분포.
- 특징: 특정 구간에 집중된 연속값.
PTRATIO (학생-교사 비율)
- 분포: 정규분포.
- 특징: 학교 특성에 따라 제한된 구간.
B (흑인 거주 비율)
- 분포: 정규분포 또는 감마분포.
- 특징: 특정 값 주변으로 퍼져 있음.
LSTAT (저소득층 비율)
- 분포: 지수분포 또는 정규분포.
- 특징: 비율 값으로 데이터가 고르게 분포.
MEDV (집값, 타겟)
- 분포: 정규분포 (일부 데이터는 고정값으로 제한됨).
- 특징: 주택 시장 특성에 따라 특정 값에 몰림.

모수 추정 방법론

모수의 값으로 가장 가능성이 높은 하나의 숫자를 찾아내는 작업을 모수 추정(parameter estimation) 이라고 한다. 모수 추정 방법에는 다음과 같은 방법들이 있다. 이 절에서는 우선 모멘트 방법을 공부한다. 최대가능도 추정법과 베이즈 추정법은 다음 절에서 공부한다.

모멘트 방법
최대가능도 추정법
베이즈 추정법

모멘트 방법

모멘트 방법(method of moment)은 표본자료에 대한 표본모멘트가 확률분포의 이론적 모멘트와 같다고 가정하여 모수를 구한다. $N$ 은 데이터의 개수, $x_i$ 는 표본 데이터다.

\mu = \text{E}[X] \triangleq \bar{x} = \dfrac{1}{N} \sum_{i=1}^N x_i \tag{9.1.1}

\sigma^2 = \text{E}[(X-\mu)^2] \triangleq \bar{s}^2 = \dfrac{1}{N-1} \sum_{i=1}^N (x_i - \bar{x})^2 \tag{9.1.2}

예제: 베르누이분포의 모수 추정

모멘트 방법으로 베르누이 확률변수의 모수 $\mu$ 를 구하면 다음과 같다. 이 식에서 $N_1$ 은 1의 개수이다.

\text{E}[X] = \mu \triangleq \bar{x} = \dfrac{1}{N} \sum_{i=1}^N x_i = \dfrac{N_1}{N} \tag{9.1.3}

예제: 정규분포의 모수 추정

모멘트 방법으로 정규분포의 모수 $\mu$ , $\sigma^2$ 를 구하면 다음과 같다.

\text{E}[X] = \mu \triangleq \bar{x} = \dfrac{1}{N} \sum_{i=1}^N x_i \tag{9.1.4}

\text{E}[(X-\mu)^2] = \sigma^2 \triangleq s^2 = \dfrac{1}{N-1} \sum_{i=1}^N (x_i - \bar{x})^2 \tag{9.1.5}

예제: 베타 분포의 모수 추정

모멘트 방법으로 베타 분포의 모수 $a$ , $b$ 를 구하면 다음과 같다. 이 경우에는 모수와 모멘트 간의 관계를 이용하여 비선형 연립 방정식을 풀어야 한다.

\text{E}[X] = \dfrac{a}{a+b} \triangleq \bar{x} \tag{9.1.6}

\text{E}[(X-\mu)^2] = \dfrac{ab}{(a+b)^2(a+b+1)} \triangleq s^2 \tag{9.1.7}

a = \bar{x} \left( \frac{\bar{x} (1 - \bar{x})}{s^2} - 1 \right) \tag{9.1.8}

b = (1 - \bar{x}) \left( \frac{\bar{x} (1 - \bar{x})}{s^2} - 1 \right) \tag{9.1.9}

예를 들어 다음과 같은 데이터 100개가 있다고 하자. 값이 항상 0과 1사이에 있으므로 베타분포를 따른다고 가정하다. 사실 이 데이터는 $a=15, b=12$ 인 베타분포에서 생성한 것이다.

np.random.seed(0)
x = sp.stats.beta(15, 12).rvs(10000)

sns.distplot(x, kde=False, norm_hist=True)
plt.title("베타 분포를 따르는 표본의 히스토그램")
plt.show()

모멘트 방법으로 모수를 계산하면 원래의 모수와 비슷한 값을 구할 수 있다.

def estimate_beta(x):
    x_bar = x.mean()
    s2 = x.var()
    a = x_bar * (x_bar * (1 - x_bar) / s2 - 1)
    b = (1 - x_bar) * (x_bar * (1 - x_bar) / s2 - 1)
    return a, b

params = estimate_beta(x)
print(params)

(15.346682046700685, 12.2121537049535)

추정된 모수값으로 확률밀도분포를 그리면 히스토그램과 일치하는 것을 볼 수 있다.

xx = np.linspace(0, 1, 1000)
sns.distplot(x, kde=False, norm_hist=True)
plt.plot(xx, sp.stats.beta(params[0], params[1]).pdf(xx))
plt.xlim(0, 1)
plt.title("베타 분포를 따르는 표본의 히스토그램과 추정된 확률밀도함수")
plt.show()

연속확률분포의 히스토그램을 그릴 수 있는 시본의 distplot() 함수에는 사실 모수 추정기능이 포함되어 있다. fit인수로 사이파이의 확률변수 명령을 넣으면 이 명령을 사용하여 모수를 추정한 뒤에 해당 확률밀도함수 그래프를 히스토그램과 함께 보여준다.

sns.distplot(x, kde=False, norm_hist=True, fit=sp.stats.beta)
plt.xlim(0, 1)
plt.title("distplot으로 추정한 확률밀도함수")
plt.show()

연습 문제 9.1.2

위 연습 문제에서 나온 보스턴 집값 데이터 각각에 대해 시본의 distplot() 함수로 히스토그램을 그려라. 그리고 distplot() 함수의 모수 추정 기능을 사용하여 각각의 데이터에 적합한 확률분포의 확률밀도함수를 그려라. (범주형 데이터는 제외한다.)

✏️
(범주형 데이터는 제외한다.) 라는 조건을 늦게 봐서 억지로 코드를 짰더니, 범주형 데이터의 경우 결과가 잘 안나왔다.

import numpy as np
import pandas as pd
import scipy.stats as stats
import seaborn as sns
import matplotlib.pyplot as plt

# 예시 데이터셋 (df는 이미 존재한다고 가정)
# df = pd.DataFrame(...)

# 각 컬럼에 대해 적합한 분포 리스트
distributions = [
    stats.expon,      # CRIM: 지수분포
    stats.bernoulli,  # ZN: 베르누이 분포
    stats.norm,       # INDUS: 정규분포
    stats.bernoulli,  # CHAS: 베르누이 분포
    stats.beta,       # NOX: 베타분포
    stats.norm,       # RM: 정규분포
    stats.beta,       # AGE: 베타분포
    stats.lognorm,    # DIS: 로그 정규분포
    stats.poisson,    # RAD: 포아송 분포
    stats.gamma,      # TAX: 감마분포
    stats.norm,       # PTRATIO: 정규분포
    stats.gamma,      # B: 감마분포
    stats.expon,      # LSTAT: 지수분포
    stats.norm        # MEDV: 정규분포
]

# 서브플롯 설정
fig, axes = plt.subplots(7, 2, figsize=(14, 20))
axes = axes.flatten()

# 각 컬럼에 대해 분포 피팅 및 시각화
for i, (column, dist) in enumerate(zip(df.columns, distributions)):
    data = df[column]
    
    # 데이터의 히스토그램
    sns.histplot(data, kde=True, bins=30, ax=axes[i], color='lightblue', stat='density')
    
    # 이산형 분포 (베르누이, 포아송 등) 경우 파라미터 추정
    if dist == stats.bernoulli:
        # 베르누이 분포의 경우 성공 확률(p) 추정
        p = np.mean(data)
        params = (p,)
        # 베르누이 분포는 PMF를 사용해야 합니다.
        x = np.arange(0, 2)  # 베르누이 분포의 경우 가능한 값은 0, 1입니다.
        p_values = dist.pmf(x, *params)
    elif dist == stats.poisson:
        # 포아송 분포의 경우 평균(lambda) 추정
        lambda_ = np.mean(data)
        params = (lambda_,)
        # 포아송 분포는 PMF를 사용합니다.
        x = np.arange(0, np.max(data) + 1)
        p_values = dist.pmf(x, *params)
    else:
        # 연속형 분포의 경우 fit 사용
        params = dist.fit(data)
        # 연속형 분포는 PDF를 사용합니다.
        xmin, xmax = axes[i].get_xlim()
        x = np.linspace(xmin, xmax, 100)
        p_values = dist.pdf(x, *params)
    
    axes[i].plot(x, p_values, 'r-', lw=2, label=f'{dist.name} fit')
    axes[i].set_title(f'{column} - {dist.name} fit')
    axes[i].legend()

plt.tight_layout()
plt.show()