Chapter_2.py - master - vasilenko/Chapter_1

Chapter_2.py
187 строк9 KB
vasilenko
Create: Chapter_2.py
11 июн 2026, 15:59
11 июн 2026, 15:59c57c29e
# -*- coding: utf-8 -*-
"""
Раздел 2 курсового проекта: первичный анализ временных рядов.
Набор данных: Household Electric Power Consumption (UCI, id 235).
Скрипт воспроизводит этапы 1-7 учебно-методического материала.
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.tsa.seasonal import seasonal_decompose

# Путь к файлу с данными (скачать: https://archive.ics.uci.edu/dataset/235)
DATA_PATH = "ryadi.txt"

# Русские названия каналов — для подписей на графиках и в таблицах
NAMES = {
    "Global_active_power": "Активная мощность",
    "Global_reactive_power": "Реактивная мощность",
    "Voltage": "Напряжение",
    "Global_intensity": "Сила тока",
    "Sub_metering_1": "Счётчик 1 (кухня)",
    "Sub_metering_2": "Счётчик 2 (прачечная)",
    "Sub_metering_3": "Счётчик 3 (ВО/конд.)",
}

# ============================================
# ЧАСТЬ 1 — Загрузка временного ряда (этап 1)
# ============================================
# В файле разделитель «;», а пропуски записаны знаком «?».
df = pd.read_csv(DATA_PATH, sep=";", na_values="?", low_memory=False)

# Объединяем столбцы даты и времени в один индекс типа datetime
df["datetime"] = pd.to_datetime(df["Date"] + " " + df["Time"],
                                format="%d/%m/%Y %H:%M:%S")
df = df.drop(columns=["Date", "Time"]).set_index("datetime").sort_index()

# Все каналы должны быть числовыми
df = df.astype(float)

print("=== Этап 1. Загрузка и первичное знакомство ===")
print("Исходные данные:", df.shape[0], "записей,", df.shape[1], "каналов")
print(df.head(), "\n")
print(df.info(), "\n")

# ============================================
# ЧАСТЬ 2 — Предобработка данных
# ============================================
# Минутных данных очень много (более 2 млн строк), поэтому
# агрегируем (усредняем) их до одного значения в час.
df = df.resample("1h").mean()

# После усреднения единичные пропуски исчезают,
# оставшиеся неполные строки удаляем.
df = df.dropna()
print("После агрегации до часа:", df.shape[0], "наблюдений\n")

# ============================================
# ЧАСТЬ 3 — Анализ ряда (этапы 2-7)
# ============================================

# --- Этап 2. Визуализация с границей train/test ---
split = int(len(df) * 0.8)  # 80% — обучение, 20% — тест
border = df.index[split]

fig, axes = plt.subplots(len(df.columns), 1, figsize=(12, 14), sharex=True)
for ax, col in zip(axes, df.columns):
    ax.plot(df.index, df[col], linewidth=0.4)
    ax.axvline(border, color="red", linestyle="--", label="граница train/test")
    ax.set_ylabel(NAMES[col], fontsize=8)
    ax.grid(True, alpha=0.3)
axes[0].legend(loc="upper right")
fig.suptitle("Каналы временного ряда с разделением train/test")
plt.tight_layout()
plt.savefig("ris25_kanaly.png", dpi=150)

# --- Этап 3. Статистический анализ ---
print("=== Этап 3. Описательные статистики ===")
desc = df.describe().T.rename(index=NAMES)
print(desc.round(4), "\n")

# Частота дискретизации = типичный шаг между соседними отметками времени
step = df.index.to_series().diff().mode()[0]
print("Частота дискретизации:", step, "\n")

# --- Этап 4. Пропуски и выбросы ---
print("=== Этап 4. Пропуски и выбросы ===")
missing = df.isna().mean() * 100
print("Доля пропусков, %:")
print(missing.rename(NAMES).round(2), "\n")

# Правило трёх сигм: выброс — значение, отклоняющееся от среднего
# больше чем на три стандартных отклонения.
print("Выбросы по правилу трёх сигм:")
for col in df.columns:
    z = np.abs((df[col] - df[col].mean()) / df[col].std())
    n = int((z > 3).sum())
    print(f"  {NAMES[col]:25s}: {n:5d} шт. ({n / len(df) * 100:.4f} %)")
print()

# Диаграммы размаха по каждому каналу
fig, axes = plt.subplots(2, 4, figsize=(14, 6))
for ax, col in zip(axes.ravel(), df.columns):
    ax.boxplot(df[col].dropna())
    ax.set_title(NAMES[col], fontsize=8)
axes.ravel()[-1].axis("off")
fig.suptitle("Диаграммы размаха (boxplot) по каждому каналу")
plt.tight_layout()
plt.savefig("ris26_boxplot.png", dpi=150)

# --- Этап 5. Сравнение диапазонов значений ---
plt.figure(figsize=(12, 5))
plt.boxplot([df[c].dropna() for c in df.columns])
plt.xticks(range(1, len(df.columns) + 1), [NAMES[c] for c in df.columns],
           rotation=30, ha="right", fontsize=8)
plt.title("Сравнение диапазонов значений всех каналов")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ris27_diapazony.png", dpi=150)

# --- Этап 6. Корреляционный анализ ---
corr = df.corr(method="pearson")
plt.figure(figsize=(8, 6))
im = plt.imshow(corr, cmap="coolwarm", vmin=-1, vmax=1)
plt.colorbar(im)
labels = [NAMES[c] for c in corr.columns]
plt.xticks(range(len(labels)), labels, rotation=45, ha="right", fontsize=8)
plt.yticks(range(len(labels)), labels, fontsize=8)
# подписываем каждую ячейку числом
for i in range(len(corr)):
    for j in range(len(corr)):
        plt.text(j, i, f"{corr.iloc[i, j]:.2f}",
                 ha="center", va="center", fontsize=7)
plt.title("Корреляционная матрица Пирсона")
plt.tight_layout()
plt.savefig("ris28_korrelyaciya.png", dpi=150)
print("=== Этап 6. Максимальная корреляция ===")
print("Активная мощность и сила тока:",
      round(corr.loc["Global_active_power", "Global_intensity"], 4), "\n")

# --- Этап 7. Декомпозиция и анализ шума ---
# Ключевой канал — активная мощность. Аддитивная модель, период 24 часа
# (одни сутки), потому что у данных суточная сезонность.
key = df["Global_active_power"]
dec = seasonal_decompose(key, model="additive", period=24)

fig, axes = plt.subplots(4, 1, figsize=(12, 9), sharex=True)
for ax, part, title in zip(
        axes,
        [key, dec.trend, dec.seasonal, dec.resid],
        ["Исходный ряд", "Тренд", "Сезонность", "Остатки (шум)"]):
    ax.plot(part, linewidth=0.4)
    ax.set_ylabel(title, fontsize=8)
    ax.grid(True, alpha=0.3)
fig.suptitle("Декомпозиция временного ряда (активная мощность)")
plt.tight_layout()
plt.savefig("ris29_dekompoziciya.png", dpi=150)

# Сигнал = тренд + сезонность (аддитивная модель).
signal = (dec.trend + dec.seasonal).dropna()
noise = dec.resid.dropna()

# SNR (отношение сигнал/шум) в децибелах:
# чем выше, тем меньше шум мешает видеть закономерности.
snr = 10 * np.log10(signal.var() / noise.var())
print("=== Этап 7. Анализ шума ===")
print(f"SNR = {snr:.2f} дБ "
      "(меньше 0 — шум сильнее сигнала, оценка «плохо»)")
print(f"Асимметрия остатков (skewness): {stats.skew(noise):.2f}")
print(f"Эксцесс остатков (kurtosis):    {stats.kurtosis(noise):.2f}")

# Гистограмма остатков — чтобы оценить форму распределения шума
plt.figure(figsize=(8, 5))
plt.hist(noise, bins=100, edgecolor="black", linewidth=0.3)
plt.title("Гистограмма распределения остатков (шума)")
plt.xlabel("Значение остатка")
plt.ylabel("Частота")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ris30_ostatki.png", dpi=150)

# ============================================
# ЧАСТЬ 4 — Итоговый вывод
# ============================================
print("\nГотово! Графики сохранены в файлы ris25...ris30 (.png).")
plt.show()