5 P1 프로젝트: 자전거 수요

import sys
import sklearn
from packaging.version import Version

# 파이썬 버전과 scikit-learn 버전 확인
assert sys.version_info >= (3, 10)
assert Version(sklearn.__version__) >= Version("1.6.1")

# 필요한 라이브러리 임포트
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 그래프 스타일 설정
plt.rc('font', size=12)
plt.rc('axes', labelsize=14, titlesize=14)
plt.rc('legend', fontsize=12)
plt.rc('xtick', labelsize=10)
plt.rc('ytick', labelsize=10)

# UCI Bike Sharing 데이터 불러오기
# 처음 실행할 때 인터넷 연결이 필요하다.

from ucimlrepo import fetch_ucirepo

bike = fetch_ucirepo(id=275)

parts = []

if bike.data.features is not None:
    parts.append(bike.data.features.copy())

if bike.data.targets is not None:
    parts.append(bike.data.targets.copy())

if getattr(bike.data, "ids", None) is not None:
    parts.append(bike.data.ids.copy())

df = pd.concat(parts, axis=1)
df = df.loc[:, ~df.columns.duplicated()]

# 날짜는 이후 시간 순서 정렬에 사용한다.
df["dteday"] = pd.to_datetime(df["dteday"])

print("데이터 크기:", df.shape)
df.head()

required_cols = {
    "instant", "dteday", "season", "yr", "mnth", "hr",
    "holiday", "weekday", "workingday", "weathersit",
    "temp", "atemp", "hum", "windspeed", "cnt"
}

missing_cols = sorted(required_cols - set(df.columns))

if missing_cols:
    raise ValueError(f"필요한 열이 없습니다: {missing_cols}")

print("기본 열 확인:")
print(df.columns.tolist())

df.info()

df.isna().sum().sort_values(ascending=False).head(20)

# UCI 변수 설명에서 대여량 관련 변수를 확인한다.
bike.variables.loc[
    bike.variables["name"].isin(["casual", "registered", "cnt"]),
    ["name", "role", "type", "description"]
]

# 현재 불러온 데이터에 세 열이 실제로 있는지 확인한다.
rental_cols = [c for c in ["casual", "registered", "cnt"] if c in df.columns]

print("현재 데이터에 포함된 대여량 관련 열:", rental_cols)

if {"casual", "registered", "cnt"}.issubset(df.columns):
    display(df[["casual", "registered", "cnt"]].head())
    print(
        "cnt = casual + registered가 모든 행에서 성립:",
        (df["cnt"] == df["casual"] + df["registered"]).all()
    )
else:
    print("casual/registered는 현재 분석 데이터에 포함되어 있지 않습니다.")
    print("다른 방식으로 원자료를 불러올 경우에는 두 열을 입력 특성에서 반드시 제외합니다.")

df_2011 = df.loc[df["yr"] == 0].copy()
df_2012 = df.loc[df["yr"] == 1].copy()

print("2011년 데이터:", df_2011.shape)
print("2012년 데이터:", df_2012.shape)

print(
    "2011년 기간:",
    df_2011["dteday"].min().date(),
    "~",
    df_2011["dteday"].max().date()
)
print(
    "2012년 기간:",
    df_2012["dteday"].min().date(),
    "~",
    df_2012["dteday"].max().date()
)

df_2011["cnt"].hist(bins=30)

plt.xlabel("Hourly bike rentals")
plt.ylabel("Frequency")
plt.title("Distribution of hourly bike rentals: 2011")
plt.show()

hourly = df_2011.groupby("hr")["cnt"].mean()

hourly.plot(marker="o")

plt.xlabel("Hour")
plt.ylabel("Mean rentals")
plt.title("Mean bike rentals by hour: 2011")
plt.show()

hour_work = (
    df_2011.groupby(["hr", "workingday"])["cnt"]
           .mean()
           .unstack()
)

hour_work.columns = ["Non-working day", "Working day"]

hour_work.plot()

plt.xlabel("Hour")
plt.ylabel("Mean rentals")
plt.title("Hourly rentals by working-day status: 2011")
plt.show()

# TODO
# 위에서 정한 EDA 질문을 확인하는 코드를 작성한다.
#

target = "cnt"

exclude_cols = [
    "instant", "dteday", "yr", target,
    "casual", "registered"
]

feature_cols = [
    c for c in df.columns
    if c not in exclude_cols
]

categorical_cols = [
    c for c in [
        "season", "mnth", "hr", "holiday",
        "weekday", "workingday", "weathersit"
    ]
    if c in feature_cols
]

numeric_cols = [
    c for c in feature_cols
    if c not in categorical_cols
]

print("타깃:", target)
print("사용 특성:", feature_cols)
print("범주형:", categorical_cols)
print("수치형:", numeric_cols)

df_2011_ordered = (
    df_2011
    .sort_values(["dteday", "hr"])
    .reset_index(drop=True)
)

split_idx = int(len(df_2011_ordered) * 0.8)

dev_df = df_2011_ordered.iloc[:split_idx].copy()
val_df = df_2011_ordered.iloc[split_idx:].copy()

X_dev = dev_df[feature_cols]
y_dev = dev_df[target]

X_val = val_df[feature_cols]
y_val = val_df[target]

print("모델 학습용:", X_dev.shape)
print(
    "기간:",
    dev_df["dteday"].min().date(),
    "~",
    dev_df["dteday"].max().date()
)

print()

print("모델 선택용 검증:", X_val.shape)
print(
    "기간:",
    val_df["dteday"].min().date(),
    "~",
    val_df["dteday"].max().date()
)

from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder, StandardScaler

numeric_pipe = Pipeline([
    ("imputer", SimpleImputer(strategy="median")),
    ("scaler", StandardScaler())
])

categorical_pipe = Pipeline([
    ("imputer", SimpleImputer(strategy="most_frequent")),
    ("onehot", OneHotEncoder(handle_unknown="ignore"))
])

preprocess = ColumnTransformer([
    ("num", numeric_pipe, numeric_cols),
    ("cat", categorical_pipe, categorical_cols)
])

from sklearn.base import clone
from sklearn.dummy import DummyRegressor
from sklearn.linear_model import LinearRegression, Ridge

def make_pipeline(estimator):
    return Pipeline([
        ("prep", clone(preprocess)),
        ("model", estimator)
    ])

models = {
    "Baseline(mean)": make_pipeline(
        DummyRegressor(strategy="mean")
    ),
    "LinearRegression": make_pipeline(
        LinearRegression()
    )
}

alpha_candidates = [0.01, 0.1, 1, 10, 100]

for alpha in alpha_candidates:
    models[f"Ridge(alpha={alpha})"] = make_pipeline(
        Ridge(alpha=alpha)
    )

print("비교할 모델:")
for name in models:
    print("-", name)

from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    r2_score
)

def regression_scores(model, X, y):
    pred = model.predict(X)

    return {
        "MAE": mean_absolute_error(y, pred),
        "RMSE": np.sqrt(mean_squared_error(y, pred)),
        "R2": r2_score(y, pred)
    }

validation_rows = []

for name, model in models.items():
    model.fit(X_dev, y_dev)

    train_scores = regression_scores(model, X_dev, y_dev)
    val_scores = regression_scores(model, X_val, y_val)

    validation_rows.append({
        "Model": name,
        "Train MAE": train_scores["MAE"],
        "Validation MAE": val_scores["MAE"],
        "Train RMSE": train_scores["RMSE"],
        "Validation RMSE": val_scores["RMSE"],
        "Train R2": train_scores["R2"],
        "Validation R2": val_scores["R2"]
    })

validation_results = (
    pd.DataFrame(validation_rows)
      .set_index("Model")
      .sort_values("Validation RMSE")
)

validation_results["RMSE 차이(검증-훈련)"] = (
    validation_results["Validation RMSE"]
    - validation_results["Train RMSE"]
)

validation_results

# 검증 RMSE만 기준으로 하면 가장 낮은 모델은 무엇인지 확인한다.
suggested_model_name = validation_results["Validation RMSE"].idxmin()

print("검증 RMSE 기준 후보:", suggested_model_name)

chosen_model_name = suggested_model_name

# 다른 모델을 선택할 근거가 있다면 위 값을 직접 수정한다.
# 예: chosen_model_name = "Ridge(alpha=10)"

assert chosen_model_name in models

print("최종 선택 모델:", chosen_model_name)

X_2011 = df_2011[feature_cols]
y_2011 = df_2011[target]

X_2012 = df_2012[feature_cols]
y_2012 = df_2012[target]

final_model = clone(models[chosen_model_name])
final_model.fit(X_2011, y_2011)

train_2011_scores = regression_scores(
    final_model, X_2011, y_2011
)

test_2012_scores = regression_scores(
    final_model, X_2012, y_2012
)

final_results = pd.DataFrame({
    "2011 train": train_2011_scores,
    "2012 final test": test_2012_scores
})

final_results

pred_2012 = final_model.predict(X_2012)

plt.scatter(y_2012, pred_2012, alpha=0.15)

low = min(y_2012.min(), pred_2012.min())
high = max(y_2012.max(), pred_2012.max())

plt.plot([low, high], [low, high])

plt.xlabel("Actual rentals in 2012")
plt.ylabel("Predicted rentals")
plt.title("Actual vs predicted bike rentals: 2012")
plt.show()

errors_2012 = pd.DataFrame({
    "hr": X_2012["hr"].to_numpy(),
    "workingday": X_2012["workingday"].to_numpy(),
    "weathersit": X_2012["weathersit"].to_numpy(),
    "actual": y_2012.to_numpy(),
    "predicted": pred_2012
})

errors_2012["error"] = (
    errors_2012["actual"]
    - errors_2012["predicted"]
)

errors_2012["abs_error"] = (
    errors_2012["error"].abs()
)

hour_error = (
    errors_2012
    .groupby("hr")["abs_error"]
    .mean()
)

hour_error.plot(marker="o")

plt.xlabel("Hour")
plt.ylabel("Mean absolute error")
plt.title("Prediction error by hour: 2012")
plt.show()

analysis_col = "workingday"  # TODO: 필요하면 "weathersit"으로 바꾼다.

condition_error = (
    errors_2012
    .groupby(analysis_col)["abs_error"]
    .agg(["mean", "median", "count"])
)

condition_error

df.groupby("yr")["cnt"].agg(
    ["mean", "median", "std", "min", "max"]
)

monthly = (
    df.groupby(["yr", "mnth"])["cnt"]
      .mean()
      .unstack(0)
)

monthly.columns = ["2011", "2012"]

monthly.plot(marker="o")

plt.xlabel("Month")
plt.ylabel("Mean hourly rentals")
plt.title("Monthly bike demand: 2011 vs 2012")
plt.show()