476 lines
12 KiB
Python
476 lines
12 KiB
Python
import pandas as pd
|
|
import numpy as np
|
|
import joblib
|
|
import os
|
|
import matplotlib.pyplot as plt
|
|
from sqlalchemy import create_engine
|
|
from xgboost import XGBRegressor
|
|
from sklearn.tree import DecisionTreeRegressor
|
|
from sklearn.metrics import mean_absolute_error
|
|
from sklearn.metrics import mean_squared_error
|
|
from sklearn.metrics import r2_score
|
|
from sklearn.metrics import mean_absolute_percentage_error
|
|
from sklearn.ensemble import RandomForestRegressor
|
|
from sklearn.ensemble import GradientBoostingRegressor
|
|
from sklearn.linear_model import LinearRegression
|
|
|
|
user = "root"
|
|
password = ""
|
|
host = "localhost"
|
|
db = "udd_pmi_module"
|
|
engine = create_engine(f"mysql+pymysql://{user}:{password}@{host}/{db}")
|
|
|
|
# =========================
|
|
# FOLDER OUTPUT VISUALISASI
|
|
# =========================
|
|
output_folder = "machine_learning"
|
|
os.makedirs(output_folder, exist_ok=True)
|
|
|
|
query = """
|
|
SELECT
|
|
m.kode,
|
|
m.tahun,
|
|
m.bulan,
|
|
m.total_pemakaian
|
|
FROM item_usage_monthly m
|
|
JOIN items i ON m.kode = i.kode
|
|
WHERE i.kategori = 0
|
|
ORDER BY m.kode, m.tahun, m.bulan;
|
|
"""
|
|
|
|
# =========================
|
|
# LOAD DATA
|
|
# =========================
|
|
df = pd.read_sql(query, engine)
|
|
df = df.sort_values(by=["kode", "tahun", "bulan"]).reset_index(drop=True)
|
|
|
|
# =========================
|
|
# VALIDASI DATA
|
|
# =========================
|
|
df = df.drop_duplicates()
|
|
df = df[df["bulan"].between(1, 12)]
|
|
df = df[df["tahun"] > 2000]
|
|
df = df[df["total_pemakaian"].fillna(0) >= 0]
|
|
df = df.dropna(subset=["kode", "tahun", "bulan"])
|
|
|
|
# =========================
|
|
# HANDLE MISSING VALUE
|
|
# =========================
|
|
df["total_pemakaian"] = df["total_pemakaian"].fillna(0)
|
|
|
|
# =========================
|
|
# HANDLE OUTLIER (VALIDASI NILAI MINIMUM)
|
|
# =========================
|
|
df["total_pemakaian"] = df["total_pemakaian"].clip(lower=0)
|
|
|
|
# =========================
|
|
# TYPE CASTING & ROUNDING
|
|
# =========================
|
|
df["total_pemakaian"] = df["total_pemakaian"].round(0).astype(int)
|
|
|
|
# =========================
|
|
# FILTER BARANG STABIL
|
|
# =========================
|
|
# std_per_barang = df.groupby("kode")["total_pemakaian"].std()
|
|
# mean_per_barang = df.groupby("kode")["total_pemakaian"].mean()
|
|
# cv = (std_per_barang / mean_per_barang).fillna(0)
|
|
# barang_stabil = cv[cv < 1.0].index
|
|
# df = df[df["kode"].isin(barang_stabil)]
|
|
|
|
# =========================
|
|
# LAG FEATURES
|
|
# =========================
|
|
df["lag_1"] = df.groupby("kode")["total_pemakaian"].shift(1)
|
|
df["lag_2"] = df.groupby("kode")["total_pemakaian"].shift(2)
|
|
df["lag_3"] = df.groupby("kode")["total_pemakaian"].shift(3)
|
|
df["lag_4"] = df.groupby("kode")["total_pemakaian"].shift(4)
|
|
df["lag_5"] = df.groupby("kode")["total_pemakaian"].shift(5)
|
|
df["lag_6"] = df.groupby("kode")["total_pemakaian"].shift(6)
|
|
df["quarter"] = ((df["bulan"] - 1) // 3) + 1
|
|
df["is_awal_tahun"] = df["bulan"].isin([1, 2, 3]).astype(int)
|
|
df["is_tengah_tahun"] = df["bulan"].isin([6, 7, 8]).astype(int)
|
|
df["is_akhir_tahun"] = df["bulan"].isin([10, 11, 12]).astype(int)
|
|
|
|
# =========================
|
|
# DROP NAN
|
|
# =========================
|
|
df = df[df["lag_6"].notna()]
|
|
|
|
# =========================
|
|
# MOVING AVERAGE
|
|
# =========================
|
|
df["ma_6"] = df[["lag_1", "lag_2", "lag_3", "lag_4", "lag_5", "lag_6"]].mean(axis=1)
|
|
|
|
# =========================
|
|
# TREND FEATURE
|
|
# =========================
|
|
df["trend_3"] = df["lag_1"] - df[["lag_2", "lag_3", "lag_4"]].mean(axis=1)
|
|
|
|
# =========================
|
|
# MOMENTUM FEATURE
|
|
# =========================
|
|
df["momentum_1"] = df["lag_1"] - df["lag_2"]
|
|
df["momentum_2"] = df["lag_2"] - df["lag_3"]
|
|
df["rolling_std_6"] = df[["lag_1", "lag_2", "lag_3", "lag_4", "lag_5", "lag_6"]].std(
|
|
axis=1
|
|
)
|
|
df["max_6"] = df[["lag_1", "lag_2", "lag_3", "lag_4", "lag_5", "lag_6"]].max(axis=1)
|
|
df["min_6"] = df[["lag_1", "lag_2", "lag_3", "lag_4", "lag_5", "lag_6"]].min(axis=1)
|
|
|
|
# =========================
|
|
# FITUR TAMBAHAN
|
|
# =========================
|
|
|
|
df["range_6"] = df["max_6"] - df["min_6"]
|
|
|
|
df["cv_6"] = (
|
|
df["rolling_std_6"] /
|
|
(df["ma_6"] + 1)
|
|
)
|
|
|
|
df["growth_rate"] = (
|
|
(df["lag_1"] - df["lag_2"]) /
|
|
(df["lag_2"] + 1)
|
|
)
|
|
|
|
df["usage_level"] = pd.qcut(
|
|
df["ma_6"],
|
|
q=4,
|
|
labels=False,
|
|
duplicates="drop"
|
|
)
|
|
|
|
# =========================
|
|
# SPLIT TRAINING & TEST (TIME-BASED)
|
|
# =========================
|
|
df["date"] = pd.to_datetime(
|
|
df["tahun"].astype(str) + "-" + df["bulan"].astype(str) + "-01"
|
|
)
|
|
|
|
df = df.sort_values(["date","kode"])
|
|
|
|
split_date = df["date"].quantile(0.8)
|
|
|
|
train_df = df[df["date"] <= split_date]
|
|
test_df = df[df["date"] > split_date]
|
|
|
|
# =========================
|
|
# PREDIKSI
|
|
# =========================
|
|
features = [
|
|
"lag_1",
|
|
"lag_2",
|
|
"lag_3",
|
|
"lag_4",
|
|
"lag_5",
|
|
"lag_6",
|
|
"ma_6",
|
|
"trend_3",
|
|
"momentum_1",
|
|
"momentum_2",
|
|
"rolling_std_6",
|
|
"max_6",
|
|
"min_6",
|
|
|
|
"range_6",
|
|
"cv_6",
|
|
"growth_rate",
|
|
"usage_level",
|
|
|
|
"bulan",
|
|
"quarter",
|
|
"is_awal_tahun",
|
|
"is_tengah_tahun",
|
|
"is_akhir_tahun",
|
|
]
|
|
X_train = train_df[features]
|
|
y_train = np.log1p(train_df["total_pemakaian"])
|
|
X_test = test_df[features]
|
|
y_test = test_df["total_pemakaian"]
|
|
|
|
# =========================
|
|
# MODEL RANDOM FOREST
|
|
# =========================
|
|
rf_model = RandomForestRegressor(
|
|
n_estimators=300,
|
|
max_depth=15,
|
|
min_samples_split=5,
|
|
min_samples_leaf=2,
|
|
max_features="sqrt",
|
|
bootstrap=True,
|
|
random_state=42,
|
|
n_jobs=-1,
|
|
)
|
|
rf_model.fit(X_train, y_train)
|
|
rf_pred_log = rf_model.predict(X_test)
|
|
rf_pred = np.expm1(rf_pred_log)
|
|
rf_pred_real = np.round(rf_pred).astype(int)
|
|
|
|
# =========================
|
|
# MODEL GRADIENT BOOSTING
|
|
# =========================
|
|
gb_model = GradientBoostingRegressor(
|
|
n_estimators=300,
|
|
learning_rate=0.03,
|
|
max_depth=5,
|
|
random_state=42,
|
|
)
|
|
gb_model.fit(X_train, y_train)
|
|
gb_pred_log = gb_model.predict(X_test)
|
|
gb_pred = np.expm1(gb_pred_log)
|
|
gb_pred_real = np.round(gb_pred).astype(int)
|
|
|
|
# =========================
|
|
# MODEL LINEAR REGRESSION
|
|
# =========================
|
|
lr_model = LinearRegression()
|
|
lr_model.fit(X_train, y_train)
|
|
lr_pred_log = lr_model.predict(X_test)
|
|
lr_pred = np.expm1(lr_pred_log)
|
|
lr_pred_real = np.round(lr_pred).astype(int)
|
|
|
|
# =========================
|
|
# MODEL XGBOOST
|
|
# =========================
|
|
xgb_model = XGBRegressor(
|
|
n_estimators=300,
|
|
learning_rate=0.03,
|
|
max_depth=5,
|
|
subsample=0.8,
|
|
colsample_bytree=0.8,
|
|
objective="reg:squarederror",
|
|
random_state=42,
|
|
)
|
|
xgb_model.fit(X_train, y_train)
|
|
xgb_pred_log = xgb_model.predict(X_test)
|
|
xgb_pred = np.expm1(xgb_pred_log)
|
|
xgb_pred_real = np.round(xgb_pred).astype(int)
|
|
|
|
# =========================
|
|
# MODEL DECISION TREE
|
|
# =========================
|
|
dt_model = DecisionTreeRegressor(
|
|
max_depth=10, min_samples_split=12, min_samples_leaf=6, random_state=42
|
|
)
|
|
dt_model.fit(X_train, y_train)
|
|
dt_pred_log = dt_model.predict(X_test)
|
|
dt_pred = np.expm1(dt_pred_log)
|
|
dt_pred_real = np.round(dt_pred).astype(int)
|
|
|
|
# =========================
|
|
# CLIP NILAI MINIMUM
|
|
# =========================
|
|
rf_pred_real = np.clip(rf_pred_real, 0, None)
|
|
gb_pred_real = np.clip(gb_pred_real, 0, None)
|
|
lr_pred_real = np.clip(lr_pred_real, 0, None)
|
|
xgb_pred_real = np.clip(xgb_pred_real, 0, None)
|
|
dt_pred_real = np.clip(dt_pred_real, 0, None)
|
|
|
|
# =========================
|
|
# SAVE LOAD MODEL
|
|
# =========================
|
|
# joblib.dump(rf_model, f"{output_folder}/model.pkl")
|
|
|
|
# =========================
|
|
# FUNCTION EVALUASI
|
|
# =========================
|
|
def smape(y_true,y_pred):
|
|
|
|
y_true=np.array(y_true)
|
|
y_pred=np.array(y_pred)
|
|
|
|
return np.mean(
|
|
(
|
|
np.abs(y_true-y_pred)
|
|
/
|
|
(
|
|
np.abs(y_true)
|
|
+
|
|
np.abs(y_pred)
|
|
+1
|
|
)
|
|
)
|
|
)*100
|
|
|
|
|
|
def evaluate_model(name,y_true,y_pred):
|
|
|
|
mae=mean_absolute_error(
|
|
y_true,
|
|
y_pred
|
|
)
|
|
|
|
mse=mean_squared_error(
|
|
y_true,
|
|
y_pred
|
|
)
|
|
|
|
rmse=np.sqrt(mse)
|
|
|
|
r2=r2_score(
|
|
y_true,
|
|
y_pred
|
|
)
|
|
|
|
mape=mean_absolute_percentage_error(
|
|
y_true.replace(0,1),
|
|
y_pred
|
|
)*100
|
|
|
|
smape_value=smape(
|
|
y_true,
|
|
y_pred
|
|
)
|
|
|
|
print(f"\n=== {name} ===")
|
|
print("MAE :",round(mae,2))
|
|
# print("RMSE :",round(rmse,2))
|
|
print("R2 :",round(r2,2))
|
|
print("MAPE :",round(mape,2),"%")
|
|
print("SMAPE :",round(smape_value,2),"%")
|
|
|
|
# =========================
|
|
# EVALUASI SEMUA MODEL
|
|
# =========================
|
|
evaluate_model("RANDOM FOREST", y_test, rf_pred_real)
|
|
evaluate_model("GRADIENT BOOSTING", y_test, gb_pred_real)
|
|
evaluate_model("LINEAR REGRESSION", y_test, lr_pred_real)
|
|
evaluate_model("XGBOOST", y_test, xgb_pred_real)
|
|
evaluate_model("DECISION TREE", y_test, dt_pred_real)
|
|
|
|
# =========================
|
|
# HASIL PREDIKSI RANDOM FOREST
|
|
# =========================
|
|
limit = min(10, len(rf_pred_real))
|
|
hasil_tabel = pd.DataFrame(
|
|
{
|
|
"Data": [f"Data {i+1}" for i in range(limit)],
|
|
"Actual": y_test.iloc[:limit].values,
|
|
"RandomForest": rf_pred_real[:limit],
|
|
"GradientBoost": gb_pred_real[:limit],
|
|
"LinearRegression": lr_pred_real[:limit],
|
|
"XGBOOST": xgb_pred_real[:limit],
|
|
"DecisionTree": dt_pred_real[:limit],
|
|
}
|
|
)
|
|
|
|
# =========================
|
|
# FUNCTION METRIK VISUAL
|
|
# =========================
|
|
def get_metrics(y_true, y_pred):
|
|
mae = mean_absolute_error(y_true, y_pred)
|
|
r2 = r2_score(y_true, y_pred)
|
|
mape = mean_absolute_percentage_error(y_true.replace(0, 1), y_pred) * 100
|
|
return mae, mape, r2
|
|
|
|
# =========================
|
|
# AMBIL HASIL METRIK
|
|
# =========================
|
|
rf_mae, rf_mape, rf_r2 = get_metrics(y_test, rf_pred_real)
|
|
lr_mae, lr_mape, lr_r2 = get_metrics(y_test, lr_pred_real)
|
|
dt_mae, dt_mape, dt_r2 = get_metrics(y_test, dt_pred_real)
|
|
models = ["Random Forest", "Linear Regression", "Decision Tree"]
|
|
mae_values = [rf_mae, lr_mae, dt_mae]
|
|
mape_values = [rf_mape, lr_mape, dt_mape]
|
|
r2_values = [rf_r2, lr_r2, dt_r2]
|
|
|
|
# =========================
|
|
# VISUALISASI KOMPARASI
|
|
# =========================
|
|
fig, axes = plt.subplots(
|
|
1,
|
|
3,
|
|
figsize=(18, 5.5),
|
|
gridspec_kw={"wspace": 0.35}
|
|
)
|
|
fig.suptitle(
|
|
"Evaluasi dan Komparasi Model Prediksi",
|
|
fontsize=22,
|
|
fontweight="bold",
|
|
y=1.02
|
|
)
|
|
|
|
# =========================
|
|
# STYLE UMUM
|
|
# =========================
|
|
for ax in axes:
|
|
ax.set_facecolor("#f8fafc")
|
|
for spine in ax.spines.values():
|
|
spine.set_edgecolor("#d1d5db")
|
|
spine.set_linewidth(1)
|
|
ax.tick_params(axis="x", labelsize=10, rotation=0, pad=6)
|
|
ax.tick_params(axis="y", labelsize=10)
|
|
|
|
# =========================
|
|
# MAE
|
|
# =========================
|
|
bars = axes[0].bar(models, mae_values, width=0.55)
|
|
axes[0].set_title("Perbandingan MAE", fontsize=15, pad=12)
|
|
axes[0].set_ylabel("Nilai MAE", fontsize=12)
|
|
for bar in bars:
|
|
yval = bar.get_height()
|
|
axes[0].text(
|
|
bar.get_x() + bar.get_width() / 2,
|
|
yval + (max(mae_values) * 0.01),
|
|
round(yval, 2),
|
|
ha="center",
|
|
va="bottom",
|
|
fontsize=10
|
|
)
|
|
|
|
# =========================
|
|
# MAPE
|
|
# =========================
|
|
bars = axes[1].bar(models, mape_values, width=0.55)
|
|
axes[1].set_title("Perbandingan MAPE", fontsize=15, pad=12)
|
|
axes[1].set_ylabel("Persentase (%)", fontsize=12)
|
|
for bar in bars:
|
|
yval = bar.get_height()
|
|
axes[1].text(
|
|
bar.get_x() + bar.get_width() / 2,
|
|
yval + (max(mape_values) * 0.02),
|
|
round(yval, 2),
|
|
ha="center",
|
|
va="bottom",
|
|
fontsize=10
|
|
)
|
|
|
|
# =========================
|
|
# R2
|
|
# =========================
|
|
r2_visual = [max(0, value) for value in r2_values]
|
|
bars = axes[2].bar(models, r2_visual, width=0.55)
|
|
axes[2].set_title("Perbandingan R²", fontsize=15, pad=12)
|
|
axes[2].set_ylabel("R² Score", fontsize=12)
|
|
axes[2].set_ylim(0, max(r2_visual) + 0.1)
|
|
for i, bar in enumerate(bars):
|
|
visual_height = bar.get_height()
|
|
real_value = r2_values[i]
|
|
axes[2].text(
|
|
bar.get_x() + bar.get_width() / 2,
|
|
visual_height + 0.02,
|
|
round(real_value, 2),
|
|
ha="center",
|
|
va="bottom",
|
|
fontsize=10
|
|
)
|
|
|
|
# =========================
|
|
# RAPATKAN & AMANKAN LABEL
|
|
# =========================
|
|
plt.tight_layout(rect=[0.02, 0.06, 1, 0.93])
|
|
|
|
# =========================
|
|
# SIMPAN GAMBAR
|
|
# =========================
|
|
# plt.savefig(
|
|
# f"{output_folder}/evaluasi_komparasi_model.png",
|
|
# dpi=300,
|
|
# bbox_inches="tight"
|
|
# )
|
|
# plt.close()
|
|
# print("\nVisualisasi berhasil dibuat")
|
|
# print(f"Lokasi file: {output_folder}/evaluasi_komparasi_model.png")
|
|
|
|
print("\n=== PERBANDINGAN HASIL PREDIKSI ===\n")
|
|
print(hasil_tabel.to_string(index=False)) |