교통사고 데이터 분석 종합 스터디 가이드

# ============================================================
# [Step 1] Google Drive 마운트
# ============================================================
# Google Colab 환경에서 Google Drive에 저장된 데이터를
# 불러오기 위해 Drive를 마운트합니다.
# 실행 시 Google 계정 인증 팝업이 나타나며,
# 인증 완료 후 /content/drive/MyDrive 경로로 접근할 수 있습니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

# Google Drive 마운트 라이브러리 임포트
from google.colab import drive

# Google Drive를 /content/drive 경로에 마운트
drive.mount('/content/drive')

# 마운트 완료 확인 메시지 출력
print("Google Drive 마운트가 완료되었습니다.")
print("데이터 경로: /content/drive/MyDrive/Data/")

# ============================================================
# [Step 2] 데이터 불러오기
# ============================================================
# Google Drive의 Data 폴더에 저장된
# 교통사고 데이터셋(dataset_traffic_accident.csv)을 불러옵니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

# 데이터 분석 라이브러리 임포트
import pandas as pd

# 데이터 파일 경로 설정
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident.csv'

# CSV 파일 불러오기
df = pd.read_csv(data_path)

# 데이터 기본 정보 확인
print("=" * 60)
print("데이터 불러오기 완료")
print("=" * 60)
print(f"행 수: {df.shape[0]}")
print(f"열 수: {df.shape[1]}")
print("=" * 60)

# 상위 5개 행 미리보기
print("\n[상위 5개 행 미리보기]")
df.head()

# ============================================================
# [Step 3] 변수 항목 정보 확인
# ============================================================
# 데이터셋의 변수(컬럼) 정보를 상세히 확인합니다.
# - 컬럼명, 데이터 타입, 결측치 수, 고유값 수 등을 파악합니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

# ----------------------------------------------------------
# 1) 전체 컬럼명 목록 확인
# ----------------------------------------------------------
print("=" * 60)
print("[1] 전체 컬럼명 목록")
print("=" * 60)
for i, col in enumerate(df.columns, 1):
    print(f"  {i:>3}. {col}")
print(f"\n총 컬럼 수: {len(df.columns)}개")

# ----------------------------------------------------------
# 2) 데이터 타입 및 결측치 정보 (df.info())
# ----------------------------------------------------------
print("\n" + "=" * 60)
print("[2] 데이터 타입 및 결측치 정보")
print("=" * 60)
df.info()

# ----------------------------------------------------------
# 3) 기초 통계량 확인 (수치형 변수)
# ----------------------------------------------------------
print("\n" + "=" * 60)
print("[3] 수치형 변수 기초 통계량")
print("=" * 60)
df.describe()

# ============================================================
# [Step 4] 데이터셋 상세 검토
# ============================================================
# 각 변수의 고유값 분포, 결측치 현황, 이상치 등을
# 더 상세히 확인합니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

# ----------------------------------------------------------
# 1) 결측치 현황 요약
# ----------------------------------------------------------
print("=" * 60)
print("[1] 결측치 현황")
print("=" * 60)
missing = df.isnull().sum()  # 각 컬럼별 결측치 수 계산
missing_pct = (df.isnull().sum() / len(df) * 100).round(2)  # 결측치 비율(%) 계산
missing_df = pd.DataFrame({
    '결측치 수': missing,
    '결측치 비율(%)': missing_pct
})
# 결측치가 있는 변수만 필터링하여 출력
print(missing_df[missing_df['결측치 수'] > 0])
print(f"\n전체 변수 중 결측치 포함 변수: {(missing > 0).sum()}개")

# ----------------------------------------------------------
# 2) 범주형 변수 고유값 분포 확인
# ----------------------------------------------------------
print("\n" + "=" * 60)
print("[2] 범주형 변수 고유값 분포")
print("=" * 60)

# object 타입인 컬럼만 선택
cat_cols = df.select_dtypes(include='object').columns.tolist()

for col in cat_cols:
    print(f"\n--- {col} ---")
    print(f"고유값 수: {df[col].nunique()}개")
    print(df[col].value_counts())  # 각 고유값의 빈도수 출력

# ----------------------------------------------------------
# 3) 수치형 변수 고유값 및 이상치 확인
# ----------------------------------------------------------
print("\n" + "=" * 60)
print("[3] 수치형 변수 고유값 및 이상치 확인")
print("=" * 60)

# 수치형 컬럼 선택
num_cols = df.select_dtypes(include=['int64', 'float64']).columns.tolist()

for col in num_cols:
    print(f"\n--- {col} ---")
    print(f"고유값 수: {df[col].nunique()}개")
    print(f"최솟값: {df[col].min()}, 최댓값: {df[col].max()}")
    print(f"평균: {df[col].mean():.2f}, 중앙값: {df[col].median():.2f}")

    # IQR 기반 이상치 탐지
    Q1 = df[col].quantile(0.25)  # 1사분위수
    Q3 = df[col].quantile(0.75)  # 3사분위수
    IQR = Q3 - Q1                # 사분위 범위
    lower = Q1 - 1.5 * IQR       # 이상치 하한
    upper = Q3 + 1.5 * IQR       # 이상치 상한
    outliers = df[(df[col] < lower) | (df[col] > upper)]  # 이상치 필터링
    print(f"IQR 기반 이상치 수: {len(outliers)}개 (하한: {lower}, 상한: {upper})")

# ----------------------------------------------------------
# 4) 타겟 변수(Accident_Severity) 분포 확인
# ----------------------------------------------------------
print("\n" + "=" * 60)
print("[4] 타겟 변수(Accident_Severity) 분포")
print("=" * 60)
print(df['Accident_Severity'].value_counts())  # 각 클래스별 빈도수
print(f"\n사고 심각(1) 비율: {df['Accident_Severity'].mean() * 100:.2f}%")

# ============================================================
# [Step 5] 수치형 변수 기초통계 분석 및 표 정리
# ============================================================
# 수치형 변수에 대해 최소값, 평균, 중위값, 최대값, 표준편차,
# 왜곡도(Skewness), 결측치 수, IQR 기반 이상치 수를
# 하나의 표로 정리합니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np
from scipy import stats  # 왜곡도 계산을 위한 라이브러리

# ----------------------------------------------------------
# 1) 수치형 변수 선택
# ----------------------------------------------------------
num_cols = df.select_dtypes(include=['int64', 'float64']).columns.tolist()

# ----------------------------------------------------------
# 2) 각 수치형 변수별 기초통계 계산
# ----------------------------------------------------------
# 결과를 저장할 빈 리스트
result_list = []

for col in num_cols:
    # 기초통계량 계산
    col_min = df[col].min()                          # 최소값
    col_mean = df[col].mean()                        # 평균
    col_median = df[col].median()                    # 중위값
    col_max = df[col].max()                          # 최대값
    col_std = df[col].std()                          # 표준편차
    col_skew = df[col].skew()                        # 왜곡도 (Skewness)
    col_missing = df[col].isnull().sum()             # 결측치 수

    # IQR 기반 이상치 수 계산
    Q1 = df[col].quantile(0.25)                      # 1사분위수
    Q3 = df[col].quantile(0.75)                      # 3사분위수
    IQR = Q3 - Q1                                    # 사분위 범위
    lower_bound = Q1 - 1.5 * IQR                     # 이상치 하한
    upper_bound = Q3 + 1.5 * IQR                     # 이상치 상한
    col_outliers = df[(df[col] < lower_bound) | (df[col] > upper_bound)].shape[0]  # 이상치 수

    # 결과 리스트에 추가
    result_list.append({
        '변수명': col,
        '최소값': round(col_min, 2),
        '평균': round(col_mean, 2),
        '중위값': round(col_median, 2),
        '최대값': round(col_max, 2),
        '표준편차': round(col_std, 2),
        '왜곡도': round(col_skew, 2),
        '결측치': col_missing,
        '이상치': col_outliers
    })

# ----------------------------------------------------------
# 3) 표(DataFrame) 형태로 정리하여 출력
# ----------------------------------------------------------
stats_df = pd.DataFrame(result_list)
stats_df = stats_df.set_index('변수명')  # 변수명을 인덱스로 설정

print("=" * 80)
print("수치형 변수 기초통계 요약표")
print("=" * 80)
print(stats_df.to_string())

# ----------------------------------------------------------
# 4) 왜곡도 해석 기준 출력
# ----------------------------------------------------------
print("\n" + "=" * 80)
print("왜곡도(Skewness) 해석 기준")
print("=" * 80)
print("  |왜곡도| < 0.5  → 대체로 대칭 분포")
print("  0.5 ≤ |왜곡도| < 1.0 → 약간 비대칭")
print("  |왜곡도| ≥ 1.0  → 강한 비대칭 (치우침이 큼)")

# ----------------------------------------------------------
# 5) 변수별 해석 출력
# ----------------------------------------------------------
print("\n" + "=" * 80)
print("변수별 분석 결과 해석")
print("=" * 80)

for _, row in stats_df.iterrows():
    name = row.name  # 변수명 (인덱스)
    print(f"\n▶ {name}")

    # 평균과 중위값 비교 → 분포 치우침 방향 판단
    if row['평균'] > row['중위값']:
        direction = "오른쪽으로 치우침 (양의 왜곡)"
    elif row['평균'] < row['중위값']:
        direction = "왼쪽으로 치우침 (음의 왜곡)"
    else:
        direction = "대칭에 가까움"

    # 왜곡도 크기에 따른 비대칭 정도 판단
    abs_skew = abs(row['왜곡도'])
    if abs_skew < 0.5:
        skew_level = "대체로 대칭"
    elif abs_skew < 1.0:
        skew_level = "약간 비대칭"
    else:
        skew_level = "강한 비대칭"

    print(f"  - 범위: {row['최소값']} ~ {row['최대값']}, 평균: {row['평균']}, 중위값: {row['중위값']}")
    print(f"  - 분포: {direction} | 왜곡도 {row['왜곡도']} → {skew_level}")
    print(f"  - 결측치: {int(row['결측치'])}건, 이상치(IQR): {int(row['이상치'])}건")

# ============================================================
# [Step 6] 범주형 변수 기초통계 분석 및 표 정리
# ============================================================
# 범주형 변수에 대해 각 고유값의 빈도, 비율(%), 결측치 수를
# 변수별로 표 형태로 정리합니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

import pandas as pd

# ----------------------------------------------------------
# 1) 범주형 변수 선택
# ----------------------------------------------------------
cat_cols = df.select_dtypes(include='object').columns.tolist()

# ----------------------------------------------------------
# 2) 변수별 빈도/비율/결측치 표 출력
# ----------------------------------------------------------
for col in cat_cols:
    print("=" * 60)
    print(f"▶ {col}")
    print("=" * 60)

    # 각 고유값의 빈도수 계산 (결측치 제외)
    freq = df[col].value_counts()

    # 각 고유값의 비율(%) 계산 (전체 840건 기준)
    ratio = (df[col].value_counts() / len(df) * 100).round(2)

    # 결측치 수 계산
    missing_count = df[col].isnull().sum()

    # 표(DataFrame) 형태로 결합
    cat_table = pd.DataFrame({
        '빈도': freq,
        '비율(%)': ratio
    })

    # 합계 행 추가
    cat_table.loc['합계(결측 제외)'] = [
        cat_table['빈도'].sum(),       # 빈도 합계
        cat_table['비율(%)'].sum()     # 비율 합계
    ]

    print(cat_table.to_string())
    print(f"\n결측치: {missing_count}건 ({(missing_count / len(df) * 100):.2f}%)")
    print(f"고유값 수: {df[col].nunique()}개")
    print()

# ----------------------------------------------------------
# 3) 범주형 변수 요약 총괄표
# ----------------------------------------------------------
print("=" * 60)
print("범주형 변수 요약 총괄표")
print("=" * 60)

# 총괄 요약 정보를 리스트로 수집
summary_list = []
for col in cat_cols:
    missing = df[col].isnull().sum()                    # 결측치 수
    n_unique = df[col].nunique()                        # 고유값 수
    top_value = df[col].mode()[0]                       # 최빈값
    top_freq = df[col].value_counts().iloc[0]           # 최빈값 빈도
    top_ratio = round(top_freq / len(df) * 100, 2)     # 최빈값 비율(%)

    summary_list.append({
        '변수명': col,
        '고유값 수': n_unique,
        '최빈값': top_value,
        '최빈값 빈도': top_freq,
        '최빈값 비율(%)': top_ratio,
        '결측치 수': missing,
        '결측치 비율(%)': round(missing / len(df) * 100, 2)
    })

# DataFrame으로 변환 후 출력
summary_df = pd.DataFrame(summary_list).set_index('변수명')
print(summary_df.to_string())

# ============================================================
# [Step 7] 전처리 이슈 식별 및 해결 방안 정리
# ============================================================
# 앞서 수행한 기초통계 분석 결과를 바탕으로
# 전처리가 필요한 이슈를 식별하고, 해결 방안과
# 구현 순서를 표 형태로 정리합니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

import pandas as pd

# ----------------------------------------------------------
# 1) 전처리 이슈 식별 및 해결 방안 테이블 구성
# ----------------------------------------------------------

# 각 전처리 단계를 딕셔너리 리스트로 정의
preprocessing_plan = [
    {
        '구현 단계': '1단계',
        '이슈': '결측치',
        '변수': 'Weather',
        '이슈 설명': '30건(3.57%) 결측 — 범주형 변수이며 최빈값(Clear)이 58.33%로 지배적',
        '해결 방안': '최빈값(Clear)으로 대체'
    },
    {
        '구현 단계': '1단계',
        '이슈': '결측치',
        '변수': 'Traffic_Density',
        '이슈 설명': '30건(3.57%) 결측 — 0/1/2 세 값만 존재하는 순서형 변수, 중위값 1.0',
        '해결 방안': '중위값(1.0)으로 대체 후 정수 변환'
    },
    {
        '구현 단계': '1단계',
        '이슈': '결측치',
        '변수': 'Driver_Experience',
        '이슈 설명': '10건(1.19%) 결측 — 연속형 변수, 평균(19.25)과 중위값(19.0)이 유사하며 대칭 분포',
        '해결 방안': '중위값(19.0)으로 대체'
    },
    {
        '구현 단계': '2단계',
        '이슈': '이상치',
        '변수': 'Speed_Limit',
        '이슈 설명': 'IQR 기반 이상치 29건 — 최대값 300은 일반 제한속도 범위를 크게 초과 (왜곡도 3.29)',
        '해결 방안': 'IQR 상한(110) 초과 값을 상한값(110)으로 클리핑(Capping)'
    },
    {
        '구현 단계': '2단계',
        '이슈': '이상치',
        '변수': 'Driver_Age',
        '이슈 설명': '최솟값 5세 — 통계적 이상치는 아니나 운전 불가능 연령(도메인 이상치)',
        '해결 방안': '만 18세 미만 값을 18세로 클리핑(Capping) — 운전면허 취득 최소 연령 기준'
    },
    {
        '구현 단계': '3단계',
        '이슈': '데이터 타입',
        '변수': 'Traffic_Density',
        '이슈 설명': 'float64 타입이나 실제 값은 0, 1, 2 정수만 존재 — 순서형 범주 변수',
        '해결 방안': 'int 타입으로 변환'
    },
    {
        '구현 단계': '3단계',
        '이슈': '데이터 타입',
        '변수': 'Driver_Experience',
        '이슈 설명': 'float64 타입이나 결측치 대체 후 정수로 표현 가능',
        '해결 방안': 'int 타입으로 변환'
    },
    {
        '구현 단계': '4단계',
        '이슈': '전처리 검증',
        '변수': '전체 변수',
        '이슈 설명': '전처리 적용 후 결측치·이상치·데이터 타입이 올바르게 처리되었는지 확인 필요',
        '해결 방안': '결측치 재확인, 이상치 재검사, 데이터 타입 확인, 기초통계 재출력'
    },
    {
        '구현 단계': '5단계',
        '이슈': '데이터 저장',
        '변수': '전체 변수',
        '이슈 설명': '전처리 완료된 데이터셋을 원본과 분리하여 보존 필요',
        '해결 방안': 'Google Drive의 data 폴더에 dataset_traffic_accident_cleaned.csv로 저장'
    }
]

# ----------------------------------------------------------
# 2) DataFrame으로 변환 및 출력
# ----------------------------------------------------------
plan_df = pd.DataFrame(preprocessing_plan)

print("=" * 120)
print("전처리 이슈 식별 및 구현 순서")
print("=" * 120)

# 보기 좋게 전체 컬럼이 출력되도록 설정
pd.set_option('display.max_colwidth', None)  # 컬럼 내용 잘림 방지
pd.set_option('display.width', 120)          # 출력 폭 설정

print(plan_df.to_string(index=False))

# ----------------------------------------------------------
# 3) 전처리 전후 비교를 위한 현황 요약
# ----------------------------------------------------------
print("\n" + "=" * 120)
print("전처리 전 현황 요약")
print("=" * 120)
print(f"전체 데이터: {df.shape[0]}행 × {df.shape[1]}열")
print(f"결측치 보유 변수: 3개 (Weather 30건, Traffic_Density 30건, Driver_Experience 10건)")
print(f"이상치 보유 변수: 1개 (Speed_Limit 29건, IQR 기준)")
print(f"도메인 이상치 변수: 1개 (Driver_Age, 18세 미만)")
print(f"타입 변환 필요 변수: 2개 (Traffic_Density, Driver_Experience → float에서 int로)")

# ============================================================
# [Step 8] 전처리 1단계 — 결측치 처리
# ============================================================
# 결측치가 존재하는 3개 변수에 대해 적절한 대체값을 적용합니다.
#   - Weather (범주형, 30건) → 최빈값(Clear)으로 대체
#   - Traffic_Density (순서형, 30건) → 중위값(1.0)으로 대체
#   - Driver_Experience (연속형, 10건) → 중위값(19.0)으로 대체
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

# ----------------------------------------------------------
# 처리 전 결측치 현황 확인
# ----------------------------------------------------------
print("=" * 60)
print("[1단계] 결측치 처리")
print("=" * 60)
print("\n▶ 처리 전 결측치 현황:")
print(df.isnull().sum()[df.isnull().sum() > 0])

# ----------------------------------------------------------
# Weather: 최빈값(Clear)으로 대체
# ----------------------------------------------------------
weather_mode = df['Weather'].mode()[0]  # 최빈값 계산
df['Weather'] = df['Weather'].fillna(weather_mode)
print(f"\n✔ Weather → 최빈값('{weather_mode}')으로 {30}건 대체 완료")

# ----------------------------------------------------------
# Traffic_Density: 중위값(1.0)으로 대체
# ----------------------------------------------------------
td_median = df['Traffic_Density'].median()  # 중위값 계산
df['Traffic_Density'] = df['Traffic_Density'].fillna(td_median)
print(f"✔ Traffic_Density → 중위값({td_median})으로 {30}건 대체 완료")

# ----------------------------------------------------------
# Driver_Experience: 중위값(19.0)으로 대체
# ----------------------------------------------------------
de_median = df['Driver_Experience'].median()  # 중위값 계산
df['Driver_Experience'] = df['Driver_Experience'].fillna(de_median)
print(f"✔ Driver_Experience → 중위값({de_median})으로 {10}건 대체 완료")

# ----------------------------------------------------------
# 처리 후 결측치 현황 확인
# ----------------------------------------------------------
print(f"\n▶ 처리 후 전체 결측치 수: {df.isnull().sum().sum()}건")

# ============================================================
# [Step 9] 전처리 2단계 — 이상치 처리
# ============================================================
# 이상치가 확인된 2개 변수에 대해 클리핑(Capping)을 적용합니다.
#   - Speed_Limit → IQR 상한(110) 초과 값을 110으로 클리핑
#   - Driver_Age → 18세 미만 값을 18로 클리핑 (도메인 기준)
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

print("=" * 60)
print("[2단계] 이상치 처리 (클리핑)")
print("=" * 60)

# ----------------------------------------------------------
# Speed_Limit: IQR 기반 상한 클리핑
# ----------------------------------------------------------
# IQR 계산
Q1_speed = df['Speed_Limit'].quantile(0.25)           # 1사분위수
Q3_speed = df['Speed_Limit'].quantile(0.75)           # 3사분위수
IQR_speed = Q3_speed - Q1_speed                       # 사분위 범위
upper_speed = Q3_speed + 1.5 * IQR_speed              # 이상치 상한

# 처리 전 이상치 수 확인
outlier_count_before = (df['Speed_Limit'] > upper_speed).sum()
print(f"\n▶ Speed_Limit")
print(f"  - IQR 상한: {upper_speed}")
print(f"  - 처리 전 이상치 수: {outlier_count_before}건")
print(f"  - 처리 전 최대값: {df['Speed_Limit'].max()}")

# 상한 초과 값을 상한값으로 클리핑
df['Speed_Limit'] = df['Speed_Limit'].clip(upper=upper_speed)

print(f"  - 처리 후 최대값: {df['Speed_Limit'].max()}")
print(f"  ✔ {outlier_count_before}건을 {upper_speed}으로 클리핑 완료")

# ----------------------------------------------------------
# Driver_Age: 도메인 기준 하한 클리핑 (18세 미만 → 18세)
# ----------------------------------------------------------
# 처리 전 18세 미만 수 확인
under18_count = (df['Driver_Age'] < 18).sum()
print(f"\n▶ Driver_Age")
print(f"  - 도메인 하한: 18세 (운전면허 취득 최소 연령)")
print(f"  - 처리 전 18세 미만: {under18_count}건")
print(f"  - 처리 전 최솟값: {df['Driver_Age'].min()}")

# 18세 미만 값을 18로 클리핑
df['Driver_Age'] = df['Driver_Age'].clip(lower=18)

print(f"  - 처리 후 최솟값: {df['Driver_Age'].min()}")
print(f"  ✔ {under18_count}건을 18로 클리핑 완료")

# ============================================================
# [Step 10] 전처리 3단계 — 데이터 타입 변환
# ============================================================
# 결측치 대체 후 float64로 남아 있는 변수를 int로 변환합니다.
#   - Traffic_Density: float64 → int (0, 1, 2 순서형)
#   - Driver_Experience: float64 → int (운전 경력 년수)
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

print("=" * 60)
print("[3단계] 데이터 타입 변환")
print("=" * 60)

# ----------------------------------------------------------
# 처리 전 데이터 타입 확인
# ----------------------------------------------------------
print("\n▶ 처리 전 데이터 타입:")
print(f"  - Traffic_Density: {df['Traffic_Density'].dtype}")
print(f"  - Driver_Experience: {df['Driver_Experience'].dtype}")

# ----------------------------------------------------------
# Traffic_Density: float64 → int 변환
# ----------------------------------------------------------
df['Traffic_Density'] = df['Traffic_Density'].astype(int)

# ----------------------------------------------------------
# Driver_Experience: float64 → int 변환
# ----------------------------------------------------------
df['Driver_Experience'] = df['Driver_Experience'].astype(int)

# ----------------------------------------------------------
# 처리 후 데이터 타입 확인
# ----------------------------------------------------------
print("\n▶ 처리 후 데이터 타입:")
print(f"  - Traffic_Density: {df['Traffic_Density'].dtype}")
print(f"  - Driver_Experience: {df['Driver_Experience'].dtype}")
print("\n✔ 데이터 타입 변환 완료")

# ============================================================
# [Step 11] 전처리 4단계 — 전처리 검증
# ============================================================
# 1~3단계 전처리가 올바르게 적용되었는지 검증합니다.
#   - 결측치 잔존 여부 확인
#   - 이상치 재검사 (Speed_Limit, Driver_Age)
#   - 데이터 타입 확인
#   - 수치형 변수 기초통계 재출력
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np

print("=" * 60)
print("[4단계] 전처리 검증")
print("=" * 60)

# ----------------------------------------------------------
# 1) 결측치 잔존 여부 확인
# ----------------------------------------------------------
print("\n[검증 1] 결측치 잔존 여부")
print("-" * 40)
total_missing = df.isnull().sum().sum()  # 전체 결측치 합계
print(f"전체 결측치 수: {total_missing}건")
if total_missing == 0:
    print("✔ 결측치 처리 완료 — 잔존 결측치 없음")
else:
    print("✘ 결측치가 남아 있습니다:")
    print(df.isnull().sum()[df.isnull().sum() > 0])

# ----------------------------------------------------------
# 2) 이상치 재검사
# ----------------------------------------------------------
print("\n[검증 2] 이상치 재검사")
print("-" * 40)

# Speed_Limit: IQR 상한 초과 여부 확인
Q1_s = df['Speed_Limit'].quantile(0.25)
Q3_s = df['Speed_Limit'].quantile(0.75)
IQR_s = Q3_s - Q1_s
upper_s = Q3_s + 1.5 * IQR_s
outlier_speed = (df['Speed_Limit'] > upper_s).sum()
print(f"Speed_Limit — 범위: {df['Speed_Limit'].min()} ~ {df['Speed_Limit'].max()}, IQR 이상치: {outlier_speed}건")

# Driver_Age: 18세 미만 여부 확인
under18 = (df['Driver_Age'] < 18).sum()
print(f"Driver_Age — 범위: {df['Driver_Age'].min()} ~ {df['Driver_Age'].max()}, 18세 미만: {under18}건")

if outlier_speed == 0 and under18 == 0:
    print("✔ 이상치 처리 완료 — 잔존 이상치 없음")

# ----------------------------------------------------------
# 3) 데이터 타입 확인
# ----------------------------------------------------------
print("\n[검증 3] 데이터 타입 확인")
print("-" * 40)
print(df.dtypes)

# Traffic_Density, Driver_Experience가 int인지 확인
td_ok = df['Traffic_Density'].dtype in ['int64', 'int32']
de_ok = df['Driver_Experience'].dtype in ['int64', 'int32']
if td_ok and de_ok:
    print("\n✔ 데이터 타입 변환 완료 — Traffic_Density, Driver_Experience 모두 int")

# ----------------------------------------------------------
# 4) 수치형 변수 기초통계 재출력
# ----------------------------------------------------------
print("\n[검증 4] 전처리 후 수치형 변수 기초통계")
print("-" * 40)
print(df.describe().round(2).to_string())

# ============================================================
# [Step 12] 전처리 5단계 — 전처리 완료 데이터 저장
# ============================================================
# 전처리가 완료된 데이터셋을 원본과 분리하여
# Google Drive의 data 폴더에 새로운 CSV 파일로 저장합니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

# ----------------------------------------------------------
# 저장 경로 설정
# ----------------------------------------------------------
save_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'

# ----------------------------------------------------------
# CSV 파일로 저장 (index 제외)
# ----------------------------------------------------------
df.to_csv(save_path, index=False)

print("=" * 60)
print("[5단계] 전처리 완료 데이터 저장")
print("=" * 60)
print(f"\n✔ 저장 완료: {save_path}")
print(f"  - 행 수: {df.shape[0]}")
print(f"  - 열 수: {df.shape[1]}")
print(f"  - 결측치: {df.isnull().sum().sum()}건")

# ============================================================
# [Step 13] 전처리 전·후 비교 요약
# ============================================================
# 원본 데이터를 다시 불러와 전처리된 데이터와 비교합니다.
# 변경 사항을 항목별로 정리하여 표 형태로 출력합니다.
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np

# ----------------------------------------------------------
# 1) 원본 데이터 다시 불러오기 (비교용)
# ----------------------------------------------------------
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident.csv'
df_original = pd.read_csv(data_path)  # 원본 데이터
df_cleaned = df.copy()                # 전처리된 현재 데이터

# ----------------------------------------------------------
# 2) 전처리 전·후 결측치 비교
# ----------------------------------------------------------
print("=" * 70)
print("[비교 1] 결측치 변화")
print("=" * 70)

missing_compare = pd.DataFrame({
    '전처리 전': df_original.isnull().sum(),
    '전처리 후': df_cleaned.isnull().sum(),
    '변화': df_original.isnull().sum() - df_cleaned.isnull().sum()
})
# 변화가 있는 변수만 필터링
missing_changed = missing_compare[missing_compare['변화'] > 0]
print(missing_changed.to_string())
print(f"\n총 처리된 결측치: {missing_changed['변화'].sum()}건 → 0건")

# ----------------------------------------------------------
# 3) 전처리 전·후 수치형 기초통계 비교
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[비교 2] 수치형 변수 기초통계 변화")
print("=" * 70)

num_cols = ['Traffic_Density', 'Speed_Limit', 'Driver_Age',
            'Driver_Experience', 'Accident_Severity']

for col in num_cols:
    # 전처리 전·후 값이 달라진 변수만 상세 출력
    before_min = df_original[col].min()
    after_min = df_cleaned[col].min()
    before_max = df_original[col].max()
    after_max = df_cleaned[col].max()
    before_mean = df_original[col].mean()
    after_mean = df_cleaned[col].mean()
    before_std = df_original[col].std()
    after_std = df_cleaned[col].std()

    # 변화 여부 판단 (소수점 2자리 기준)
    changed = (round(before_min, 2) != round(after_min, 2) or
               round(before_max, 2) != round(after_max, 2) or
               round(before_mean, 2) != round(after_mean, 2))

    if changed:
        print(f"\n▶ {col} — 변화 있음")
        compare_stats = pd.DataFrame({
            '항목': ['최소값', '평균', '중위값', '최대값', '표준편차'],
            '전처리 전': [
                round(before_min, 2),
                round(before_mean, 2),
                round(df_original[col].median(), 2),
                round(before_max, 2),
                round(before_std, 2)
            ],
            '전처리 후': [
                round(after_min, 2),
                round(after_mean, 2),
                round(df_cleaned[col].median(), 2),
                round(after_max, 2),
                round(after_std, 2)
            ]
        })
        compare_stats['변화'] = compare_stats['전처리 후'] - compare_stats['전처리 전']
        print(compare_stats.to_string(index=False))
    else:
        print(f"\n▶ {col} — 변화 없음")

# ----------------------------------------------------------
# 4) 전처리 전·후 데이터 타입 비교
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[비교 3] 데이터 타입 변화")
print("=" * 70)

dtype_compare = pd.DataFrame({
    '전처리 전': df_original.dtypes.astype(str),
    '전처리 후': df_cleaned.dtypes.astype(str)
})
# 타입이 변경된 변수만 필터링
dtype_changed = dtype_compare[dtype_compare['전처리 전'] != dtype_compare['전처리 후']]
if len(dtype_changed) > 0:
    print(dtype_changed.to_string())
else:
    print("변경된 데이터 타입 없음")

# ----------------------------------------------------------
# 5) 전체 변경 사항 종합 요약표
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[종합] 전처리 변경 사항 요약표")
print("=" * 70)

summary_data = [
    {
        '단계': '1단계',
        '처리 유형': '결측치 대체',
        '변수': 'Weather',
        '전처리 전': f'결측 {df_original["Weather"].isnull().sum()}건',
        '전처리 후': f'결측 {df_cleaned["Weather"].isnull().sum()}건',
        '처리 내용': f'최빈값({df_original["Weather"].mode()[0]})으로 대체'
    },
    {
        '단계': '1단계',
        '처리 유형': '결측치 대체',
        '변수': 'Traffic_Density',
        '전처리 전': f'결측 {df_original["Traffic_Density"].isnull().sum()}건',
        '전처리 후': f'결측 {df_cleaned["Traffic_Density"].isnull().sum()}건',
        '처리 내용': f'중위값({df_original["Traffic_Density"].median()})으로 대체'
    },
    {
        '단계': '1단계',
        '처리 유형': '결측치 대체',
        '변수': 'Driver_Experience',
        '전처리 전': f'결측 {df_original["Driver_Experience"].isnull().sum()}건',
        '전처리 후': f'결측 {df_cleaned["Driver_Experience"].isnull().sum()}건',
        '처리 내용': f'중위값({df_original["Driver_Experience"].median()})으로 대체'
    },
    {
        '단계': '2단계',
        '처리 유형': '이상치 클리핑',
        '변수': 'Speed_Limit',
        '전처리 전': f'범위: {df_original["Speed_Limit"].min()}~{df_original["Speed_Limit"].max()}',
        '전처리 후': f'범위: {df_cleaned["Speed_Limit"].min()}~{df_cleaned["Speed_Limit"].max()}',
        '처리 내용': 'IQR 상한(110) 초과 → 110으로 클리핑'
    },
    {
        '단계': '2단계',
        '처리 유형': '이상치 클리핑',
        '변수': 'Driver_Age',
        '전처리 전': f'범위: {df_original["Driver_Age"].min()}~{df_original["Driver_Age"].max()}',
        '전처리 후': f'범위: {df_cleaned["Driver_Age"].min()}~{df_cleaned["Driver_Age"].max()}',
        '처리 내용': '18세 미만 → 18로 클리핑'
    },
    {
        '단계': '3단계',
        '처리 유형': '타입 변환',
        '변수': 'Traffic_Density',
        '전처리 전': str(df_original['Traffic_Density'].dtype),
        '전처리 후': str(df_cleaned['Traffic_Density'].dtype),
        '처리 내용': 'float64 → int 변환'
    },
    {
        '단계': '3단계',
        '처리 유형': '타입 변환',
        '변수': 'Driver_Experience',
        '전처리 전': str(df_original['Driver_Experience'].dtype),
        '전처리 후': str(df_cleaned['Driver_Experience'].dtype),
        '처리 내용': 'float64 → int 변환'
    }
]

summary_df = pd.DataFrame(summary_data)
pd.set_option('display.max_colwidth', None)
pd.set_option('display.width', 150)
print(summary_df.to_string(index=False))

# ----------------------------------------------------------
# 6) 전체 데이터 형태 비교
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[최종] 데이터셋 크기 비교")
print("=" * 70)
print(f"전처리 전: {df_original.shape[0]}행 × {df_original.shape[1]}열")
print(f"전처리 후: {df_cleaned.shape[0]}행 × {df_cleaned.shape[1]}열")
print("✔ 행 삭제 없이 전체 840건 유지")

# ============================================================
# [Step 14] 변수 간 상관관계 분석
# ============================================================
# 전처리 완료된 데이터셋(dataset_traffic_accident_cleaned.csv)을
# 기반으로 변수 간 상관관계를 분석합니다.
#   - 수치형 변수 간 피어슨 상관계수 행렬
#   - 상관관계 히트맵 시각화
#   - 타겟 변수(Accident_Severity)와의 상관관계 정리
#   - 범주형 변수와 타겟 변수 간 관계 분석 (카이제곱 검정)
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

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

# ----------------------------------------------------------
# 0) 한글 폰트 설정 (Colab 환경)
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'  # Colab 기본 폰트
plt.rcParams['axes.unicode_minus'] = False     # 마이너스 기호 깨짐 방지

# ----------------------------------------------------------
# 1) 전처리 완료 데이터 불러오기
# ----------------------------------------------------------
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'
df = pd.read_csv(data_path)

print("=" * 70)
print("[Step 14] 변수 간 상관관계 분석")
print("=" * 70)
print(f"데이터: {df.shape[0]}행 × {df.shape[1]}열\n")

# ----------------------------------------------------------
# 2) 수치형 변수 피어슨 상관계수 행렬
# ----------------------------------------------------------
num_cols = ['Traffic_Density', 'Speed_Limit', 'Driver_Age',
            'Driver_Experience', 'Accident_Severity']

corr_matrix = df[num_cols].corr().round(4)  # 피어슨 상관계수 계산

print("=" * 70)
print("[1] 수치형 변수 피어슨 상관계수 행렬")
print("=" * 70)
print(corr_matrix.to_string())

# ----------------------------------------------------------
# 3) 상관계수 해석 기준 출력
# ----------------------------------------------------------
print("\n※ 상관계수 해석 기준:")
print("  |r| < 0.1  → 무시할 수준 (거의 무관)")
print("  0.1 ≤ |r| < 0.3 → 약한 상관")
print("  0.3 ≤ |r| < 0.5 → 보통 상관")
print("  0.5 ≤ |r| < 0.7 → 강한 상관")
print("  |r| ≥ 0.7  → 매우 강한 상관")

# ----------------------------------------------------------
# 4) 타겟 변수(Accident_Severity)와의 상관관계 정리
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[2] Accident_Severity(타겟)와 수치형 변수 간 상관관계")
print("=" * 70)

# 타겟과의 상관계수를 절대값 기준 내림차순 정렬
target_corr = corr_matrix['Accident_Severity'].drop('Accident_Severity')
target_corr_abs = target_corr.abs().sort_values(ascending=False)

# 상관 강도 해석 함수
def interpret_corr(r):
    """상관계수 절대값에 따른 해석 문자열 반환"""
    abs_r = abs(r)
    if abs_r < 0.1:
        return "무시할 수준"
    elif abs_r < 0.3:
        return "약한 상관"
    elif abs_r < 0.5:
        return "보통 상관"
    elif abs_r < 0.7:
        return "강한 상관"
    else:
        return "매우 강한 상관"

# 방향 해석 함수
def interpret_direction(r):
    """상관계수 부호에 따른 방향 문자열 반환"""
    if r > 0:
        return "양(+)의 상관"
    elif r < 0:
        return "음(-)의 상관"
    else:
        return "무상관"

# 결과 테이블 생성
target_result = []
for col in target_corr_abs.index:
    r = target_corr[col]
    target_result.append({
        '변수': col,
        '상관계수': round(r, 4),
        '|상관계수|': round(abs(r), 4),
        '방향': interpret_direction(r),
        '강도': interpret_corr(r)
    })

target_df = pd.DataFrame(target_result)
print(target_df.to_string(index=False))

# ----------------------------------------------------------
# 5) 상관관계 히트맵 시각화
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[3] 상관관계 히트맵 시각화")
print("=" * 70)

fig, ax = plt.subplots(figsize=(10, 8))

# 히트맵 생성 (상관계수 값 표시, 색상 맵 적용)
sns.heatmap(
    corr_matrix,
    annot=True,             # 셀에 상관계수 값 표시
    fmt='.4f',              # 소수점 4자리까지 표시
    cmap='RdBu_r',          # 빨강(양) ~ 파랑(음) 색상 맵
    center=0,               # 0을 기준으로 색상 대칭
    vmin=-1, vmax=1,        # 색상 범위 -1 ~ 1
    square=True,            # 정사각형 셀
    linewidths=0.5,         # 셀 구분선 두께
    cbar_kws={'shrink': 0.8, 'label': 'Pearson Correlation'}
)

ax.set_title('Pearson Correlation Heatmap (Numeric Variables)',
             fontsize=14, fontweight='bold', pad=15)
ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha='right')
ax.set_yticklabels(ax.get_yticklabels(), rotation=0)

plt.tight_layout()
plt.show()

print("✔ 히트맵 출력 완료")

# ----------------------------------------------------------
# 6) 범주형 변수와 타겟 간 관계 (카이제곱 검정)
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[4] 범주형 변수와 Accident_Severity 간 관계 (카이제곱 검정)")
print("=" * 70)

cat_cols = ['Weather', 'Road_Type', 'Time_of_Day',
            'Road_Condition', 'Vehicle_Type', 'Road_Light_Condition']

chi2_results = []
for col in cat_cols:
    # 교차표(크로스탭) 생성
    contingency = pd.crosstab(df[col], df['Accident_Severity'])

    # 카이제곱 검정 수행
    chi2, p_value, dof, expected = stats.chi2_contingency(contingency)

    # 크래머의 V 계산 (효과 크기)
    n = contingency.sum().sum()              # 전체 관측수
    min_dim = min(contingency.shape) - 1     # min(행수, 열수) - 1
    cramers_v = np.sqrt(chi2 / (n * min_dim))  # 크래머의 V 공식

    # 유의성 판단 (유의수준 0.05 기준)
    significance = "유의함 (p < 0.05)" if p_value < 0.05 else "유의하지 않음 (p ≥ 0.05)"

    # 효과 크기 해석
    def interpret_cramers_v(v):
        """크래머의 V 값에 따른 효과 크기 해석"""
        if v < 0.1:
            return "약함"
        elif v < 0.3:
            return "보통"
        elif v < 0.5:
            return "강함"
        else:
            return "매우 강함"

    chi2_results.append({
        '변수': col,
        '카이제곱(χ²)': round(chi2, 4),
        'p-value': f"{p_value:.6f}" if p_value >= 0.000001 else f"{p_value:.2e}",
        '자유도': dof,
        "Cramer's V": round(cramers_v, 4),
        '효과 크기': interpret_cramers_v(cramers_v),
        '유의성': significance
    })

chi2_df = pd.DataFrame(chi2_results)
pd.set_option('display.max_colwidth', None)
pd.set_option('display.width', 150)
print(chi2_df.to_string(index=False))

# ----------------------------------------------------------
# 7) 범주형 변수별 사고 심각도 비율표
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("[5] 범주형 변수별 Accident_Severity=1 비율")
print("=" * 70)

for col in cat_cols:
    print(f"\n▶ {col}")
    # 각 범주별 사고 심각도 평균 (= 심각 사고 비율)
    severity_rate = df.groupby(col)['Accident_Severity'].agg(
        ['count', 'sum', 'mean']
    ).round(4)
    severity_rate.columns = ['전체 건수', '심각(1) 건수', '심각 비율']
    severity_rate = severity_rate.sort_values('심각 비율', ascending=False)
    severity_rate['심각 비율(%)'] = (severity_rate['심각 비율'] * 100).round(2)
    print(severity_rate[['전체 건수', '심각(1) 건수', '심각 비율(%)']].to_string())

# ============================================================
# [Step 15] EDA 시각화 Part 1 — 핵심 변수 분포
# ============================================================
# 상관관계 분석에서 타겟(Accident_Severity)과 유의미한
# 관계가 확인된 핵심 변수들의 분포를 시각화합니다.
#   - 핵심 수치형: Traffic_Density, Speed_Limit
#   - 핵심 범주형: Weather, Road_Condition
#   - 타겟 변수: Accident_Severity (클래스 분포)
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

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

# ----------------------------------------------------------
# 0) 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 100

# 전처리 완료 데이터 불러오기
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'
df = pd.read_csv(data_path)

# 색상 팔레트 정의
severity_colors = {0: '#3498db', 1: '#e74c3c'}  # 파랑: 비심각, 빨강: 심각
severity_labels = {0: 'Non-Severe (0)', 1: 'Severe (1)'}

# ============================================================
# [시각화 1] 타겟 변수 클래스 분포
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# (1-1) 타겟 변수 바 차트
counts = df['Accident_Severity'].value_counts().sort_index()
bars = axes[0].bar(
    [severity_labels[i] for i in counts.index],
    counts.values,
    color=[severity_colors[i] for i in counts.index],
    edgecolor='black', linewidth=0.5
)
# 바 위에 건수 및 비율 표시
for bar, count in zip(bars, counts.values):
    pct = count / len(df) * 100
    axes[0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 10,
                 f'{count}\n({pct:.1f}%)', ha='center', va='bottom', fontsize=11)
axes[0].set_title('Accident Severity Distribution', fontsize=13, fontweight='bold')
axes[0].set_ylabel('Count')
axes[0].set_ylim(0, max(counts.values) * 1.25)

# (1-2) 타겟 변수 파이 차트
axes[1].pie(
    counts.values,
    labels=[severity_labels[i] for i in counts.index],
    colors=[severity_colors[i] for i in counts.index],
    autopct='%1.1f%%', startangle=90,
    explode=[0, 0.05], shadow=True,
    textprops={'fontsize': 11}
)
axes[1].set_title('Accident Severity Proportion', fontsize=13, fontweight='bold')

plt.suptitle('Figure 1. Target Variable (Accident_Severity) Distribution',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# ============================================================
# [시각화 2] 핵심 범주형 변수 — Weather & Road_Condition
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# (2-1) Weather별 사고 심각도 비율
weather_order = df.groupby('Weather')['Accident_Severity'].mean().sort_values(ascending=False).index
weather_ct = pd.crosstab(df['Weather'], df['Accident_Severity'], normalize='index') * 100
weather_ct = weather_ct.loc[weather_order]

weather_ct.plot(
    kind='barh', stacked=True, ax=axes[0],
    color=[severity_colors[0], severity_colors[1]],
    edgecolor='black', linewidth=0.5
)
# 각 바에 심각 사고 비율(%) 표시
for i, (idx, row) in enumerate(weather_ct.iterrows()):
    severe_pct = row[1]
    axes[0].text(severe_pct / 2 + row[0], i, f'{severe_pct:.1f}%',
                 ha='center', va='center', fontsize=10, fontweight='bold', color='white')
axes[0].set_title('Severity Rate by Weather', fontsize=13, fontweight='bold')
axes[0].set_xlabel('Percentage (%)')
axes[0].set_ylabel('')
axes[0].legend(['Non-Severe (0)', 'Severe (1)'], loc='lower right')

# (2-2) Road_Condition별 사고 심각도 비율
road_order = df.groupby('Road_Condition')['Accident_Severity'].mean().sort_values(ascending=False).index
road_ct = pd.crosstab(df['Road_Condition'], df['Accident_Severity'], normalize='index') * 100
road_ct = road_ct.loc[road_order]

road_ct.plot(
    kind='barh', stacked=True, ax=axes[1],
    color=[severity_colors[0], severity_colors[1]],
    edgecolor='black', linewidth=0.5
)
for i, (idx, row) in enumerate(road_ct.iterrows()):
    severe_pct = row[1]
    axes[1].text(severe_pct / 2 + row[0], i, f'{severe_pct:.1f}%',
                 ha='center', va='center', fontsize=10, fontweight='bold', color='white')
axes[1].set_title('Severity Rate by Road Condition', fontsize=13, fontweight='bold')
axes[1].set_xlabel('Percentage (%)')
axes[1].set_ylabel('')
axes[1].legend(['Non-Severe (0)', 'Severe (1)'], loc='lower right')

plt.suptitle('Figure 2. Key Categorical Variables vs Accident Severity',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# ============================================================
# [시각화 3] 핵심 수치형 변수 — Traffic_Density & Speed_Limit
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# (3-1) Traffic_Density별 사고 심각도 비율
td_ct = pd.crosstab(df['Traffic_Density'], df['Accident_Severity'])
td_pct = pd.crosstab(df['Traffic_Density'], df['Accident_Severity'], normalize='index') * 100

td_labels = {0: 'Low (0)', 1: 'Medium (1)', 2: 'High (2)'}
x_pos = range(len(td_ct.index))

# 그룹 바 차트
width = 0.35
bars0 = axes[0].bar([x - width/2 for x in x_pos], td_ct[0], width,
                     label='Non-Severe (0)', color=severity_colors[0], edgecolor='black', linewidth=0.5)
bars1 = axes[0].bar([x + width/2 for x in x_pos], td_ct[1], width,
                     label='Severe (1)', color=severity_colors[1], edgecolor='black', linewidth=0.5)

# 각 바 위에 심각 비율 표시
for x, density in zip(x_pos, td_ct.index):
    severe_pct = td_pct.loc[density, 1]
    total = td_ct.loc[density].sum()
    axes[0].text(x, total * 0.55, f'Severe\n{severe_pct:.1f}%',
                 ha='center', va='bottom', fontsize=10, fontweight='bold')

axes[0].set_xticks(x_pos)
axes[0].set_xticklabels([td_labels.get(i, str(i)) for i in td_ct.index])
axes[0].set_title('Severity by Traffic Density', fontsize=13, fontweight='bold')
axes[0].set_xlabel('Traffic Density')
axes[0].set_ylabel('Count')
axes[0].legend()

# (3-2) Speed_Limit 분포 — 사고 심각도별
for sev in [0, 1]:
    subset = df[df['Accident_Severity'] == sev]['Speed_Limit']
    axes[1].hist(subset, bins=15, alpha=0.6,
                 label=severity_labels[sev], color=severity_colors[sev],
                 edgecolor='black', linewidth=0.5)
    # 각 그룹의 평균선 표시
    axes[1].axvline(subset.mean(), color=severity_colors[sev],
                    linestyle='--', linewidth=2, alpha=0.8)
    axes[1].text(subset.mean() + 1, axes[1].get_ylim()[1] * (0.9 - sev * 0.1),
                 f'Mean: {subset.mean():.1f}', color=severity_colors[sev],
                 fontsize=10, fontweight='bold')

axes[1].set_title('Speed Limit Distribution by Severity', fontsize=13, fontweight='bold')
axes[1].set_xlabel('Speed Limit')
axes[1].set_ylabel('Count')
axes[1].legend()

plt.suptitle('Figure 3. Key Numeric Variables vs Accident Severity',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

print("✔ Part 1 시각화 완료 (Figure 1~3)")

# ============================================================
# [Step 16] EDA 시각화 Part 2 — 변수 간 상호작용
# ============================================================
# 타겟과 유의미한 관계가 있는 핵심 변수들 사이의
# 상호작용(Interaction)을 시각화합니다.
#   - Weather × Road_Condition 교차 분석
#   - Traffic_Density × Speed_Limit 상호작용
#   - Driver_Age × Driver_Experience 다중공선성 확인
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

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

# ----------------------------------------------------------
# 0) 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 100

severity_colors = {0: '#3498db', 1: '#e74c3c'}

# ============================================================
# [시각화 4] Weather × Road_Condition 교차 히트맵
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# (4-1) Weather × Road_Condition 교차 빈도
cross_count = pd.crosstab(df['Weather'], df['Road_Condition'])

# 빈도가 높은 순서로 정렬
weather_order = ['Rainy', 'Stormy', 'Snowy', 'Foggy', 'Clear']
road_order = ['Wet', 'Icy', 'Under Construction', 'Dry']
cross_count = cross_count.reindex(index=weather_order, columns=road_order)

sns.heatmap(
    cross_count, annot=True, fmt='d', cmap='YlOrRd',
    linewidths=0.5, ax=axes[0],
    cbar_kws={'label': 'Count'}
)
axes[0].set_title('Weather × Road Condition (Frequency)',
                   fontsize=13, fontweight='bold')
axes[0].set_xlabel('Road Condition')
axes[0].set_ylabel('Weather')

# (4-2) Weather × Road_Condition별 심각 사고 비율 히트맵
cross_severity = df.groupby(['Weather', 'Road_Condition'])['Accident_Severity'].mean() * 100
cross_severity = cross_severity.unstack()
cross_severity = cross_severity.reindex(index=weather_order, columns=road_order)

sns.heatmap(
    cross_severity, annot=True, fmt='.1f', cmap='RdYlGn_r',
    linewidths=0.5, ax=axes[1],
    vmin=0, vmax=100,
    cbar_kws={'label': 'Severe Rate (%)'}
)
axes[1].set_title('Weather × Road Condition (Severe Rate %)',
                   fontsize=13, fontweight='bold')
axes[1].set_xlabel('Road Condition')
axes[1].set_ylabel('Weather')

plt.suptitle('Figure 4. Interaction: Weather × Road Condition',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# ============================================================
# [시각화 5] Traffic_Density × Speed_Limit 상호작용
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

td_labels = {0: 'Low (0)', 1: 'Medium (1)', 2: 'High (2)'}
td_colors = {0: '#2ecc71', 1: '#f39c12', 2: '#e74c3c'}

# (5-1) Traffic_Density별 Speed_Limit 박스플롯 (심각도 구분)
# 심각도별로 분리하여 Traffic_Density × Speed_Limit 패턴 확인
df_plot = df.copy()
df_plot['TD_Label'] = df_plot['Traffic_Density'].map(td_labels)
df_plot['Severity_Label'] = df_plot['Accident_Severity'].map(
    {0: 'Non-Severe (0)', 1: 'Severe (1)'}
)

sns.boxplot(
    data=df_plot,
    x='TD_Label', y='Speed_Limit', hue='Severity_Label',
    palette=[severity_colors[0], severity_colors[1]],
    ax=axes[0], order=['Low (0)', 'Medium (1)', 'High (2)']
)
axes[0].set_title('Speed Limit by Traffic Density & Severity',
                   fontsize=13, fontweight='bold')
axes[0].set_xlabel('Traffic Density')
axes[0].set_ylabel('Speed Limit')
axes[0].legend(title='Severity')

# (5-2) Traffic_Density × Speed_Limit 구간별 심각 사고 비율
# Speed_Limit을 3개 구간으로 분할
df_plot['Speed_Group'] = pd.cut(
    df_plot['Speed_Limit'],
    bins=[0, 50, 80, 110],
    labels=['Low (≤50)', 'Medium (51-80)', 'High (81-110)']
)

# 교차표 생성: Traffic_Density × Speed_Group별 심각 비율
interaction_rate = df_plot.groupby(
    ['Traffic_Density', 'Speed_Group']
)['Accident_Severity'].mean() * 100

interaction_pivot = interaction_rate.unstack(level=0)
interaction_pivot.columns = [td_labels.get(c, c) for c in interaction_pivot.columns]

interaction_pivot.plot(
    kind='bar', ax=axes[1],
    color=[td_colors[0], td_colors[1], td_colors[2]],
    edgecolor='black', linewidth=0.5
)
# 각 바 위에 비율 표시
for container in axes[1].containers:
    axes[1].bar_label(container, fmt='%.1f%%', fontsize=8, padding=2)

axes[1].set_title('Severe Rate by Speed Group × Traffic Density',
                   fontsize=13, fontweight='bold')
axes[1].set_xlabel('Speed Limit Group')
axes[1].set_ylabel('Severe Accident Rate (%)')
axes[1].set_xticklabels(axes[1].get_xticklabels(), rotation=0)
axes[1].legend(title='Traffic Density')
axes[1].set_ylim(0, 100)

plt.suptitle('Figure 5. Interaction: Traffic Density × Speed Limit',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# ============================================================
# [시각화 6] Driver_Age × Driver_Experience — 다중공선성 확인
# ============================================================
fig, ax = plt.subplots(figsize=(10, 7))

# 산점도: 심각도별 색상 구분
for sev in [0, 1]:
    subset = df[df['Accident_Severity'] == sev]
    ax.scatter(
        subset['Driver_Age'], subset['Driver_Experience'],
        c=severity_colors[sev],
        label=f'{"Severe (1)" if sev == 1 else "Non-Severe (0)"}',
        alpha=0.4, s=30, edgecolors='white', linewidth=0.3
    )

# 회귀선 추가
z = np.polyfit(df['Driver_Age'], df['Driver_Experience'], 1)  # 1차 회귀
p = np.poly1d(z)
age_range = np.linspace(df['Driver_Age'].min(), df['Driver_Age'].max(), 100)
ax.plot(age_range, p(age_range), '--', color='black', linewidth=2, alpha=0.7,
        label=f'Regression (r = {df["Driver_Age"].corr(df["Driver_Experience"]):.4f})')

ax.set_title('Figure 6. Driver Age vs Experience (Multicollinearity Check)',
             fontsize=14, fontweight='bold')
ax.set_xlabel('Driver Age', fontsize=12)
ax.set_ylabel('Driver Experience (years)', fontsize=12)
ax.legend(fontsize=10)

plt.tight_layout()
plt.show()

print("✔ Part 2 시각화 완료 (Figure 4~6)")

# ============================================================
# [Step 17] EDA 시각화 Part 3 — 종합 대시보드
# ============================================================
# 타겟 변수를 중심으로 모든 핵심 변수를 한눈에 비교할 수 있는
# 종합 시각화를 생성합니다.
#   - 전체 변수 중요도 바 차트 (상관계수 + Cramer's V)
#   - 핵심 4개 변수의 심각 사고 조건 종합
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

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

# ----------------------------------------------------------
# 0) 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 100

# ============================================================
# [시각화 7] 전체 변수 — 타겟과의 연관 강도 종합 비교
# ============================================================

# 수치형 변수: 피어슨 상관계수의 절대값
num_associations = {
    'Traffic_Density': abs(df['Traffic_Density'].corr(df['Accident_Severity'])),
    'Speed_Limit': abs(df['Speed_Limit'].corr(df['Accident_Severity'])),
    'Driver_Age': abs(df['Driver_Age'].corr(df['Accident_Severity'])),
    'Driver_Experience': abs(df['Driver_Experience'].corr(df['Accident_Severity']))
}

# 범주형 변수: Cramer's V 계산
cat_cols = ['Weather', 'Road_Type', 'Time_of_Day',
            'Road_Condition', 'Vehicle_Type', 'Road_Light_Condition']
cat_associations = {}

for col in cat_cols:
    contingency = pd.crosstab(df[col], df['Accident_Severity'])
    chi2, p_val, dof, expected = stats.chi2_contingency(contingency)
    n = contingency.sum().sum()
    min_dim = min(contingency.shape) - 1
    cramers_v = np.sqrt(chi2 / (n * min_dim))
    cat_associations[col] = cramers_v

# 전체 합치기
all_associations = {**num_associations, **cat_associations}
assoc_df = pd.DataFrame({
    'Variable': list(all_associations.keys()),
    'Association': list(all_associations.values()),
    'Type': (['Numeric'] * len(num_associations) +
             ['Categorical'] * len(cat_associations))
}).sort_values('Association', ascending=True)

# 바 차트 생성
fig, ax = plt.subplots(figsize=(12, 7))

colors = ['#e74c3c' if t == 'Categorical' else '#3498db'
          for t in assoc_df['Type']]

bars = ax.barh(
    assoc_df['Variable'], assoc_df['Association'],
    color=colors, edgecolor='black', linewidth=0.5, height=0.6
)

# 각 바 옆에 수치 표시
for bar, val in zip(bars, assoc_df['Association']):
    ax.text(bar.get_width() + 0.01, bar.get_y() + bar.get_height()/2,
            f'{val:.4f}', ha='left', va='center', fontsize=10, fontweight='bold')

# 유의성 기준선
ax.axvline(x=0.1, color='gray', linestyle=':', linewidth=1, alpha=0.7)
ax.text(0.105, len(assoc_df) - 0.5, 'Weak (0.1)', fontsize=9, color='gray')
ax.axvline(x=0.3, color='gray', linestyle=':', linewidth=1, alpha=0.7)
ax.text(0.305, len(assoc_df) - 0.5, 'Moderate (0.3)', fontsize=9, color='gray')

# 범례
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#3498db', edgecolor='black', label='Numeric (|Pearson r|)'),
    Patch(facecolor='#e74c3c', edgecolor='black', label="Categorical (Cramer's V)")
]
ax.legend(handles=legend_elements, loc='lower right', fontsize=10)

ax.set_title("Figure 7. Variable Association Strength with Accident Severity",
             fontsize=14, fontweight='bold')
ax.set_xlabel('Association Strength', fontsize=12)
ax.set_xlim(0, 0.7)

plt.tight_layout()
plt.show()

# ============================================================
# [시각화 8] 고위험 조건 조합 — Top 10 심각 사고 비율
# ============================================================
fig, ax = plt.subplots(figsize=(14, 7))

# Weather × Road_Condition × Traffic_Density 조합별 심각 비율 계산
td_map = {0: 'Low', 1: 'Medium', 2: 'High'}
df_temp = df.copy()
df_temp['TD_Label'] = df_temp['Traffic_Density'].map(td_map)

combo = df_temp.groupby(['Weather', 'Road_Condition', 'TD_Label']).agg(
    total=('Accident_Severity', 'count'),
    severe=('Accident_Severity', 'sum'),
    severe_rate=('Accident_Severity', 'mean')
).reset_index()

# 최소 10건 이상인 조합만 필터링 (신뢰성 확보)
combo = combo[combo['total'] >= 10]
combo['severe_rate_pct'] = (combo['severe_rate'] * 100).round(1)
combo['label'] = combo['Weather'] + ' + ' + combo['Road_Condition'] + ' + TD:' + combo['TD_Label']

# 심각 비율 상위 10개 조합
top10 = combo.nlargest(10, 'severe_rate_pct').sort_values('severe_rate_pct', ascending=True)

# 색상: 심각도에 따른 그라데이션
norm = plt.Normalize(vmin=0, vmax=100)
cmap = plt.cm.RdYlGn_r
bar_colors = [cmap(norm(val)) for val in top10['severe_rate_pct']]

bars = ax.barh(
    top10['label'], top10['severe_rate_pct'],
    color=bar_colors, edgecolor='black', linewidth=0.5, height=0.6
)

# 각 바 옆에 비율과 건수 표시
for bar, (_, row) in zip(bars, top10.iterrows()):
    ax.text(bar.get_width() + 1, bar.get_y() + bar.get_height()/2,
            f'{row["severe_rate_pct"]}% ({int(row["severe"])}/{int(row["total"])}건)',
            ha='left', va='center', fontsize=10, fontweight='bold')

ax.set_title('Figure 8. Top 10 High-Risk Condition Combinations (Severe Rate)',
             fontsize=14, fontweight='bold')
ax.set_xlabel('Severe Accident Rate (%)', fontsize=12)
ax.set_xlim(0, 110)

plt.tight_layout()
plt.show()

# ============================================================
# 분석 결과 텍스트 요약 출력
# ============================================================
print("=" * 70)
print("EDA 시각화 분석 결과 요약")
print("=" * 70)

print("""
[Figure 1] 타겟 변수 분포
  - Non-Severe(0) 71.5% vs Severe(1) 28.5%로 클래스 불균형 존재

[Figure 2] 핵심 범주형 변수
  - Weather: Rainy 66.8% > Stormy 44.4% > Clear 10.4%
  - Road_Condition: Wet 65.6% > Icy 35.6% > Dry 8.8%
  → 비/습한 조건에서 사고 심각도가 급격히 증가

[Figure 3] 핵심 수치형 변수
  - Traffic_Density: 밀도가 높을수록 심각 사고 비율 증가 (r=0.44)
  - Speed_Limit: 제한속도가 높을수록 심각 사고 소폭 증가 (r=0.25)

[Figure 4] Weather × Road_Condition 상호작용
  - Rainy+Wet 조합이 가장 높은 심각 비율 → 두 변수의 강한 연관성 확인

[Figure 5] Traffic_Density × Speed_Limit 상호작용
  - 교통 밀도 High + 제한속도 High 구간에서 심각 비율 극대화

[Figure 6] Driver_Age × Driver_Experience 다중공선성
  - r=0.9453의 매우 강한 선형 관계 → 모델링 시 하나 제거 필요

[Figure 7] 전체 변수 연관 강도 종합
  - Road_Condition(0.57) > Weather(0.55) > Traffic_Density(0.44) 순

[Figure 8] 고위험 조건 조합
  - 가장 위험한 조합: Rainy + Wet + 높은 교통 밀도
""")

# ============================================================
# [Step 18] 포괄적 EDA 시각화 세트 — 단일 PDF 저장
# ============================================================
# 전처리된 데이터셋의 주요 분포, 패턴, 관계, 특성을
# 8페이지 시각화로 구성하여 하나의 PDF 파일로 저장합니다.
#
#   Page 1: 타겟 변수(Accident_Severity) 분포
#   Page 2: 핵심 범주형 변수 — Weather & Road_Condition
#   Page 3: 핵심 수치형 변수 — Traffic_Density & Speed_Limit
#   Page 4: 전체 수치형 변수 분포 (히스토그램 + 박스플롯)
#   Page 5: Weather × Road_Condition 상호작용 히트맵
#   Page 6: Traffic_Density × Speed_Limit 상호작용
#   Page 7: 다중공선성 확인 & 상관계수 히트맵
#   Page 8: 변수 중요도 종합 & 고위험 조합 Top 10
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from matplotlib.backends.backend_pdf import PdfPages  # PDF 저장용
from matplotlib.patches import Patch

# ----------------------------------------------------------
# 0) 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

# 전처리 완료 데이터 불러오기
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'
df = pd.read_csv(data_path)

# 공통 색상 팔레트
SEV_COLORS = {0: '#3498db', 1: '#e74c3c'}    # 파랑: 비심각, 빨강: 심각
SEV_LABELS = {0: 'Non-Severe (0)', 1: 'Severe (1)'}
TD_LABELS = {0: 'Low (0)', 1: 'Medium (1)', 2: 'High (2)'}
TD_COLORS = {0: '#2ecc71', 1: '#f39c12', 2: '#e74c3c'}

# PDF 저장 경로
pdf_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/traffic_accident_visualization.pdf'

# ----------------------------------------------------------
# 공통 유틸 함수
# ----------------------------------------------------------
def add_page_header(fig, title, subtitle=''):
    """각 페이지 상단에 제목과 부제 추가"""
    fig.suptitle(title, fontsize=16, fontweight='bold', y=0.98)
    if subtitle:
        fig.text(0.5, 0.94, subtitle, ha='center', fontsize=10,
                 style='italic', color='gray')

# ============================================================
# PDF 생성 시작
# ============================================================
with PdfPages(pdf_path) as pdf:

    # ========================================================
    # PAGE 1: 타겟 변수 분포
    # ========================================================
    fig, axes = plt.subplots(1, 2, figsize=(14, 7))
    add_page_header(fig, 'Page 1. Target Variable Distribution',
                    'Accident_Severity: Binary Classification (0 = Non-Severe, 1 = Severe)')

    # (1-1) 바 차트
    counts = df['Accident_Severity'].value_counts().sort_index()
    bars = axes[0].bar(
        [SEV_LABELS[i] for i in counts.index], counts.values,
        color=[SEV_COLORS[i] for i in counts.index],
        edgecolor='black', linewidth=0.5, width=0.5
    )
    for bar, count in zip(bars, counts.values):
        pct = count / len(df) * 100
        axes[0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 10,
                     f'{count}\n({pct:.1f}%)', ha='center', va='bottom',
                     fontsize=11, fontweight='bold')
    axes[0].set_title('Frequency Distribution', fontsize=13, fontweight='bold')
    axes[0].set_ylabel('Count', fontsize=11)
    axes[0].set_ylim(0, max(counts.values) * 1.25)

    # (1-2) 파이 차트
    axes[1].pie(
        counts.values,
        labels=[SEV_LABELS[i] for i in counts.index],
        colors=[SEV_COLORS[i] for i in counts.index],
        autopct='%1.1f%%', startangle=90, explode=[0, 0.05],
        shadow=True, textprops={'fontsize': 12}
    )
    axes[1].set_title('Class Proportion', fontsize=13, fontweight='bold')

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 1 저장 완료")

    # ========================================================
    # PAGE 2: 핵심 범주형 변수 — Weather & Road_Condition
    # ========================================================
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    add_page_header(fig, 'Page 2. Key Categorical Variables vs Severity',
                    "Weather (Cramer's V=0.55) & Road_Condition (Cramer's V=0.57)")

    # (2-1) Weather 빈도
    w_order = df.groupby('Weather')['Accident_Severity'].mean().sort_values(ascending=False).index
    w_counts = df['Weather'].value_counts().loc[w_order]
    axes[0, 0].barh(w_counts.index, w_counts.values, color='#5dade2',
                     edgecolor='black', linewidth=0.5)
    for i, v in enumerate(w_counts.values):
        axes[0, 0].text(v + 5, i, str(v), va='center', fontsize=10)
    axes[0, 0].set_title('Weather — Frequency', fontsize=12, fontweight='bold')
    axes[0, 0].invert_yaxis()

    # (2-2) Weather별 심각 비율
    w_ct = pd.crosstab(df['Weather'], df['Accident_Severity'], normalize='index') * 100
    w_ct = w_ct.loc[w_order]
    w_ct.plot(kind='barh', stacked=True, ax=axes[0, 1],
              color=[SEV_COLORS[0], SEV_COLORS[1]], edgecolor='black', linewidth=0.5)
    for i, (idx, row) in enumerate(w_ct.iterrows()):
        axes[0, 1].text(row[0] + row[1]/2, i, f'{row[1]:.1f}%',
                         ha='center', va='center', fontsize=10,
                         fontweight='bold', color='white')
    axes[0, 1].set_title('Weather — Severe Rate (%)', fontsize=12, fontweight='bold')
    axes[0, 1].set_xlabel('Percentage (%)')
    axes[0, 1].set_ylabel('')
    axes[0, 1].legend(['Non-Severe', 'Severe'], loc='lower right', fontsize=9)
    axes[0, 1].invert_yaxis()

    # (2-3) Road_Condition 빈도
    r_order = df.groupby('Road_Condition')['Accident_Severity'].mean().sort_values(ascending=False).index
    r_counts = df['Road_Condition'].value_counts().loc[r_order]
    axes[1, 0].barh(r_counts.index, r_counts.values, color='#58d68d',
                     edgecolor='black', linewidth=0.5)
    for i, v in enumerate(r_counts.values):
        axes[1, 0].text(v + 5, i, str(v), va='center', fontsize=10)
    axes[1, 0].set_title('Road Condition — Frequency', fontsize=12, fontweight='bold')
    axes[1, 0].invert_yaxis()

    # (2-4) Road_Condition별 심각 비율
    r_ct = pd.crosstab(df['Road_Condition'], df['Accident_Severity'], normalize='index') * 100
    r_ct = r_ct.loc[r_order]
    r_ct.plot(kind='barh', stacked=True, ax=axes[1, 1],
              color=[SEV_COLORS[0], SEV_COLORS[1]], edgecolor='black', linewidth=0.5)
    for i, (idx, row) in enumerate(r_ct.iterrows()):
        axes[1, 1].text(row[0] + row[1]/2, i, f'{row[1]:.1f}%',
                         ha='center', va='center', fontsize=10,
                         fontweight='bold', color='white')
    axes[1, 1].set_title('Road Condition — Severe Rate (%)', fontsize=12, fontweight='bold')
    axes[1, 1].set_xlabel('Percentage (%)')
    axes[1, 1].set_ylabel('')
    axes[1, 1].legend(['Non-Severe', 'Severe'], loc='lower right', fontsize=9)
    axes[1, 1].invert_yaxis()

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 2 저장 완료")

    # ========================================================
    # PAGE 3: 핵심 수치형 변수 — Traffic_Density & Speed_Limit
    # ========================================================
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    add_page_header(fig, 'Page 3. Key Numeric Variables vs Severity',
                    'Traffic_Density (r=0.44) & Speed_Limit (r=0.25)')

    # (3-1) Traffic_Density별 심각도 그룹 바
    td_ct = pd.crosstab(df['Traffic_Density'], df['Accident_Severity'])
    x_pos = range(len(td_ct.index))
    width = 0.35
    axes[0, 0].bar([x - width/2 for x in x_pos], td_ct[0], width,
                    label='Non-Severe', color=SEV_COLORS[0], edgecolor='black', linewidth=0.5)
    axes[0, 0].bar([x + width/2 for x in x_pos], td_ct[1], width,
                    label='Severe', color=SEV_COLORS[1], edgecolor='black', linewidth=0.5)
    axes[0, 0].set_xticks(x_pos)
    axes[0, 0].set_xticklabels([TD_LABELS[i] for i in td_ct.index])
    axes[0, 0].set_title('Traffic Density — Count by Severity', fontsize=12, fontweight='bold')
    axes[0, 0].set_ylabel('Count')
    axes[0, 0].legend(fontsize=9)

    # (3-2) Traffic_Density별 심각 비율
    td_pct = pd.crosstab(df['Traffic_Density'], df['Accident_Severity'], normalize='index') * 100
    bars = axes[0, 1].bar(
        [TD_LABELS[i] for i in td_pct.index], td_pct[1],
        color=[TD_COLORS[i] for i in td_pct.index],
        edgecolor='black', linewidth=0.5, width=0.5
    )
    for bar, val in zip(bars, td_pct[1]):
        axes[0, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                         f'{val:.1f}%', ha='center', fontsize=11, fontweight='bold')
    axes[0, 1].set_title('Traffic Density — Severe Rate (%)', fontsize=12, fontweight='bold')
    axes[0, 1].set_ylabel('Severe Rate (%)')
    axes[0, 1].set_ylim(0, td_pct[1].max() * 1.3)

    # (3-3) Speed_Limit 히스토그램 (심각도별)
    for sev in [0, 1]:
        subset = df[df['Accident_Severity'] == sev]['Speed_Limit']
        axes[1, 0].hist(subset, bins=15, alpha=0.6, label=SEV_LABELS[sev],
                         color=SEV_COLORS[sev], edgecolor='black', linewidth=0.5)
        axes[1, 0].axvline(subset.mean(), color=SEV_COLORS[sev],
                            linestyle='--', linewidth=2, alpha=0.8)
    axes[1, 0].set_title('Speed Limit — Distribution by Severity', fontsize=12, fontweight='bold')
    axes[1, 0].set_xlabel('Speed Limit')
    axes[1, 0].set_ylabel('Count')
    axes[1, 0].legend(fontsize=9)

    # (3-4) Speed_Limit 박스플롯 (심각도별)
    df_plot_temp = df.copy()
    df_plot_temp['Severity_Label'] = df_plot_temp['Accident_Severity'].map(SEV_LABELS)
    sns.boxplot(data=df_plot_temp, x='Severity_Label', y='Speed_Limit',
                palette=[SEV_COLORS[0], SEV_COLORS[1]], ax=axes[1, 1])
    # 각 그룹의 평균 표시
    for sev in [0, 1]:
        mean_val = df[df['Accident_Severity'] == sev]['Speed_Limit'].mean()
        axes[1, 1].text(sev, mean_val + 2, f'Mean: {mean_val:.1f}',
                         ha='center', fontsize=10, fontweight='bold', color=SEV_COLORS[sev])
    axes[1, 1].set_title('Speed Limit — Box Plot by Severity', fontsize=12, fontweight='bold')
    axes[1, 1].set_xlabel('')
    axes[1, 1].set_ylabel('Speed Limit')

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 3 저장 완료")

    # ========================================================
    # PAGE 4: 전체 수치형 변수 분포 (히스토그램 + 박스플롯)
    # ========================================================
    num_cols = ['Traffic_Density', 'Speed_Limit', 'Driver_Age',
                'Driver_Experience', 'Accident_Severity']

    fig, axes = plt.subplots(2, 5, figsize=(18, 8))
    add_page_header(fig, 'Page 4. All Numeric Variables — Distribution Overview',
                    'Top: Histogram | Bottom: Box Plot (by Severity)')

    for i, col in enumerate(num_cols):
        # 상단: 히스토그램
        axes[0, i].hist(df[col], bins=20, color='#5dade2', edgecolor='black',
                         linewidth=0.5, alpha=0.8)
        axes[0, i].axvline(df[col].mean(), color='red', linestyle='--',
                            linewidth=1.5, label=f'Mean: {df[col].mean():.1f}')
        axes[0, i].axvline(df[col].median(), color='green', linestyle='-.',
                            linewidth=1.5, label=f'Med: {df[col].median():.1f}')
        axes[0, i].set_title(col, fontsize=10, fontweight='bold')
        axes[0, i].legend(fontsize=7)
        axes[0, i].tick_params(labelsize=8)

        # 하단: 박스플롯 (심각도별)
        data_0 = df[df['Accident_Severity'] == 0][col]
        data_1 = df[df['Accident_Severity'] == 1][col]
        bp = axes[1, i].boxplot([data_0, data_1], labels=['Sev=0', 'Sev=1'],
                                 patch_artist=True, widths=0.5)
        bp['boxes'][0].set_facecolor(SEV_COLORS[0])
        bp['boxes'][1].set_facecolor(SEV_COLORS[1])
        for box in bp['boxes']:
            box.set_alpha(0.6)
        axes[1, i].set_title(f'{col}\nby Severity', fontsize=9, fontweight='bold')
        axes[1, i].tick_params(labelsize=8)

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 4 저장 완료")

    # ========================================================
    # PAGE 5: Weather × Road_Condition 상호작용
    # ========================================================
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))
    add_page_header(fig, 'Page 5. Interaction: Weather × Road Condition',
                    'Frequency (left) & Severe Accident Rate (right)')

    weather_order = ['Rainy', 'Stormy', 'Snowy', 'Foggy', 'Clear']
    road_order = ['Wet', 'Icy', 'Under Construction', 'Dry']

    # (5-1) 교차 빈도 히트맵
    cross_count = pd.crosstab(df['Weather'], df['Road_Condition'])
    cross_count = cross_count.reindex(index=weather_order, columns=road_order)
    sns.heatmap(cross_count, annot=True, fmt='d', cmap='YlOrRd',
                linewidths=0.5, ax=axes[0], cbar_kws={'label': 'Count'})
    axes[0].set_title('Cross-tabulation (Frequency)', fontsize=12, fontweight='bold')

    # (5-2) 교차 심각 비율 히트맵
    cross_sev = df.groupby(['Weather', 'Road_Condition'])['Accident_Severity'].mean() * 100
    cross_sev = cross_sev.unstack().reindex(index=weather_order, columns=road_order)
    sns.heatmap(cross_sev, annot=True, fmt='.1f', cmap='RdYlGn_r',
                linewidths=0.5, ax=axes[1], vmin=0, vmax=100,
                cbar_kws={'label': 'Severe Rate (%)'})
    axes[1].set_title('Severe Accident Rate (%)', fontsize=12, fontweight='bold')

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 5 저장 완료")

    # ========================================================
    # PAGE 6: Traffic_Density × Speed_Limit 상호작용
    # ========================================================
    fig, axes = plt.subplots(1, 2, figsize=(14, 7))
    add_page_header(fig, 'Page 6. Interaction: Traffic Density × Speed Limit',
                    'Box Plot by Severity (left) & Severe Rate by Speed Group (right)')

    # (6-1) Traffic_Density별 Speed_Limit 박스플롯 (심각도 구분)
    df_p = df.copy()
    df_p['TD_Label'] = df_p['Traffic_Density'].map(TD_LABELS)
    df_p['Sev_Label'] = df_p['Accident_Severity'].map(SEV_LABELS)
    sns.boxplot(data=df_p, x='TD_Label', y='Speed_Limit', hue='Sev_Label',
                palette=[SEV_COLORS[0], SEV_COLORS[1]], ax=axes[0],
                order=['Low (0)', 'Medium (1)', 'High (2)'])
    axes[0].set_title('Speed Limit by Density & Severity', fontsize=12, fontweight='bold')
    axes[0].set_xlabel('Traffic Density')
    axes[0].legend(title='Severity', fontsize=9)

    # (6-2) Speed 구간 × Traffic_Density별 심각 비율
    df_p['Speed_Group'] = pd.cut(df_p['Speed_Limit'], bins=[0, 50, 80, 110],
                                  labels=['Low (<=50)', 'Mid (51-80)', 'High (81-110)'])
    interact_rate = df_p.groupby(['Speed_Group', 'Traffic_Density'])['Accident_Severity'].mean() * 100
    interact_pivot = interact_rate.unstack()
    interact_pivot.columns = [TD_LABELS[c] for c in interact_pivot.columns]
    interact_pivot.plot(kind='bar', ax=axes[1],
                        color=[TD_COLORS[0], TD_COLORS[1], TD_COLORS[2]],
                        edgecolor='black', linewidth=0.5)
    for container in axes[1].containers:
        axes[1].bar_label(container, fmt='%.1f%%', fontsize=8, padding=2)
    axes[1].set_title('Severe Rate: Speed Group × Density', fontsize=12, fontweight='bold')
    axes[1].set_xlabel('Speed Limit Group')
    axes[1].set_ylabel('Severe Rate (%)')
    axes[1].set_xticklabels(axes[1].get_xticklabels(), rotation=0)
    axes[1].legend(title='Traffic Density', fontsize=9)
    axes[1].set_ylim(0, 100)

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 6 저장 완료")

    # ========================================================
    # PAGE 7: 다중공선성 & 상관계수 히트맵
    # ========================================================
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))
    add_page_header(fig, 'Page 7. Correlation Analysis',
                    'Multicollinearity Check (left) & Pearson Correlation Heatmap (right)')

    # (7-1) Driver_Age vs Driver_Experience 산점도
    for sev in [0, 1]:
        subset = df[df['Accident_Severity'] == sev]
        axes[0].scatter(subset['Driver_Age'], subset['Driver_Experience'],
                         c=SEV_COLORS[sev], label=SEV_LABELS[sev],
                         alpha=0.35, s=25, edgecolors='white', linewidth=0.3)
    # 회귀선
    z = np.polyfit(df['Driver_Age'], df['Driver_Experience'], 1)
    p_line = np.poly1d(z)
    age_range = np.linspace(df['Driver_Age'].min(), df['Driver_Age'].max(), 100)
    r_val = df['Driver_Age'].corr(df['Driver_Experience'])
    axes[0].plot(age_range, p_line(age_range), '--', color='black', linewidth=2,
                  alpha=0.7, label=f'r = {r_val:.4f}')
    axes[0].set_title('Driver Age vs Experience\n(Multicollinearity)', fontsize=12, fontweight='bold')
    axes[0].set_xlabel('Driver Age')
    axes[0].set_ylabel('Driver Experience (years)')
    axes[0].legend(fontsize=9)

    # (7-2) 피어슨 상관계수 히트맵
    corr_cols = ['Traffic_Density', 'Speed_Limit', 'Driver_Age',
                 'Driver_Experience', 'Accident_Severity']
    corr_matrix = df[corr_cols].corr().round(4)
    sns.heatmap(corr_matrix, annot=True, fmt='.4f', cmap='RdBu_r',
                center=0, vmin=-1, vmax=1, square=True, linewidths=0.5,
                ax=axes[1], cbar_kws={'shrink': 0.8, 'label': 'Pearson r'})
    axes[1].set_title('Pearson Correlation Heatmap', fontsize=12, fontweight='bold')
    axes[1].set_xticklabels(axes[1].get_xticklabels(), rotation=45, ha='right')

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 7 저장 완료")

    # ========================================================
    # PAGE 8: 변수 중요도 종합 & 고위험 조합 Top 10
    # ========================================================
    fig, axes = plt.subplots(1, 2, figsize=(16, 8))
    add_page_header(fig, 'Page 8. Variable Importance & High-Risk Combinations',
                    "Association with Severity (left) & Top 10 Risk Combos (right)")

    # (8-1) 전체 변수 연관 강도
    # 수치형: Pearson r 절대값
    num_assoc = {
        'Traffic_Density': abs(df['Traffic_Density'].corr(df['Accident_Severity'])),
        'Speed_Limit': abs(df['Speed_Limit'].corr(df['Accident_Severity'])),
        'Driver_Age': abs(df['Driver_Age'].corr(df['Accident_Severity'])),
        'Driver_Experience': abs(df['Driver_Experience'].corr(df['Accident_Severity']))
    }
    # 범주형: Cramer's V
    cat_cols = ['Weather', 'Road_Type', 'Time_of_Day',
                'Road_Condition', 'Vehicle_Type', 'Road_Light_Condition']
    cat_assoc = {}
    for col in cat_cols:
        ct = pd.crosstab(df[col], df['Accident_Severity'])
        chi2, _, _, _ = stats.chi2_contingency(ct)
        n = ct.sum().sum()
        min_d = min(ct.shape) - 1
        cat_assoc[col] = np.sqrt(chi2 / (n * min_d))

    all_assoc = {**num_assoc, **cat_assoc}
    assoc_df = pd.DataFrame({
        'Variable': list(all_assoc.keys()),
        'Association': list(all_assoc.values()),
        'Type': (['Numeric'] * 4 + ['Categorical'] * 6)
    }).sort_values('Association', ascending=True)

    bar_colors = ['#e74c3c' if t == 'Categorical' else '#3498db' for t in assoc_df['Type']]
    bars = axes[0].barh(assoc_df['Variable'], assoc_df['Association'],
                         color=bar_colors, edgecolor='black', linewidth=0.5, height=0.6)
    for bar, val in zip(bars, assoc_df['Association']):
        axes[0].text(bar.get_width() + 0.008, bar.get_y() + bar.get_height()/2,
                     f'{val:.4f}', ha='left', va='center', fontsize=9, fontweight='bold')
    axes[0].axvline(0.1, color='gray', linestyle=':', linewidth=1, alpha=0.6)
    axes[0].axvline(0.3, color='gray', linestyle=':', linewidth=1, alpha=0.6)
    axes[0].text(0.105, len(assoc_df) - 0.5, '0.1', fontsize=8, color='gray')
    axes[0].text(0.305, len(assoc_df) - 0.5, '0.3', fontsize=8, color='gray')
    legend_el = [Patch(facecolor='#3498db', edgecolor='black', label='Numeric (|r|)'),
                 Patch(facecolor='#e74c3c', edgecolor='black', label="Categorical (Cramer's V)")]
    axes[0].legend(handles=legend_el, loc='lower right', fontsize=9)
    axes[0].set_title('Variable Association Strength', fontsize=12, fontweight='bold')
    axes[0].set_xlabel('Association Strength')
    axes[0].set_xlim(0, 0.7)

    # (8-2) 고위험 조합 Top 10
    td_map = {0: 'Low', 1: 'Medium', 2: 'High'}
    df_t = df.copy()
    df_t['TD_Label'] = df_t['Traffic_Density'].map(td_map)
    combo = df_t.groupby(['Weather', 'Road_Condition', 'TD_Label']).agg(
        total=('Accident_Severity', 'count'),
        severe=('Accident_Severity', 'sum'),
        severe_rate=('Accident_Severity', 'mean')
    ).reset_index()
    combo = combo[combo['total'] >= 10]  # 최소 10건 이상 필터링
    combo['pct'] = (combo['severe_rate'] * 100).round(1)
    combo['label'] = combo['Weather'] + ' + ' + combo['Road_Condition'] + '\n(TD: ' + combo['TD_Label'] + ')'
    top10 = combo.nlargest(10, 'pct').sort_values('pct', ascending=True)

    norm = plt.Normalize(vmin=0, vmax=100)
    cmap = plt.cm.RdYlGn_r
    t_bars = axes[1].barh(top10['label'], top10['pct'],
                           color=[cmap(norm(v)) for v in top10['pct']],
                           edgecolor='black', linewidth=0.5, height=0.6)
    for bar, (_, row) in zip(t_bars, top10.iterrows()):
        axes[1].text(bar.get_width() + 1, bar.get_y() + bar.get_height()/2,
                     f'{row["pct"]}% ({int(row["severe"])}/{int(row["total"])})',
                     ha='left', va='center', fontsize=9, fontweight='bold')
    axes[1].set_title('Top 10 High-Risk Combinations', fontsize=12, fontweight='bold')
    axes[1].set_xlabel('Severe Rate (%)')
    axes[1].set_xlim(0, 110)

    plt.tight_layout(rect=[0, 0, 1, 0.92])
    pdf.savefig(fig)
    plt.close()
    print("✔ Page 8 저장 완료")

# ============================================================
# 완료 메시지
# ============================================================
print("\n" + "=" * 70)
print("PDF 저장 완료!")
print("=" * 70)
print(f"파일 경로: {pdf_path}")
print(f"총 페이지: 8페이지")
print("=" * 70)

# ============================================================
# [시각화 2] 교통 밀도의 증폭 효과 — 패싯 히트맵
# ============================================================
# 분석 질문: 교통 혼잡도(Traffic_Density)는 다른 위험 요인의
#            효과를 얼마나 증폭시키는가?
# 관련 변수: Traffic_Density, Weather, Road_Condition, Accident_Severity
# 시각화 유형: 패싯 히트맵 (Traffic_Density 수준별 3개 히트맵)
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

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

# ----------------------------------------------------------
# 0) 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 120

# 전처리 완료 데이터 불러오기
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'
df = pd.read_csv(data_path)

# Traffic_Density 라벨 매핑
td_map = {0: 'Low (0)', 1: 'Medium (1)', 2: 'High (2)'}
td_levels = [0, 1, 2]

# 정렬 순서 정의 (심각 비율 높은 순)
weather_order = ['Rainy', 'Stormy', 'Snowy', 'Foggy', 'Clear']
road_order = ['Wet', 'Icy', 'Under Construction', 'Dry']

# ----------------------------------------------------------
# 1) 패싯 히트맵: Traffic_Density 수준별 심각 비율
# ----------------------------------------------------------
fig, axes = plt.subplots(1, 3, figsize=(20, 6))

fig.suptitle('Traffic Density Amplification Effect\nSevere Accident Rate (%) by Weather × Road Condition at Each Density Level',
             fontsize=14, fontweight='bold', y=1.05)

# 전체 데이터의 심각 비율 (비교 기준)
overall_rate = df['Accident_Severity'].mean() * 100

for idx, td in enumerate(td_levels):
    # 해당 Traffic_Density 수준의 데이터만 필터링
    df_sub = df[df['Traffic_Density'] == td]

    # Weather × Road_Condition 교차 심각 비율 계산
    cross_rate = df_sub.groupby(['Weather', 'Road_Condition'])['Accident_Severity'].mean() * 100
    cross_rate = cross_rate.unstack()

    # 순서 재정렬 (존재하는 컬럼/인덱스만)
    avail_weather = [w for w in weather_order if w in cross_rate.index]
    avail_road = [r for r in road_order if r in cross_rate.columns]
    cross_rate = cross_rate.reindex(index=avail_weather, columns=avail_road)

    # 건수 행렬도 계산 (셀 내 표기용)
    cross_count = df_sub.groupby(['Weather', 'Road_Condition'])['Accident_Severity'].count()
    cross_count = cross_count.unstack().reindex(index=avail_weather, columns=avail_road)

    # annot 텍스트: 비율% + (건수) 형태
    annot_text = cross_rate.copy().astype(str)
    for w in avail_weather:
        for r in avail_road:
            rate_val = cross_rate.loc[w, r] if not pd.isna(cross_rate.loc[w, r]) else np.nan
            count_val = cross_count.loc[w, r] if not pd.isna(cross_count.loc[w, r]) else 0
            if pd.isna(rate_val) or count_val == 0:
                annot_text.loc[w, r] = '-'
            else:
                annot_text.loc[w, r] = f'{rate_val:.0f}%\n(n={int(count_val)})'

    # 히트맵 그리기
    sns.heatmap(
        cross_rate, annot=annot_text, fmt='',
        cmap='RdYlGn_r', vmin=0, vmax=100,
        linewidths=0.8, linecolor='white',
        ax=axes[idx],
        cbar=True if idx == 2 else False,  # 마지막 축에만 컬러바
        cbar_kws={'label': 'Severe Rate (%)', 'shrink': 0.8} if idx == 2 else {},
        annot_kws={'fontsize': 9}
    )

    # 해당 밀도 수준의 전체 심각 비율
    sub_rate = df_sub['Accident_Severity'].mean() * 100
    sub_count = len(df_sub)

    axes[idx].set_title(
        f'Traffic Density: {td_map[td]}\n(n={sub_count}, Overall Severe Rate: {sub_rate:.1f}%)',
        fontsize=11, fontweight='bold', pad=10
    )
    axes[idx].set_xlabel('Road Condition', fontsize=10)
    axes[idx].set_ylabel('Weather' if idx == 0 else '', fontsize=10)

    # Y축 라벨 첫 번째만 표시
    if idx > 0:
        axes[idx].set_yticklabels([])

plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 2) 증폭 효과 요약표 출력
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("교통 밀도별 심각 사고 비율 증폭 효과 요약")
print("=" * 70)

# 주요 조합별 밀도에 따른 심각 비율 변화
key_combos = [
    ('Rainy', 'Wet'),
    ('Clear', 'Dry'),
    ('Snowy', 'Icy'),
    ('Clear', 'Wet')
]

print(f"\n{'날씨':<10} {'도로상태':<20} {'Low(0)':<12} {'Medium(1)':<12} {'High(2)':<12} {'증폭배수':<10}")
print("-" * 76)

for weather, road in key_combos:
    rates = []
    for td in td_levels:
        sub = df[(df['Traffic_Density'] == td) &
                 (df['Weather'] == weather) &
                 (df['Road_Condition'] == road)]
        if len(sub) >= 3:  # 최소 3건 이상
            rate = sub['Accident_Severity'].mean() * 100
            rates.append(f'{rate:.1f}%')
        else:
            rates.append('N/A')

    # 증폭 배수 계산 (Low 대비 High)
    sub_low = df[(df['Traffic_Density'] == 0) & (df['Weather'] == weather) & (df['Road_Condition'] == road)]
    sub_high = df[(df['Traffic_Density'] == 2) & (df['Weather'] == weather) & (df['Road_Condition'] == road)]
    if len(sub_low) >= 3 and len(sub_high) >= 3:
        r_low = sub_low['Accident_Severity'].mean()
        r_high = sub_high['Accident_Severity'].mean()
        if r_low > 0:
            amplify = f'{r_high / r_low:.1f}x'
        else:
            amplify = '-'
    else:
        amplify = 'N/A'

    print(f'{weather:<10} {road:<20} {rates[0]:<12} {rates[1]:<12} {rates[2]:<12} {amplify:<10}')

print("\n※ 증폭배수 = High(2) 심각비율 / Low(0) 심각비율")
print("※ N/A = 해당 조합의 데이터가 3건 미만으로 신뢰도 부족")

# ============================================================
# [시각화 5] 차량 유형별 숨겨진 패턴 — 조건부 스택 바 차트
# ============================================================
# 분석 질문: 차량 유형(Car/Truck/Bus)은 전체적으로는 사고 심각도와
#            유의하지 않지만, 혼잡한 고속도로에서는 대형 차량이
#            더 위험한가?
# 관련 변수: Vehicle_Type, Road_Type, Traffic_Density, Accident_Severity
# 시각화 유형: 조건별 그룹 바 차트 + 전체 vs 특정 조건 비교
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

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

# ----------------------------------------------------------
# 0) 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 120

# 전처리 완료 데이터 불러오기
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'
df = pd.read_csv(data_path)

# 색상 정의
vehicle_colors = {'Car': '#3498db', 'Truck': '#e67e22', 'Bus': '#e74c3c'}
vehicle_order = ['Car', 'Truck', 'Bus']

# ----------------------------------------------------------
# 1) 4가지 조건에서 차량 유형별 심각 비율 비교
# ----------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(16, 11))

fig.suptitle('Vehicle Type Hidden Patterns — Severe Rate by Conditions',
             fontsize=15, fontweight='bold', y=1.02)

# 분석할 4가지 조건 정의
conditions = [
    {
        'title': 'Overall (All Data)',
        'subtitle': 'n=840',
        'filter': df,
        'ax': axes[0, 0]
    },
    {
        'title': 'Highway Only',
        'subtitle': '',
        'filter': df[df['Road_Type'] == 'Highway'],
        'ax': axes[0, 1]
    },
    {
        'title': 'Highway + High Traffic Density',
        'subtitle': '',
        'filter': df[(df['Road_Type'] == 'Highway') & (df['Traffic_Density'] == 2)],
        'ax': axes[1, 0]
    },
    {
        'title': 'Highway + High Density + Wet/Rainy',
        'subtitle': '',
        'filter': df[(df['Road_Type'] == 'Highway') &
                     (df['Traffic_Density'] == 2) &
                     ((df['Road_Condition'] == 'Wet') | (df['Weather'] == 'Rainy'))],
        'ax': axes[1, 1]
    }
]

for cond in conditions:
    ax = cond['ax']
    df_sub = cond['filter']

    # 차량 유형별 심각 비율 및 건수 계산
    results = []
    for vtype in vehicle_order:
        sub = df_sub[df_sub['Vehicle_Type'] == vtype]
        total = len(sub)
        if total > 0:
            severe = sub['Accident_Severity'].sum()
            rate = severe / total * 100
        else:
            severe = 0
            rate = 0
        results.append({
            'type': vtype,
            'total': total,
            'severe': int(severe),
            'rate': rate
        })

    # 막대 그래프 그리기
    x_pos = range(len(results))
    bars = ax.bar(
        [r['type'] for r in results],
        [r['rate'] for r in results],
        color=[vehicle_colors[r['type']] for r in results],
        edgecolor='black', linewidth=0.5, width=0.55
    )

    # 각 바 위에 비율과 건수 표시
    for bar, r in zip(bars, results):
        if r['total'] > 0:
            ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1.5,
                    f'{r["rate"]:.1f}%\n({r["severe"]}/{r["total"]})',
                    ha='center', va='bottom', fontsize=10, fontweight='bold')

    # 해당 조건의 전체 심각 비율 기준선
    if len(df_sub) > 0:
        overall = df_sub['Accident_Severity'].mean() * 100
        ax.axhline(y=overall, color='gray', linestyle='--', linewidth=1.5, alpha=0.7)
        ax.text(len(results) - 0.5, overall + 1, f'Avg: {overall:.1f}%',
                fontsize=9, color='gray', fontweight='bold', ha='right')

    # 제목 설정
    sub_n = len(df_sub)
    ax.set_title(f'{cond["title"]}\n(n={sub_n})', fontsize=12, fontweight='bold')
    ax.set_ylabel('Severe Accident Rate (%)', fontsize=10)
    ax.set_ylim(0, min(max([r['rate'] for r in results]) * 1.4, 105))

plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 2) Road_Type × Vehicle_Type × Traffic_Density 종합 히트맵
# ----------------------------------------------------------
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

fig.suptitle('Severe Rate (%) by Vehicle Type × Road Type at Each Traffic Density',
             fontsize=14, fontweight='bold', y=1.05)

td_map = {0: 'Low (0)', 1: 'Medium (1)', 2: 'High (2)'}
road_type_order = ['Highway', 'City Road', 'Rural Road', 'Mountain Road']

for idx, td in enumerate([0, 1, 2]):
    df_sub = df[df['Traffic_Density'] == td]

    # 교차 심각 비율 계산
    cross = df_sub.groupby(['Vehicle_Type', 'Road_Type']).agg(
        rate=('Accident_Severity', 'mean'),
        count=('Accident_Severity', 'count')
    ).reset_index()

    # 피벗 테이블
    pivot_rate = cross.pivot(index='Vehicle_Type', columns='Road_Type', values='rate') * 100
    pivot_count = cross.pivot(index='Vehicle_Type', columns='Road_Type', values='count')

    # 순서 정렬
    avail_vtype = [v for v in vehicle_order if v in pivot_rate.index]
    avail_road = [r for r in road_type_order if r in pivot_rate.columns]
    pivot_rate = pivot_rate.reindex(index=avail_vtype, columns=avail_road)
    pivot_count = pivot_count.reindex(index=avail_vtype, columns=avail_road)

    # annot 텍스트
    annot = pivot_rate.copy().astype(str)
    for v in avail_vtype:
        for r in avail_road:
            rv = pivot_rate.loc[v, r] if not pd.isna(pivot_rate.loc[v, r]) else np.nan
            cv = pivot_count.loc[v, r] if not pd.isna(pivot_count.loc[v, r]) else 0
            if pd.isna(rv) or cv == 0:
                annot.loc[v, r] = '-'
            else:
                annot.loc[v, r] = f'{rv:.0f}%\n(n={int(cv)})'

    sns.heatmap(
        pivot_rate, annot=annot, fmt='',
        cmap='YlOrRd', vmin=0, vmax=100,
        linewidths=0.8, linecolor='white',
        ax=axes[idx],
        cbar=True if idx == 2 else False,
        cbar_kws={'label': 'Severe Rate (%)', 'shrink': 0.8} if idx == 2 else {},
        annot_kws={'fontsize': 9}
    )

    sub_rate = df_sub['Accident_Severity'].mean() * 100
    axes[idx].set_title(f'Density: {td_map[td]}\n(Avg Severe: {sub_rate:.1f}%)',
                         fontsize=11, fontweight='bold')
    axes[idx].set_xlabel('Road Type', fontsize=10)
    axes[idx].set_ylabel('Vehicle Type' if idx == 0 else '', fontsize=10)
    if idx > 0:
        axes[idx].set_yticklabels([])

plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 3) 결과 요약 출력
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("차량 유형별 심각 사고 비율 — 조건별 비교")
print("=" * 70)
print(f"\n{'조건':<40} {'Car':<15} {'Truck':<15} {'Bus':<15}")
print("-" * 85)

for cond in conditions:
    df_sub = cond['filter']
    rates = []
    for vtype in vehicle_order:
        sub = df_sub[df_sub['Vehicle_Type'] == vtype]
        if len(sub) >= 5:
            rate = sub['Accident_Severity'].mean() * 100
            rates.append(f'{rate:.1f}% (n={len(sub)})')
        else:
            rates.append(f'N/A (n={len(sub)})')
    print(f'{cond["title"]:<40} {rates[0]:<15} {rates[1]:<15} {rates[2]:<15}')

# ============================================================
# [시각화 7] 복합 위험 점수 시각화 — 위험 요인 중첩 효과
# ============================================================
# 분석 질문: 위험 요인이 몇 개 겹칠 때 심각 사고 비율이
#            급등하는 임계점이 존재하는가?
# 관련 변수: Weather, Road_Condition, Traffic_Density, Speed_Limit,
#            Accident_Severity
# 시각화 유형: 계단형 바 차트 + 누적 구성 분석
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

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

# ----------------------------------------------------------
# 0) 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 120

# 전처리 완료 데이터 불러오기
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'
df = pd.read_csv(data_path)

# ----------------------------------------------------------
# 1) 위험 요인 플래그 생성 (각 요인별 0 또는 1)
# ----------------------------------------------------------
# 위험 요인 1: 나쁜 날씨 (Clear가 아닌 모든 날씨)
df['Risk_Weather'] = (df['Weather'] != 'Clear').astype(int)

# 위험 요인 2: 위험한 도로 상태 (Dry가 아닌 모든 상태)
df['Risk_Road'] = (df['Road_Condition'] != 'Dry').astype(int)

# 위험 요인 3: 높은 교통 혼잡도 (Traffic_Density == 2)
df['Risk_Density'] = (df['Traffic_Density'] == 2).astype(int)

# 위험 요인 4: 높은 제한 속도 (80 초과)
df['Risk_Speed'] = (df['Speed_Limit'] > 80).astype(int)

# 복합 위험 점수 (0~4): 위험 요인 개수의 합
df['Risk_Score'] = (df['Risk_Weather'] + df['Risk_Road'] +
                    df['Risk_Density'] + df['Risk_Speed'])

# ----------------------------------------------------------
# 2) 위험 점수별 통계 계산
# ----------------------------------------------------------
risk_stats = df.groupby('Risk_Score').agg(
    total=('Accident_Severity', 'count'),      # 전체 건수
    severe=('Accident_Severity', 'sum'),        # 심각 사고 건수
    rate=('Accident_Severity', 'mean')          # 심각 사고 비율
).reset_index()
risk_stats['rate_pct'] = (risk_stats['rate'] * 100).round(1)
risk_stats['pct_of_total'] = (risk_stats['total'] / len(df) * 100).round(1)

# ----------------------------------------------------------
# 3) 시각화: 2×2 대시보드
# ----------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(16, 12))

fig.suptitle('Compound Risk Score Analysis\nHow Many Risk Factors Overlap → Severe Accident Rate',
             fontsize=15, fontweight='bold', y=1.02)

# --- (3-1) 메인: 위험 점수별 심각 사고 비율 (계단형 바 차트) ---
colors_gradient = ['#27ae60', '#f1c40f', '#e67e22', '#e74c3c', '#8e44ad']

bars = axes[0, 0].bar(
    risk_stats['Risk_Score'], risk_stats['rate_pct'],
    color=[colors_gradient[i] for i in risk_stats['Risk_Score']],
    edgecolor='black', linewidth=0.8, width=0.6
)

# 각 바 위에 비율과 건수 표시
for bar, (_, row) in zip(bars, risk_stats.iterrows()):
    axes[0, 0].text(
        bar.get_x() + bar.get_width()/2, bar.get_height() + 2,
        f'{row["rate_pct"]}%\n({int(row["severe"])}/{int(row["total"])}건)',
        ha='center', va='bottom', fontsize=10, fontweight='bold'
    )

# 전체 평균선
avg_rate = df['Accident_Severity'].mean() * 100
axes[0, 0].axhline(y=avg_rate, color='gray', linestyle='--', linewidth=1.5, alpha=0.7)
axes[0, 0].text(4.3, avg_rate + 1, f'Overall Avg: {avg_rate:.1f}%',
                 fontsize=9, color='gray', fontweight='bold')

# 임계점 강조 (50% 기준선)
axes[0, 0].axhline(y=50, color='red', linestyle=':', linewidth=1, alpha=0.5)
axes[0, 0].text(4.3, 51, '50% Threshold', fontsize=8, color='red', alpha=0.7)

axes[0, 0].set_title('Severe Rate by Risk Score', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('Number of Risk Factors (0~4)', fontsize=10)
axes[0, 0].set_ylabel('Severe Accident Rate (%)', fontsize=10)
axes[0, 0].set_xticks(range(5))
axes[0, 0].set_xticklabels(['0\n(Safe)', '1', '2', '3', '4\n(Most\nDangerous)'])
axes[0, 0].set_ylim(0, max(risk_stats['rate_pct']) * 1.25)

# --- (3-2) 건수 분포: 위험 점수별 데이터 비중 ---
bars2 = axes[0, 1].bar(
    risk_stats['Risk_Score'], risk_stats['total'],
    color=[colors_gradient[i] for i in risk_stats['Risk_Score']],
    edgecolor='black', linewidth=0.8, width=0.6, alpha=0.8
)
for bar, (_, row) in zip(bars2, risk_stats.iterrows()):
    axes[0, 1].text(
        bar.get_x() + bar.get_width()/2, bar.get_height() + 5,
        f'{int(row["total"])}건\n({row["pct_of_total"]}%)',
        ha='center', va='bottom', fontsize=10, fontweight='bold'
    )
axes[0, 1].set_title('Data Distribution by Risk Score', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('Number of Risk Factors (0~4)', fontsize=10)
axes[0, 1].set_ylabel('Number of Accidents', fontsize=10)
axes[0, 1].set_xticks(range(5))
axes[0, 1].set_ylim(0, max(risk_stats['total']) * 1.25)

# --- (3-3) 개별 위험 요인의 기여도 ---
risk_factors = {
    'Bad Weather\n(Not Clear)': 'Risk_Weather',
    'Bad Road\n(Not Dry)': 'Risk_Road',
    'High Density\n(TD=2)': 'Risk_Density',
    'High Speed\n(>80 km/h)': 'Risk_Speed'
}

factor_rates = []
for label, col in risk_factors.items():
    exposed = df[df[col] == 1]['Accident_Severity'].mean() * 100
    not_exposed = df[df[col] == 0]['Accident_Severity'].mean() * 100
    factor_rates.append({
        'factor': label,
        'exposed': exposed,
        'not_exposed': not_exposed,
        'diff': exposed - not_exposed
    })

factor_df = pd.DataFrame(factor_rates).sort_values('diff', ascending=True)

# 수평 쌍봉 바 차트
y_pos = range(len(factor_df))
bars_safe = axes[1, 0].barh(
    [y - 0.2 for y in y_pos], factor_df['not_exposed'],
    height=0.35, label='Without Risk Factor',
    color='#3498db', edgecolor='black', linewidth=0.5, alpha=0.7
)
bars_risk = axes[1, 0].barh(
    [y + 0.2 for y in y_pos], factor_df['exposed'],
    height=0.35, label='With Risk Factor',
    color='#e74c3c', edgecolor='black', linewidth=0.5, alpha=0.7
)

# 차이 표시
for y, (_, row) in zip(y_pos, factor_df.iterrows()):
    axes[1, 0].text(
        row['exposed'] + 1, y + 0.2,
        f'+{row["diff"]:.1f}%p',
        va='center', fontsize=9, fontweight='bold', color='#c0392b'
    )

axes[1, 0].set_yticks(list(y_pos))
axes[1, 0].set_yticklabels(factor_df['factor'])
axes[1, 0].set_title('Individual Risk Factor Impact', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('Severe Accident Rate (%)', fontsize=10)
axes[1, 0].legend(loc='lower right', fontsize=9)

# --- (3-4) 위험 점수별 위험 요인 구성 비율 (스택 바) ---
# 각 Risk Score에서 어떤 위험 요인이 얼마나 활성화되었는지
risk_composition = df.groupby('Risk_Score')[
    ['Risk_Weather', 'Risk_Road', 'Risk_Density', 'Risk_Speed']
].mean() * 100

factor_labels = ['Bad Weather', 'Bad Road', 'High Density', 'High Speed']
factor_colors = ['#3498db', '#2ecc71', '#e74c3c', '#f39c12']

bottom = np.zeros(len(risk_composition))
for i, (col, label) in enumerate(zip(
    ['Risk_Weather', 'Risk_Road', 'Risk_Density', 'Risk_Speed'], factor_labels
)):
    vals = risk_composition[col].values
    axes[1, 1].bar(
        risk_composition.index, vals, bottom=bottom,
        label=label, color=factor_colors[i],
        edgecolor='white', linewidth=0.5, width=0.6
    )
    # 각 세그먼트에 비율 표시 (충분히 큰 경우만)
    for j, (v, b) in enumerate(zip(vals, bottom)):
        if v > 15:  # 15% 이상일 때만 표시
            axes[1, 1].text(
                risk_composition.index[j], b + v/2,
                f'{v:.0f}%', ha='center', va='center',
                fontsize=8, fontweight='bold', color='white'
            )
    bottom += vals

axes[1, 1].set_title('Risk Factor Composition by Score', fontsize=12, fontweight='bold')
axes[1, 1].set_xlabel('Risk Score', fontsize=10)
axes[1, 1].set_ylabel('Factor Activation Rate (%)', fontsize=10)
axes[1, 1].set_xticks(range(5))
axes[1, 1].legend(loc='upper left', fontsize=8, ncol=2)

plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 4) 결과 요약 출력
# ----------------------------------------------------------
print("\n" + "=" * 70)
print("복합 위험 점수 분석 결과 요약")
print("=" * 70)
print("\n[위험 요인 정의]")
print("  1. Bad Weather: 맑음(Clear)이 아닌 날씨")
print("  2. Bad Road: 건조(Dry)가 아닌 도로 상태")
print("  3. High Density: 교통 혼잡도 High(2)")
print("  4. High Speed: 제한 속도 80km/h 초과")

print(f"\n{'점수':<8} {'건수':<10} {'비율':<10} {'심각 건수':<10} {'심각 비율':<12} {'위험 수준'}")
print("-" * 65)
for _, row in risk_stats.iterrows():
    score = int(row['Risk_Score'])
    level = ['안전', '주의', '경고', '위험', '극도위험'][score]
    marker = ['🟢', '🟡', '🟠', '🔴', '🚨'][score]
    print(f'  {score}     {int(row["total"]):<10} {row["pct_of_total"]}%     '
          f'{int(row["severe"]):<10} {row["rate_pct"]}%       {marker} {level}')

# 임계점 분석
score_2plus = df[df['Risk_Score'] >= 2]
score_3plus = df[df['Risk_Score'] >= 3]
print(f"\n[임계점 분석]")
print(f"  위험 요인 2개 이상: 심각 비율 {score_2plus['Accident_Severity'].mean()*100:.1f}% (n={len(score_2plus)})")
print(f"  위험 요인 3개 이상: 심각 비율 {score_3plus['Accident_Severity'].mean()*100:.1f}% (n={len(score_3plus)})")

# ============================================================
# [Full EDA Report — Part 1/3] 초기 설정 + 기초통계 + 분포
# ============================================================
# 전처리된 데이터셋(dataset_traffic_accident_cleaned.csv) 기반으로
# 가능한 모든 탐색적 데이터 분석을 수행하여
# traffic_accident_full_EDA_report.pdf 파일로 저장합니다.
#
# Part 1: 데이터 개요, 기초통계, 타겟 분포, 범주형·수치형 분포
# Part 2: 상관관계, 카이제곱 검정, 변수-타겟 관계
# Part 3: 상호작용, 복합위험, 종합 대시보드, 저장
# ============================================================

# 경고 메시지 숨기기
import warnings
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib.patches import Patch
import matplotlib.gridspec as gridspec

# ----------------------------------------------------------
# 0) 전역 설정
# ----------------------------------------------------------
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

# 데이터 불러오기
data_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/dataset_traffic_accident_cleaned.csv'
df = pd.read_csv(data_path)

# 공통 색상
SEV = {0: '#3498db', 1: '#e74c3c'}
SEV_LBL = {0: 'Non-Severe (0)', 1: 'Severe (1)'}
TD_LBL = {0: 'Low (0)', 1: 'Medium (1)', 2: 'High (2)'}
TD_CLR = {0: '#2ecc71', 1: '#f39c12', 2: '#e74c3c'}

num_cols = ['Traffic_Density', 'Speed_Limit', 'Driver_Age', 'Driver_Experience', 'Accident_Severity']
cat_cols = ['Weather', 'Road_Type', 'Time_of_Day', 'Road_Condition', 'Vehicle_Type', 'Road_Light_Condition']

# PDF 저장 경로
pdf_path = '/content/drive/MyDrive/대학원/생성형AI를활용한데이터과학실무/3rd-week/data/traffic_accident_full_EDA_report.pdf'

# 헤더 유틸 함수
def page_header(fig, title, sub=''):
    fig.suptitle(title, fontsize=15, fontweight='bold', y=0.98)
    if sub:
        fig.text(0.5, 0.94, sub, ha='center', fontsize=9, style='italic', color='gray')

# PDF 객체 생성 (Part 3에서 닫음)
pdf = PdfPages(pdf_path)

# ============================================================
# PAGE 1: 데이터 개요 및 기초통계 요약
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(16, 10))
page_header(fig, 'Page 1. Dataset Overview & Descriptive Statistics',
            f'Cleaned Dataset: {df.shape[0]} rows × {df.shape[1]} columns | No missing values')

# (1-1) 데이터 타입 분포
dtype_counts = df.dtypes.map(lambda x: 'Numeric (int64)' if 'int' in str(x) else 'Categorical (object)')
dtype_vc = dtype_counts.value_counts()
axes[0, 0].pie(dtype_vc.values, labels=dtype_vc.index, autopct='%1.0f%%',
               colors=['#3498db', '#e74c3c'], startangle=90, textprops={'fontsize': 11})
axes[0, 0].set_title(f'Variable Types\n({len(df.columns)} total)', fontsize=12, fontweight='bold')

# (1-2) 수치형 변수 기초통계 표
stats_data = []
for col in num_cols:
    stats_data.append([col, f'{df[col].min()}', f'{df[col].mean():.2f}',
                       f'{df[col].median():.1f}', f'{df[col].max()}',
                       f'{df[col].std():.2f}', f'{df[col].skew():.2f}'])
axes[0, 1].axis('off')
table = axes[0, 1].table(
    cellText=stats_data,
    colLabels=['Variable', 'Min', 'Mean', 'Median', 'Max', 'Std', 'Skew'],
    loc='center', cellLoc='center'
)
table.auto_set_font_size(False)
table.set_fontsize(8.5)
table.scale(1, 1.4)
for (r, c), cell in table.get_celld().items():
    if r == 0:
        cell.set_facecolor('#34495e')
        cell.set_text_props(color='white', fontweight='bold')
    elif r % 2 == 0:
        cell.set_facecolor('#f8f9fa')
axes[0, 1].set_title('Numeric Variables Summary', fontsize=12, fontweight='bold', pad=20)

# (1-3) 범주형 변수 고유값 수
cat_unique = [(col, df[col].nunique(), df[col].mode()[0]) for col in cat_cols]
axes[1, 0].axis('off')
table2 = axes[1, 0].table(
    cellText=[(c, str(u), m, f'{df[df[c]==m].shape[0]} ({df[df[c]==m].shape[0]/len(df)*100:.1f}%)')
              for c, u, m in cat_unique],
    colLabels=['Variable', 'Unique', 'Mode', 'Mode Count (%)'],
    loc='center', cellLoc='center'
)
table2.auto_set_font_size(False)
table2.set_fontsize(8.5)
table2.scale(1, 1.4)
for (r, c), cell in table2.get_celld().items():
    if r == 0:
        cell.set_facecolor('#34495e')
        cell.set_text_props(color='white', fontweight='bold')
    elif r % 2 == 0:
        cell.set_facecolor('#f8f9fa')
axes[1, 0].set_title('Categorical Variables Summary', fontsize=12, fontweight='bold', pad=20)

# (1-4) 타겟 변수 분포
counts = df['Accident_Severity'].value_counts().sort_index()
bars = axes[1, 1].bar([SEV_LBL[i] for i in counts.index], counts.values,
                       color=[SEV[i] for i in counts.index], edgecolor='black', linewidth=0.5, width=0.5)
for bar, cnt in zip(bars, counts.values):
    axes[1, 1].text(bar.get_x()+bar.get_width()/2, bar.get_height()+8,
                     f'{cnt} ({cnt/len(df)*100:.1f}%)', ha='center', fontsize=10, fontweight='bold')
axes[1, 1].set_title('Target: Accident_Severity', fontsize=12, fontweight='bold')
axes[1, 1].set_ylabel('Count')
axes[1, 1].set_ylim(0, max(counts)*1.25)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 1 저장 완료")

# ============================================================
# PAGE 2: 전체 수치형 변수 분포 (히스토그램 + 박스플롯)
# ============================================================
fig, axes = plt.subplots(2, 5, figsize=(20, 8))
page_header(fig, 'Page 2. Numeric Variable Distributions',
            'Top: Histogram (Mean=Red, Median=Green) | Bottom: Box Plot by Severity')

for i, col in enumerate(num_cols):
    # 히스토그램
    axes[0, i].hist(df[col], bins=20, color='#5dade2', edgecolor='black', linewidth=0.5, alpha=0.8)
    axes[0, i].axvline(df[col].mean(), color='red', linestyle='--', lw=1.5, label=f'Mean:{df[col].mean():.1f}')
    axes[0, i].axvline(df[col].median(), color='green', linestyle='-.', lw=1.5, label=f'Med:{df[col].median():.1f}')
    axes[0, i].set_title(col, fontsize=10, fontweight='bold')
    axes[0, i].legend(fontsize=6)
    axes[0, i].tick_params(labelsize=7)

    # 박스플롯
    d0 = df[df['Accident_Severity']==0][col]
    d1 = df[df['Accident_Severity']==1][col]
    bp = axes[1, i].boxplot([d0, d1], labels=['Sev=0','Sev=1'], patch_artist=True, widths=0.5)
    bp['boxes'][0].set_facecolor(SEV[0]); bp['boxes'][0].set_alpha(0.6)
    bp['boxes'][1].set_facecolor(SEV[1]); bp['boxes'][1].set_alpha(0.6)
    axes[1, i].set_title(f'{col}\nby Severity', fontsize=9, fontweight='bold')
    axes[1, i].tick_params(labelsize=7)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 2 저장 완료")

# ============================================================
# PAGE 3: 범주형 변수 빈도 분포 (6개 변수)
# ============================================================
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
page_header(fig, 'Page 3. Categorical Variable Frequency Distributions')

palette_list = ['#5dade2','#58d68d','#f0b27a','#bb8fce','#f1948a','#85c1e9']
for i, col in enumerate(cat_cols):
    ax = axes[i//3, i%3]
    vc = df[col].value_counts()
    bars = ax.barh(vc.index[::-1], vc.values[::-1], color=palette_list[i],
                    edgecolor='black', linewidth=0.5)
    for bar, v in zip(bars, vc.values[::-1]):
        ax.text(bar.get_width()+3, bar.get_y()+bar.get_height()/2,
                f'{v} ({v/len(df)*100:.1f}%)', va='center', fontsize=9)
    ax.set_title(col, fontsize=12, fontweight='bold')
    ax.set_xlabel('Count')

plt.tight_layout(rect=[0, 0, 1, 0.94])
pdf.savefig(fig); plt.close()
print("✔ Page 3 저장 완료")

# ============================================================
# PAGE 4: 범주형 변수별 심각 사고 비율
# ============================================================
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
page_header(fig, 'Page 4. Severe Accident Rate by Categorical Variables',
            'Stacked 100% bars — Red portion = Severe Rate')

for i, col in enumerate(cat_cols):
    ax = axes[i//3, i%3]
    order = df.groupby(col)['Accident_Severity'].mean().sort_values(ascending=False).index
    ct = pd.crosstab(df[col], df['Accident_Severity'], normalize='index') * 100
    ct = ct.loc[order]
    ct.plot(kind='barh', stacked=True, ax=ax, color=[SEV[0], SEV[1]],
            edgecolor='black', linewidth=0.5)
    for j, (idx, row) in enumerate(ct.iterrows()):
        if row[1] > 5:
            ax.text(row[0]+row[1]/2, j, f'{row[1]:.1f}%', ha='center', va='center',
                    fontsize=9, fontweight='bold', color='white')
    ax.set_title(col, fontsize=12, fontweight='bold')
    ax.set_xlabel('%')
    ax.set_ylabel('')
    ax.legend(['Non-Severe','Severe'], loc='lower right', fontsize=8)

plt.tight_layout(rect=[0, 0, 1, 0.94])
pdf.savefig(fig); plt.close()
print("✔ Page 4 저장 완료")

# ============================================================
# PAGE 5: 수치형 변수별 심각 사고 상세 비교
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
page_header(fig, 'Page 5. Key Numeric Variables vs Severity (Detail)')

# (5-1) Traffic_Density 그룹바
td_ct = pd.crosstab(df['Traffic_Density'], df['Accident_Severity'])
x = range(len(td_ct)); w = 0.35
axes[0,0].bar([i-w/2 for i in x], td_ct[0], w, label='Non-Severe', color=SEV[0], edgecolor='black', lw=0.5)
axes[0,0].bar([i+w/2 for i in x], td_ct[1], w, label='Severe', color=SEV[1], edgecolor='black', lw=0.5)
axes[0,0].set_xticks(list(x)); axes[0,0].set_xticklabels([TD_LBL[i] for i in td_ct.index])
axes[0,0].set_title('Traffic Density — Count by Severity', fontsize=12, fontweight='bold')
axes[0,0].legend(fontsize=9); axes[0,0].set_ylabel('Count')

# (5-2) Traffic_Density 심각비율
td_pct = pd.crosstab(df['Traffic_Density'], df['Accident_Severity'], normalize='index')*100
bars = axes[0,1].bar([TD_LBL[i] for i in td_pct.index], td_pct[1],
                      color=[TD_CLR[i] for i in td_pct.index], edgecolor='black', lw=0.5, width=0.5)
for bar, v in zip(bars, td_pct[1]):
    axes[0,1].text(bar.get_x()+bar.get_width()/2, bar.get_height()+1,
                    f'{v:.1f}%', ha='center', fontsize=11, fontweight='bold')
axes[0,1].set_title('Traffic Density — Severe Rate', fontsize=12, fontweight='bold')
axes[0,1].set_ylabel('Severe Rate (%)'); axes[0,1].set_ylim(0, td_pct[1].max()*1.3)

# (5-3) Speed_Limit 히스토그램
for s in [0, 1]:
    sub = df[df['Accident_Severity']==s]['Speed_Limit']
    axes[1,0].hist(sub, bins=15, alpha=0.6, label=SEV_LBL[s], color=SEV[s], edgecolor='black', lw=0.5)
    axes[1,0].axvline(sub.mean(), color=SEV[s], linestyle='--', lw=2, alpha=0.8)
axes[1,0].set_title('Speed Limit — Distribution by Severity', fontsize=12, fontweight='bold')
axes[1,0].set_xlabel('Speed Limit'); axes[1,0].set_ylabel('Count'); axes[1,0].legend(fontsize=9)

# (5-4) Speed_Limit 박스플롯
df_tmp = df.copy(); df_tmp['Sev'] = df_tmp['Accident_Severity'].map(SEV_LBL)
sns.boxplot(data=df_tmp, x='Sev', y='Speed_Limit', palette=[SEV[0], SEV[1]], ax=axes[1,1])
for s in [0, 1]:
    m = df[df['Accident_Severity']==s]['Speed_Limit'].mean()
    axes[1,1].text(s, m+2, f'Mean:{m:.1f}', ha='center', fontsize=9, fontweight='bold', color=SEV[s])
axes[1,1].set_title('Speed Limit — Box Plot', fontsize=12, fontweight='bold')
axes[1,1].set_xlabel(''); axes[1,1].set_ylabel('Speed Limit')

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 5 저장 완료")
print("\n>>> Part 1 완료 — Part 2 셀을 실행하세요 >>>")

# ============================================================
# [Full EDA Report — Part 2/3] 상관관계 + 카이제곱 + 상호작용
# ============================================================
# Part 1에서 생성한 pdf 객체에 이어서 페이지를 추가합니다.
# ============================================================

# ============================================================
# PAGE 6: 피어슨 상관계수 히트맵 + 다중공선성
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
page_header(fig, 'Page 6. Correlation Analysis',
            'Pearson Correlation (left) & Multicollinearity Check (right)')

# (6-1) 히트맵
corr_m = df[num_cols].corr().round(4)
sns.heatmap(corr_m, annot=True, fmt='.4f', cmap='RdBu_r', center=0, vmin=-1, vmax=1,
            square=True, linewidths=0.5, ax=axes[0],
            cbar_kws={'shrink': 0.8, 'label': 'Pearson r'})
axes[0].set_title('Pearson Correlation Heatmap', fontsize=12, fontweight='bold')
axes[0].set_xticklabels(axes[0].get_xticklabels(), rotation=45, ha='right')

# (6-2) Driver_Age vs Experience 산점도
for s in [0, 1]:
    sub = df[df['Accident_Severity']==s]
    axes[1].scatter(sub['Driver_Age'], sub['Driver_Experience'], c=SEV[s],
                     label=SEV_LBL[s], alpha=0.35, s=25, edgecolors='white', linewidth=0.3)
z = np.polyfit(df['Driver_Age'], df['Driver_Experience'], 1)
p = np.poly1d(z)
xr = np.linspace(df['Driver_Age'].min(), df['Driver_Age'].max(), 100)
rv = df['Driver_Age'].corr(df['Driver_Experience'])
axes[1].plot(xr, p(xr), '--', color='black', lw=2, alpha=0.7, label=f'r = {rv:.4f}')
axes[1].set_title('Driver Age vs Experience\n(Multicollinearity)', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Driver Age'); axes[1].set_ylabel('Driver Experience (years)')
axes[1].legend(fontsize=9)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 6 저장 완료")

# ============================================================
# PAGE 7: 타겟과의 연관 강도 종합 (Pearson r + Cramer's V)
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
page_header(fig, 'Page 7. Variable Association Strength with Target',
            "|Pearson r| for numeric, Cramer's V for categorical")

# 수치형: |Pearson r|
num_assoc = {c: abs(df[c].corr(df['Accident_Severity'])) for c in num_cols if c != 'Accident_Severity'}
# 범주형: Cramer's V
cat_assoc = {}
for col in cat_cols:
    ct = pd.crosstab(df[col], df['Accident_Severity'])
    chi2, pv, dof, _ = stats.chi2_contingency(ct)
    n = ct.sum().sum(); md = min(ct.shape) - 1
    cat_assoc[col] = np.sqrt(chi2 / (n * md))

all_a = {**num_assoc, **cat_assoc}
adf = pd.DataFrame({'Var': list(all_a.keys()), 'Assoc': list(all_a.values()),
                     'Type': ['Numeric']*len(num_assoc) + ['Categorical']*len(cat_assoc)})
adf = adf.sort_values('Assoc', ascending=True)

# (7-1) 수평 바 차트
bc = ['#e74c3c' if t == 'Categorical' else '#3498db' for t in adf['Type']]
bars = axes[0].barh(adf['Var'], adf['Assoc'], color=bc, edgecolor='black', lw=0.5, height=0.6)
for bar, v in zip(bars, adf['Assoc']):
    axes[0].text(bar.get_width()+0.008, bar.get_y()+bar.get_height()/2,
                  f'{v:.4f}', ha='left', va='center', fontsize=9, fontweight='bold')
axes[0].axvline(0.1, color='gray', ls=':', lw=1, alpha=0.6)
axes[0].axvline(0.3, color='gray', ls=':', lw=1, alpha=0.6)
axes[0].legend(handles=[Patch(fc='#3498db', ec='black', label='Numeric (|r|)'),
                         Patch(fc='#e74c3c', ec='black', label="Categorical (Cramer's V)")],
               loc='lower right', fontsize=9)
axes[0].set_title('All Variables — Association Strength', fontsize=12, fontweight='bold')
axes[0].set_xlabel('Association'); axes[0].set_xlim(0, 0.7)

# (7-2) 카이제곱 검정 결과 표
chi_data = []
for col in cat_cols:
    ct = pd.crosstab(df[col], df['Accident_Severity'])
    chi2, pv, dof, _ = stats.chi2_contingency(ct)
    n = ct.sum().sum(); md = min(ct.shape) - 1
    cv = np.sqrt(chi2 / (n * md))
    sig = '✔ Yes' if pv < 0.05 else '✘ No'
    pstr = f'{pv:.2e}' if pv < 0.001 else f'{pv:.4f}'
    chi_data.append([col, f'{chi2:.1f}', pstr, f'{cv:.4f}', sig])

axes[1].axis('off')
tbl = axes[1].table(cellText=chi_data,
                     colLabels=['Variable', 'Chi²', 'p-value', "Cramer's V", 'Significant?'],
                     loc='center', cellLoc='center')
tbl.auto_set_font_size(False); tbl.set_fontsize(9); tbl.scale(1, 1.6)
for (r, c), cell in tbl.get_celld().items():
    if r == 0:
        cell.set_facecolor('#34495e'); cell.set_text_props(color='white', fontweight='bold')
    elif r % 2 == 0:
        cell.set_facecolor('#f8f9fa')
    # 유의한 행 강조
    if r > 0 and c == 4:
        txt = cell.get_text().get_text()
        if '✔' in txt:
            cell.set_text_props(color='#c0392b', fontweight='bold')
axes[1].set_title('Chi-Square Test Results (α=0.05)', fontsize=12, fontweight='bold', pad=20)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 7 저장 완료")

# ============================================================
# PAGE 8: Weather × Road_Condition 상호작용
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
page_header(fig, 'Page 8. Interaction: Weather × Road Condition',
            'Frequency (left) & Severe Rate % (right)')

wo = ['Rainy','Stormy','Snowy','Foggy','Clear']
ro = ['Wet','Icy','Under Construction','Dry']

cc = pd.crosstab(df['Weather'], df['Road_Condition']).reindex(index=wo, columns=ro)
sns.heatmap(cc, annot=True, fmt='d', cmap='YlOrRd', linewidths=0.5, ax=axes[0],
            cbar_kws={'label': 'Count'})
axes[0].set_title('Cross-tabulation (Frequency)', fontsize=12, fontweight='bold')

cs = (df.groupby(['Weather','Road_Condition'])['Accident_Severity'].mean()*100).unstack()
cs = cs.reindex(index=wo, columns=ro)
sns.heatmap(cs, annot=True, fmt='.1f', cmap='RdYlGn_r', linewidths=0.5, ax=axes[1],
            vmin=0, vmax=100, cbar_kws={'label': 'Severe Rate (%)'})
axes[1].set_title('Severe Accident Rate (%)', fontsize=12, fontweight='bold')

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 8 저장 완료")

# ============================================================
# PAGE 9: Traffic_Density × Speed_Limit 상호작용
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
page_header(fig, 'Page 9. Interaction: Traffic Density × Speed Limit',
            'Box Plot (left) & Severe Rate by Speed Group × Density (right)')

dp = df.copy()
dp['TD'] = dp['Traffic_Density'].map(TD_LBL)
dp['Sev'] = dp['Accident_Severity'].map(SEV_LBL)
sns.boxplot(data=dp, x='TD', y='Speed_Limit', hue='Sev',
            palette=[SEV[0], SEV[1]], ax=axes[0],
            order=['Low (0)','Medium (1)','High (2)'])
axes[0].set_title('Speed Limit by Density & Severity', fontsize=12, fontweight='bold')
axes[0].legend(title='Severity', fontsize=9)

dp['SG'] = pd.cut(dp['Speed_Limit'], bins=[0,50,80,110],
                    labels=['Low(<=50)','Mid(51-80)','High(81-110)'])
ir = dp.groupby(['SG','Traffic_Density'])['Accident_Severity'].mean()*100
ip = ir.unstack(); ip.columns = [TD_LBL[c] for c in ip.columns]
ip.plot(kind='bar', ax=axes[1], color=[TD_CLR[0],TD_CLR[1],TD_CLR[2]],
        edgecolor='black', lw=0.5)
for cont in axes[1].containers:
    axes[1].bar_label(cont, fmt='%.1f%%', fontsize=8, padding=2)
axes[1].set_title('Severe Rate: Speed × Density', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Speed Group'); axes[1].set_ylabel('Severe Rate (%)')
axes[1].set_xticklabels(axes[1].get_xticklabels(), rotation=0)
axes[1].legend(title='Traffic Density', fontsize=9); axes[1].set_ylim(0, 100)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 9 저장 완료")

# ============================================================
# PAGE 10: Traffic_Density 증폭 효과 — 패싯 히트맵
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(20, 6))
page_header(fig, 'Page 10. Traffic Density Amplification Effect',
            'Severe Rate by Weather × Road Condition at Low / Medium / High Density')

for idx, td in enumerate([0, 1, 2]):
    ds = df[df['Traffic_Density']==td]
    cr = (ds.groupby(['Weather','Road_Condition'])['Accident_Severity'].mean()*100).unstack()
    cn = ds.groupby(['Weather','Road_Condition'])['Accident_Severity'].count().unstack()
    aw = [w for w in wo if w in cr.index]
    ar = [r for r in ro if r in cr.columns]
    cr = cr.reindex(index=aw, columns=ar)
    cn = cn.reindex(index=aw, columns=ar)
    ann = cr.copy().astype(str)
    for w in aw:
        for r in ar:
            rv2 = cr.loc[w,r] if not pd.isna(cr.loc[w,r]) else np.nan
            cv2 = cn.loc[w,r] if not pd.isna(cn.loc[w,r]) else 0
            ann.loc[w,r] = f'{rv2:.0f}%\n(n={int(cv2)})' if not pd.isna(rv2) and cv2>0 else '-'
    sns.heatmap(cr, annot=ann, fmt='', cmap='RdYlGn_r', vmin=0, vmax=100,
                linewidths=0.8, linecolor='white', ax=axes[idx],
                cbar=(idx==2), cbar_kws={'label':'Severe %','shrink':0.8} if idx==2 else {},
                annot_kws={'fontsize': 9})
    sr2 = ds['Accident_Severity'].mean()*100
    axes[idx].set_title(f'Density: {TD_LBL[td]}\n(n={len(ds)}, Avg:{sr2:.1f}%)',
                         fontsize=11, fontweight='bold')
    axes[idx].set_xlabel('Road Condition'); axes[idx].set_ylabel('Weather' if idx==0 else '')
    if idx > 0: axes[idx].set_yticklabels([])

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 10 저장 완료")
print("\n>>> Part 2 완료 — Part 3 셀을 실행하세요 >>>")

# ============================================================
# [Full EDA Report — Part 3/3] 차량유형 + 복합위험 + Top10 + 요약
# ============================================================
# Part 1·2에서 생성한 pdf 객체에 이어서 페이지를 추가하고,
# 최종적으로 pdf.close()로 파일을 완성합니다.
# ============================================================

# ============================================================
# PAGE 11: 차량 유형별 숨겨진 패턴 (조건별 비교)
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(16, 11))
page_header(fig, 'Page 11. Vehicle Type — Hidden Patterns by Conditions',
            'Severe Rate comparison across increasingly restrictive conditions')

# 차량 유형 순서 및 색상
vehicle_order = ['Car', 'Truck', 'Bus']
vehicle_colors = {'Car': '#3498db', 'Truck': '#e67e22', 'Bus': '#e74c3c'}

# 4가지 조건 정의
conditions = [
    {'title': 'Overall (All Data)', 'data': df},
    {'title': 'Highway Only', 'data': df[df['Road_Type'] == 'Highway']},
    {'title': 'Highway + High Density',
     'data': df[(df['Road_Type'] == 'Highway') & (df['Traffic_Density'] == 2)]},
    {'title': 'Highway + High Density\n+ Wet/Rainy',
     'data': df[(df['Road_Type'] == 'Highway') & (df['Traffic_Density'] == 2) &
               ((df['Road_Condition'] == 'Wet') | (df['Weather'] == 'Rainy'))]}
]

for ci, cond in enumerate(conditions):
    ax = axes[ci // 2, ci % 2]
    df_sub = cond['data']
    results = []
    for vtype in vehicle_order:
        sub = df_sub[df_sub['Vehicle_Type'] == vtype]
        total = len(sub)
        severe = int(sub['Accident_Severity'].sum()) if total > 0 else 0
        rate = severe / total * 100 if total > 0 else 0
        results.append({'type': vtype, 'total': total, 'severe': severe, 'rate': rate})

    bars = ax.bar([r['type'] for r in results], [r['rate'] for r in results],
                  color=[vehicle_colors[r['type']] for r in results],
                  edgecolor='black', linewidth=0.5, width=0.55)
    for bar, r in zip(bars, results):
        if r['total'] > 0:
            ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1.5,
                    f'{r["rate"]:.1f}%\n({r["severe"]}/{r["total"]})',
                    ha='center', va='bottom', fontsize=9, fontweight='bold')
    if len(df_sub) > 0:
        overall = df_sub['Accident_Severity'].mean() * 100
        ax.axhline(y=overall, color='gray', linestyle='--', linewidth=1.5, alpha=0.7)
        ax.text(2.4, overall + 1, f'Avg: {overall:.1f}%', fontsize=8, color='gray', fontweight='bold')
    ax.set_title(f'{cond["title"]}\n(n={len(df_sub)})', fontsize=11, fontweight='bold')
    ax.set_ylabel('Severe Rate (%)')
    ax.set_ylim(0, min(max([r['rate'] for r in results]) * 1.4 + 5, 105))

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 11 저장 완료")

# ============================================================
# PAGE 12: 복합 위험 점수 분석 (4-패널 대시보드)
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
page_header(fig, 'Page 12. Compound Risk Score Analysis',
            'Risk Factors: Bad Weather + Bad Road + High Density + High Speed (>80)')

# 위험 요인 플래그 생성
df['Risk_Weather'] = (df['Weather'] != 'Clear').astype(int)
df['Risk_Road'] = (df['Road_Condition'] != 'Dry').astype(int)
df['Risk_Density'] = (df['Traffic_Density'] == 2).astype(int)
df['Risk_Speed'] = (df['Speed_Limit'] > 80).astype(int)
df['Risk_Score'] = df['Risk_Weather'] + df['Risk_Road'] + df['Risk_Density'] + df['Risk_Speed']

# 위험 점수별 통계
risk_stats = df.groupby('Risk_Score').agg(
    total=('Accident_Severity', 'count'),
    severe=('Accident_Severity', 'sum'),
    rate=('Accident_Severity', 'mean')
).reset_index()
risk_stats['rate_pct'] = (risk_stats['rate'] * 100).round(1)
risk_stats['pct_of_total'] = (risk_stats['total'] / len(df) * 100).round(1)

# (12-1) 위험 점수별 심각 사고 비율
colors_gradient = ['#27ae60', '#f1c40f', '#e67e22', '#e74c3c', '#8e44ad']
bars = axes[0, 0].bar(risk_stats['Risk_Score'], risk_stats['rate_pct'],
                       color=[colors_gradient[i] for i in risk_stats['Risk_Score']],
                       edgecolor='black', linewidth=0.8, width=0.6)
for bar, (_, row) in zip(bars, risk_stats.iterrows()):
    axes[0, 0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2,
                     f'{row["rate_pct"]}%\n({int(row["severe"])}/{int(row["total"])})',
                     ha='center', va='bottom', fontsize=9, fontweight='bold')
avg_rate = df['Accident_Severity'].mean() * 100
axes[0, 0].axhline(y=avg_rate, color='gray', linestyle='--', linewidth=1.5, alpha=0.7)
axes[0, 0].text(4.3, avg_rate + 1, f'Avg: {avg_rate:.1f}%', fontsize=8, color='gray', fontweight='bold')
axes[0, 0].axhline(y=50, color='red', linestyle=':', linewidth=1, alpha=0.5)
axes[0, 0].set_title('Severe Rate by Risk Score', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('Number of Risk Factors (0~4)')
axes[0, 0].set_ylabel('Severe Accident Rate (%)')
axes[0, 0].set_xticks(range(5))
axes[0, 0].set_xticklabels(['0\n(Safe)', '1', '2', '3', '4\n(Most\nDangerous)'])
axes[0, 0].set_ylim(0, max(risk_stats['rate_pct']) * 1.25)

# (12-2) 건수 분포
bars2 = axes[0, 1].bar(risk_stats['Risk_Score'], risk_stats['total'],
                        color=[colors_gradient[i] for i in risk_stats['Risk_Score']],
                        edgecolor='black', linewidth=0.8, width=0.6, alpha=0.8)
for bar, (_, row) in zip(bars2, risk_stats.iterrows()):
    axes[0, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 5,
                     f'{int(row["total"])}\n({row["pct_of_total"]}%)',
                     ha='center', va='bottom', fontsize=9, fontweight='bold')
axes[0, 1].set_title('Data Distribution by Risk Score', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('Number of Risk Factors (0~4)')
axes[0, 1].set_ylabel('Number of Accidents')
axes[0, 1].set_xticks(range(5))
axes[0, 1].set_ylim(0, max(risk_stats['total']) * 1.25)

# (12-3) 개별 위험 요인 기여도
risk_factors = {
    'Bad Weather\n(Not Clear)': 'Risk_Weather',
    'Bad Road\n(Not Dry)': 'Risk_Road',
    'High Density\n(TD=2)': 'Risk_Density',
    'High Speed\n(>80 km/h)': 'Risk_Speed'
}
factor_rates = []
for label, col in risk_factors.items():
    exposed = df[df[col] == 1]['Accident_Severity'].mean() * 100
    not_exposed = df[df[col] == 0]['Accident_Severity'].mean() * 100
    factor_rates.append({'factor': label, 'exposed': exposed,
                         'not_exposed': not_exposed, 'diff': exposed - not_exposed})
factor_df = pd.DataFrame(factor_rates).sort_values('diff', ascending=True)

y_pos = range(len(factor_df))
axes[1, 0].barh([y - 0.2 for y in y_pos], factor_df['not_exposed'], height=0.35,
                 label='Without Factor', color='#3498db', edgecolor='black', lw=0.5, alpha=0.7)
axes[1, 0].barh([y + 0.2 for y in y_pos], factor_df['exposed'], height=0.35,
                 label='With Factor', color='#e74c3c', edgecolor='black', lw=0.5, alpha=0.7)
for y, (_, row) in zip(y_pos, factor_df.iterrows()):
    axes[1, 0].text(row['exposed'] + 1, y + 0.2, f'+{row["diff"]:.1f}%p',
                     va='center', fontsize=9, fontweight='bold', color='#c0392b')
axes[1, 0].set_yticks(list(y_pos))
axes[1, 0].set_yticklabels(factor_df['factor'])
axes[1, 0].set_title('Individual Risk Factor Impact', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('Severe Accident Rate (%)')
axes[1, 0].legend(loc='lower right', fontsize=9)

# (12-4) 위험 점수별 요인 구성 비율 (스택 바)
risk_composition = df.groupby('Risk_Score')[
    ['Risk_Weather', 'Risk_Road', 'Risk_Density', 'Risk_Speed']
].mean() * 100
factor_labels = ['Bad Weather', 'Bad Road', 'High Density', 'High Speed']
factor_colors = ['#3498db', '#2ecc71', '#e74c3c', '#f39c12']
bottom = np.zeros(len(risk_composition))
for i, (col, label) in enumerate(zip(
    ['Risk_Weather', 'Risk_Road', 'Risk_Density', 'Risk_Speed'], factor_labels
)):
    vals = risk_composition[col].values
    axes[1, 1].bar(risk_composition.index, vals, bottom=bottom,
                    label=label, color=factor_colors[i],
                    edgecolor='white', linewidth=0.5, width=0.6)
    for j, (v, b) in enumerate(zip(vals, bottom)):
        if v > 15:
            axes[1, 1].text(risk_composition.index[j], b + v/2, f'{v:.0f}%',
                             ha='center', va='center', fontsize=8, fontweight='bold', color='white')
    bottom += vals
axes[1, 1].set_title('Risk Factor Composition by Score', fontsize=12, fontweight='bold')
axes[1, 1].set_xlabel('Risk Score')
axes[1, 1].set_ylabel('Factor Activation Rate (%)')
axes[1, 1].set_xticks(range(5))
axes[1, 1].legend(loc='upper left', fontsize=8, ncol=2)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 12 저장 완료")

# ============================================================
# PAGE 13: 고위험 조합 Top 10
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
page_header(fig, 'Page 13. Top 10 High-Risk Combinations',
            'Highest Severe Rate combinations (min 5 observations)')

# 모든 범주형+Traffic_Density 조합에 대한 심각 비율 계산
combo_cols = ['Weather', 'Road_Condition', 'Traffic_Density']
combo = df.groupby(combo_cols).agg(
    total=('Accident_Severity', 'count'),
    severe=('Accident_Severity', 'sum'),
    rate=('Accident_Severity', 'mean')
).reset_index()
combo['rate_pct'] = combo['rate'] * 100
combo['TD_label'] = combo['Traffic_Density'].map(TD_LBL)
combo['label'] = combo['Weather'] + ' + ' + combo['Road_Condition'] + '\n(TD=' + combo['TD_label'] + ')'

# 최소 5건 이상인 조합만 필터링 후 상위 10개
combo_top = combo[combo['total'] >= 5].sort_values('rate_pct', ascending=False).head(10)
combo_top = combo_top.sort_values('rate_pct', ascending=True)  # 수평 바 차트용 역순

# (13-1) 수평 바 차트
cmap_risk = plt.cm.RdYlGn_r(np.linspace(0.3, 1.0, len(combo_top)))
bars = axes[0].barh(range(len(combo_top)), combo_top['rate_pct'],
                     color=cmap_risk, edgecolor='black', linewidth=0.5, height=0.6)
for i, (bar, (_, row)) in enumerate(zip(bars, combo_top.iterrows())):
    axes[0].text(bar.get_width() + 1, bar.get_y() + bar.get_height()/2,
                  f'{row["rate_pct"]:.1f}% ({int(row["severe"])}/{int(row["total"])})',
                  ha='left', va='center', fontsize=9, fontweight='bold')
axes[0].set_yticks(range(len(combo_top)))
axes[0].set_yticklabels(combo_top['label'], fontsize=8)
axes[0].set_title('Top 10 Highest Severe-Rate Combinations', fontsize=12, fontweight='bold')
axes[0].set_xlabel('Severe Accident Rate (%)')
axes[0].axvline(avg_rate, color='gray', linestyle='--', lw=1.5, alpha=0.7)
axes[0].text(avg_rate + 1, 0, f'Avg: {avg_rate:.1f}%', fontsize=8, color='gray')
axes[0].set_xlim(0, 110)

# (13-2) 상세 표
combo_bottom = combo[combo['total'] >= 5].sort_values('rate_pct', ascending=False).head(10)
tbl_data = []
for rank, (_, row) in enumerate(combo_bottom.iterrows(), 1):
    tbl_data.append([
        f'#{rank}', row['Weather'], row['Road_Condition'],
        row['TD_label'], f'{int(row["total"])}',
        f'{int(row["severe"])}', f'{row["rate_pct"]:.1f}%'
    ])

axes[1].axis('off')
tbl = axes[1].table(
    cellText=tbl_data,
    colLabels=['Rank', 'Weather', 'Road Cond.', 'Density', 'Total', 'Severe', 'Rate'],
    loc='center', cellLoc='center'
)
tbl.auto_set_font_size(False); tbl.set_fontsize(9); tbl.scale(1, 1.5)
for (r, c), cell in tbl.get_celld().items():
    if r == 0:
        cell.set_facecolor('#34495e'); cell.set_text_props(color='white', fontweight='bold')
    elif r % 2 == 0:
        cell.set_facecolor('#f8f9fa')
    # 상위 3개 강조
    if r > 0 and r <= 3:
        cell.set_facecolor('#fde8e8')
    if r > 0 and c == 6:
        cell.set_text_props(fontweight='bold', color='#c0392b')
axes[1].set_title('Detailed Breakdown', fontsize=12, fontweight='bold', pad=20)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 13 저장 완료")

# ============================================================
# PAGE 14: Driver Age × Experience 상세 분석
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(16, 11))
page_header(fig, 'Page 14. Driver Age & Experience — Detailed Analysis',
            'Age groups, Experience groups, and their interaction with Severity')

# (14-1) Driver_Age 그룹별 심각 비율
age_bins = [17, 25, 35, 45, 55, 65, 80]
age_labels = ['18-25', '26-35', '36-45', '46-55', '56-65', '66+']
df['Age_Group'] = pd.cut(df['Driver_Age'], bins=age_bins, labels=age_labels)
age_rate = df.groupby('Age_Group')['Accident_Severity'].agg(['mean', 'count']).reset_index()
age_rate['pct'] = age_rate['mean'] * 100

bars_age = axes[0, 0].bar(age_rate['Age_Group'].astype(str), age_rate['pct'],
                            color='#5dade2', edgecolor='black', lw=0.5, width=0.6)
for bar, (_, row) in zip(bars_age, age_rate.iterrows()):
    axes[0, 0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                     f'{row["pct"]:.1f}%\n(n={int(row["count"])})',
                     ha='center', fontsize=9, fontweight='bold')
axes[0, 0].axhline(avg_rate, color='gray', ls='--', lw=1.5, alpha=0.7)
axes[0, 0].set_title('Severe Rate by Age Group', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('Age Group'); axes[0, 0].set_ylabel('Severe Rate (%)')
axes[0, 0].set_ylim(0, age_rate['pct'].max() * 1.3)

# (14-2) Driver_Experience 그룹별 심각 비율
exp_bins = [-1, 5, 10, 20, 30, 50]
exp_labels = ['0-5', '6-10', '11-20', '21-30', '31+']
df['Exp_Group'] = pd.cut(df['Driver_Experience'], bins=exp_bins, labels=exp_labels)
exp_rate = df.groupby('Exp_Group')['Accident_Severity'].agg(['mean', 'count']).reset_index()
exp_rate['pct'] = exp_rate['mean'] * 100

bars_exp = axes[0, 1].bar(exp_rate['Exp_Group'].astype(str), exp_rate['pct'],
                            color='#58d68d', edgecolor='black', lw=0.5, width=0.6)
for bar, (_, row) in zip(bars_exp, exp_rate.iterrows()):
    axes[0, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                     f'{row["pct"]:.1f}%\n(n={int(row["count"])})',
                     ha='center', fontsize=9, fontweight='bold')
axes[0, 1].axhline(avg_rate, color='gray', ls='--', lw=1.5, alpha=0.7)
axes[0, 1].set_title('Severe Rate by Experience Group', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('Experience (years)'); axes[0, 1].set_ylabel('Severe Rate (%)')
axes[0, 1].set_ylim(0, exp_rate['pct'].max() * 1.3)

# (14-3) Age Group × Density 히트맵
age_density = (df.groupby(['Age_Group', 'Traffic_Density'])['Accident_Severity'].mean() * 100).unstack()
age_density.columns = [TD_LBL[c] for c in age_density.columns]
sns.heatmap(age_density, annot=True, fmt='.1f', cmap='YlOrRd', linewidths=0.5,
            ax=axes[1, 0], vmin=0, vmax=80, cbar_kws={'label': 'Severe %', 'shrink': 0.8})
axes[1, 0].set_title('Severe Rate: Age Group × Density', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('Traffic Density'); axes[1, 0].set_ylabel('Age Group')

# (14-4) Experience Group × Density 히트맵
exp_density = (df.groupby(['Exp_Group', 'Traffic_Density'])['Accident_Severity'].mean() * 100).unstack()
exp_density.columns = [TD_LBL[c] for c in exp_density.columns]
sns.heatmap(exp_density, annot=True, fmt='.1f', cmap='YlOrRd', linewidths=0.5,
            ax=axes[1, 1], vmin=0, vmax=80, cbar_kws={'label': 'Severe %', 'shrink': 0.8})
axes[1, 1].set_title('Severe Rate: Exp Group × Density', fontsize=12, fontweight='bold')
axes[1, 1].set_xlabel('Traffic Density'); axes[1, 1].set_ylabel('Experience Group')

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 14 저장 완료")

# ============================================================
# PAGE 15: Road_Type × Time_of_Day 상호작용
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
page_header(fig, 'Page 15. Interaction: Road Type × Time of Day',
            'Frequency (left) & Severe Rate % (right)')

road_order = ['Highway', 'City Road', 'Rural Road', 'Mountain Road']
time_order = ['Morning', 'Afternoon', 'Evening', 'Night']

# (15-1) 빈도 히트맵
ct_rt = pd.crosstab(df['Road_Type'], df['Time_of_Day'])
ct_rt = ct_rt.reindex(index=[r for r in road_order if r in ct_rt.index],
                       columns=[t for t in time_order if t in ct_rt.columns])
sns.heatmap(ct_rt, annot=True, fmt='d', cmap='Blues', linewidths=0.5, ax=axes[0],
            cbar_kws={'label': 'Count'})
axes[0].set_title('Cross-tabulation (Frequency)', fontsize=12, fontweight='bold')

# (15-2) 심각 비율 히트맵
sr_rt = (df.groupby(['Road_Type', 'Time_of_Day'])['Accident_Severity'].mean() * 100).unstack()
sr_rt = sr_rt.reindex(index=[r for r in road_order if r in sr_rt.index],
                       columns=[t for t in time_order if t in sr_rt.columns])
sns.heatmap(sr_rt, annot=True, fmt='.1f', cmap='RdYlGn_r', linewidths=0.5, ax=axes[1],
            vmin=0, vmax=80, cbar_kws={'label': 'Severe Rate (%)'})
axes[1].set_title('Severe Accident Rate (%)', fontsize=12, fontweight='bold')

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 15 저장 완료")

# ============================================================
# PAGE 16: Road_Light_Condition 분석
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
page_header(fig, 'Page 16. Road Light Condition Analysis',
            'Distribution and Severe Rate by Light Condition')

# (16-1) 조명 조건별 분포 (Severity별 그룹 바)
light_ct = pd.crosstab(df['Road_Light_Condition'], df['Accident_Severity'])
light_order = light_ct.index.tolist()
x = range(len(light_ct)); w = 0.35
axes[0].bar([i - w/2 for i in x], light_ct[0], w, label='Non-Severe',
             color=SEV[0], edgecolor='black', lw=0.5)
axes[0].bar([i + w/2 for i in x], light_ct[1], w, label='Severe',
             color=SEV[1], edgecolor='black', lw=0.5)
axes[0].set_xticks(list(x))
axes[0].set_xticklabels(light_order, rotation=15, ha='right')
axes[0].set_title('Count by Light Condition & Severity', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=9); axes[0].set_ylabel('Count')

# (16-2) 조명 조건별 심각 비율
light_rate = df.groupby('Road_Light_Condition')['Accident_Severity'].agg(['mean', 'count']).reset_index()
light_rate['pct'] = light_rate['mean'] * 100
light_rate = light_rate.sort_values('pct', ascending=True)

bars_lr = axes[1].barh(light_rate['Road_Light_Condition'], light_rate['pct'],
                         color='#f0b27a', edgecolor='black', lw=0.5, height=0.5)
for bar, (_, row) in zip(bars_lr, light_rate.iterrows()):
    axes[1].text(bar.get_width() + 1, bar.get_y() + bar.get_height()/2,
                  f'{row["pct"]:.1f}% (n={int(row["count"])})',
                  ha='left', va='center', fontsize=10, fontweight='bold')
axes[1].axvline(avg_rate, color='gray', ls='--', lw=1.5, alpha=0.7)
axes[1].set_title('Severe Rate by Light Condition', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Severe Rate (%)'); axes[1].set_xlim(0, light_rate['pct'].max() * 1.4)

plt.tight_layout(rect=[0, 0, 1, 0.92])
pdf.savefig(fig); plt.close()
print("✔ Page 16 저장 완료")

# ============================================================
# PAGE 17: 종합 요약 대시보드
# ============================================================
fig = plt.figure(figsize=(16, 12))
page_header(fig, 'Page 17. EDA Summary Dashboard',
            'Key Findings from the Exploratory Data Analysis')

# 그리드 레이아웃
gs = gridspec.GridSpec(3, 2, figure=fig, hspace=0.5, wspace=0.35,
                        top=0.90, bottom=0.05, left=0.08, right=0.95)

# --- (17-1) 데이터셋 개요 텍스트 ---
ax1 = fig.add_subplot(gs[0, 0])
ax1.axis('off')
overview_text = (
    f"Dataset Overview\n"
    f"{'='*40}\n"
    f"Total Records: {len(df)}\n"
    f"Variables: {df.shape[1]} ({len(num_cols)} numeric, {len(cat_cols)} categorical)\n"
    f"Target: Accident_Severity (Binary)\n"
    f"  - Non-Severe (0): {(df['Accident_Severity']==0).sum()} "
    f"({(df['Accident_Severity']==0).mean()*100:.1f}%)\n"
    f"  - Severe (1): {(df['Accident_Severity']==1).sum()} "
    f"({(df['Accident_Severity']==1).mean()*100:.1f}%)\n"
    f"Missing Values: 0 (cleaned)\n"
    f"Outliers: Treated (Speed_Limit, Driver_Age)"
)
ax1.text(0.05, 0.95, overview_text, transform=ax1.transAxes,
         fontsize=10, verticalalignment='top', fontfamily='monospace',
         bbox=dict(boxstyle='round,pad=0.5', facecolor='#eaf2f8', alpha=0.8))
ax1.set_title('Dataset Overview', fontsize=12, fontweight='bold')

# --- (17-2) 주요 발견 사항 ---
ax2 = fig.add_subplot(gs[0, 1])
ax2.axis('off')

# 상위 4개 변수와 타겟 연관 강도 계산
num_assoc_sum = {c: abs(df[c].corr(df['Accident_Severity'])) for c in num_cols if c != 'Accident_Severity'}
cat_assoc_sum = {}
for col in cat_cols:
    ct = pd.crosstab(df[col], df['Accident_Severity'])
    chi2, pv, dof, _ = stats.chi2_contingency(ct)
    n = ct.sum().sum(); md = min(ct.shape) - 1
    cat_assoc_sum[col] = np.sqrt(chi2 / (n * md))
all_assoc = {**num_assoc_sum, **cat_assoc_sum}
top4 = sorted(all_assoc.items(), key=lambda x: x[1], reverse=True)[:4]

findings_text = (
    f"Key Findings\n"
    f"{'='*40}\n"
    f"[Top 4 Variables by Association]\n"
    f"  1. {top4[0][0]}: {top4[0][1]:.4f}\n"
    f"  2. {top4[1][0]}: {top4[1][1]:.4f}\n"
    f"  3. {top4[2][0]}: {top4[2][1]:.4f}\n"
    f"  4. {top4[3][0]}: {top4[3][1]:.4f}\n\n"
    f"[Multicollinearity]\n"
    f"  Driver_Age x Experience: r={df['Driver_Age'].corr(df['Driver_Experience']):.4f}\n\n"
    f"[Class Imbalance]\n"
    f"  Severe ratio = {df['Accident_Severity'].mean()*100:.1f}%"
)
ax2.text(0.05, 0.95, findings_text, transform=ax2.transAxes,
         fontsize=10, verticalalignment='top', fontfamily='monospace',
         bbox=dict(boxstyle='round,pad=0.5', facecolor='#fef9e7', alpha=0.8))
ax2.set_title('Key Findings', fontsize=12, fontweight='bold')

# --- (17-3) 상위 5 위험 조합 차트 ---
ax3 = fig.add_subplot(gs[1, :])
combo_top5 = combo[combo['total'] >= 5].sort_values('rate_pct', ascending=False).head(5)
combo_top5 = combo_top5.sort_values('rate_pct', ascending=True)
combo_top5['short_label'] = (combo_top5['Weather'] + ' + ' +
                              combo_top5['Road_Condition'] + ' (TD=' +
                              combo_top5['Traffic_Density'].astype(str) + ')')

cmap5 = plt.cm.Reds(np.linspace(0.4, 0.9, len(combo_top5)))
bars5 = ax3.barh(range(len(combo_top5)), combo_top5['rate_pct'],
                   color=cmap5, edgecolor='black', lw=0.5, height=0.5)
for bar, (_, row) in zip(bars5, combo_top5.iterrows()):
    ax3.text(bar.get_width() + 1, bar.get_y() + bar.get_height()/2,
              f'{row["rate_pct"]:.1f}% (n={int(row["total"])})',
              ha='left', va='center', fontsize=10, fontweight='bold')
ax3.set_yticks(range(len(combo_top5)))
ax3.set_yticklabels(combo_top5['short_label'], fontsize=10)
ax3.set_title('Top 5 Most Dangerous Combinations', fontsize=12, fontweight='bold')
ax3.set_xlabel('Severe Accident Rate (%)')
ax3.axvline(avg_rate, color='gray', ls='--', lw=1.5, alpha=0.7)
ax3.set_xlim(0, 110)

# --- (17-4) 위험 점수 요약 차트 ---
ax4 = fig.add_subplot(gs[2, 0])
risk_bars = ax4.bar(risk_stats['Risk_Score'], risk_stats['rate_pct'],
                     color=[colors_gradient[i] for i in risk_stats['Risk_Score']],
                     edgecolor='black', lw=0.8, width=0.6)
for bar, (_, row) in zip(risk_bars, risk_stats.iterrows()):
    ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
              f'{row["rate_pct"]}%', ha='center', fontsize=9, fontweight='bold')
ax4.set_title('Risk Score Summary', fontsize=12, fontweight='bold')
ax4.set_xlabel('Risk Score (0-4)'); ax4.set_ylabel('Severe Rate (%)')
ax4.set_xticks(range(5))
ax4.set_ylim(0, max(risk_stats['rate_pct']) * 1.2)

# --- (17-5) 분석 결론 텍스트 ---
ax5 = fig.add_subplot(gs[2, 1])
ax5.axis('off')

# 임계점 계산
score_2plus = df[df['Risk_Score'] >= 2]
score_3plus = df[df['Risk_Score'] >= 3]

conclusion_text = (
    f"Conclusions & Recommendations\n"
    f"{'='*40}\n"
    f"1. Road Condition & Weather are the\n"
    f"   strongest predictors of severity.\n\n"
    f"2. Traffic Density amplifies all other\n"
    f"   risk factors significantly.\n\n"
    f"3. Risk Score Thresholds:\n"
    f"   - 2+ factors: {score_2plus['Accident_Severity'].mean()*100:.1f}% severe\n"
    f"   - 3+ factors: {score_3plus['Accident_Severity'].mean()*100:.1f}% severe\n\n"
    f"4. Vehicle Type shows hidden patterns\n"
    f"   under specific conditions (Highway\n"
    f"   + High Density).\n\n"
    f"5. Driver Age & Experience have high\n"
    f"   multicollinearity (r>0.94).\n"
    f"   Consider dropping one variable."
)
ax5.text(0.05, 0.95, conclusion_text, transform=ax5.transAxes,
         fontsize=9.5, verticalalignment='top', fontfamily='monospace',
         bbox=dict(boxstyle='round,pad=0.5', facecolor='#eafaf1', alpha=0.8))
ax5.set_title('Conclusions', fontsize=12, fontweight='bold')

pdf.savefig(fig); plt.close()
print("✔ Page 17 저장 완료")

# ============================================================
# PDF 파일 닫기 — 최종 저장
# ============================================================
pdf.close()
print("\n" + "=" * 60)
print("✅ Full EDA Report 저장 완료!")
print(f"📄 파일 위치: {pdf_path}")
print(f"📊 총 페이지 수: 17페이지")
print("=" * 60)
print("\n[페이지 구성]")
print("  Part 1 (Pages 1-5):  데이터 개요, 기초통계, 분포")
print("  Part 2 (Pages 6-10): 상관관계, 카이제곱, 상호작용")
print("  Part 3 (Pages 11-17): 차량유형, 복합위험, Top10, 운전자분석, 요약")

# 임시 컬럼 정리 (선택)
temp_cols = ['Risk_Weather', 'Risk_Road', 'Risk_Density', 'Risk_Speed',
             'Risk_Score', 'Age_Group', 'Exp_Group']
df.drop(columns=[c for c in temp_cols if c in df.columns], inplace=True)
print("\n🧹 분석용 임시 컬럼 정리 완료")