381 lines
16 KiB
Python
381 lines
16 KiB
Python
from google.colab import files
|
||
import matplotlib.pyplot as plt
|
||
from statsmodels.distributions.empirical_distribution import ECDF
|
||
import math
|
||
from scipy import stats as sp_stats
|
||
import numpy as np
|
||
|
||
|
||
|
||
def load_data(uploaded: dict) -> dict[str, np.ndarray]:
|
||
filename = list(uploaded.keys())[0]
|
||
data = np.loadtxt(filename, delimiter=",", skiprows=1)
|
||
return {"X3": data[:, 2]}
|
||
|
||
def compute_statistics(data: np.ndarray) -> dict:
|
||
return {
|
||
"mean": np.mean(data),
|
||
"var_b": np.var(data),
|
||
"var_ub": np.var(data, ddof=1),
|
||
"std_b": np.std(data),
|
||
"std_ub": np.std(data, ddof=1),
|
||
"median": np.median(data),
|
||
"q25": np.quantile(data, 0.25),
|
||
"q75": np.quantile(data, 0.75),
|
||
"q10": np.quantile(data, 0.10),
|
||
"q90": np.quantile(data, 0.90),
|
||
"min": np.min(data),
|
||
"max": np.max(data),
|
||
"n": len(data),
|
||
}
|
||
|
||
def print_statistics(name: str, stats: dict) -> None:
|
||
print(f"\n{'=' * 40}")
|
||
print(f" Выборка: {name} (n = {stats['n']})")
|
||
print(f"{'=' * 40}")
|
||
print(f" Среднее: {stats['mean']:.4f}")
|
||
print(f" Медиана: {stats['median']:.4f}")
|
||
print(f" Дисперсия смещённая S²: {stats['var_b']:.4f}")
|
||
print(f" Дисперсия несмещённая σ̂²: {stats['var_ub']:.4f}")
|
||
print(f" Ст. откл. смещённое S: {stats['std_b']:.4f}")
|
||
print(f" Ст. откл. несмещённое σ̂: {stats['std_ub']:.4f}")
|
||
print(f" Квартиль Q1 (0.25): {stats['q25']:.4f}")
|
||
print(f" Квартиль Q3 (0.75): {stats['q75']:.4f}")
|
||
print(f" Квантиль 0.10: {stats['q10']:.4f}")
|
||
print(f" Квантиль 0.90: {stats['q90']:.4f}")
|
||
print(f" Минимум: {stats['min']:.4f}")
|
||
print(f" Максимум: {stats['max']:.4f}")
|
||
print(f" IQR: {stats['q75'] - stats['q25']:.4f}")
|
||
|
||
def plot_variation_series(data: np.ndarray, name: str) -> None:
|
||
sorted_data = np.sort(data)
|
||
plt.figure(figsize=(8, 3))
|
||
plt.plot(sorted_data, color="steelblue", linewidth=0.8)
|
||
plt.title(f"Вариационный ряд — {name}", fontweight="bold")
|
||
plt.xlabel("Порядковый номер")
|
||
plt.ylabel("Значение")
|
||
plt.grid(True, alpha=0.3)
|
||
plt.tight_layout()
|
||
|
||
plt.show()
|
||
|
||
def plot_ecdf(data: np.ndarray, name: str) -> None:
|
||
ecdf = ECDF(data)
|
||
x = np.linspace(np.min(data) - 0.5, np.max(data) + 0.5, 1000)
|
||
|
||
plt.figure(figsize=(7, 4))
|
||
plt.step(x, ecdf(x), where="post", color="steelblue", linewidth=1.5)
|
||
plt.title(f"Эмпирическая функция распределения — {name}", fontweight="bold")
|
||
plt.xlabel("x")
|
||
plt.ylabel(r"$F_n(x)$")
|
||
plt.grid(True, alpha=0.3)
|
||
plt.tight_layout()
|
||
|
||
plt.show()
|
||
|
||
def _scott_bins(data: np.ndarray) -> int:
|
||
h = 3.5 * np.std(data, ddof=1) * len(data) ** (-1 / 3)
|
||
return math.ceil((np.max(data) - np.min(data)) / h)
|
||
|
||
|
||
def _freedman_diaconis_bins(data: np.ndarray) -> int:
|
||
iqr = np.quantile(data, 0.75) - np.quantile(data, 0.25)
|
||
h = 2 * iqr * len(data) ** (-1 / 3)
|
||
return math.ceil((np.max(data) - np.min(data)) / h)
|
||
|
||
|
||
def _sturges_bins(data: np.ndarray) -> int:
|
||
return 1 + math.floor(math.log2(len(data)))
|
||
|
||
|
||
_BIN_RULES: dict[str, callable] = {
|
||
"scott": _scott_bins,
|
||
"freedman-diaconis": _freedman_diaconis_bins,
|
||
"sturges": _sturges_bins,
|
||
}
|
||
|
||
|
||
def get_number_of_bins(data: np.ndarray, method: str = "scott") -> int:
|
||
if method not in _BIN_RULES:
|
||
raise ValueError(
|
||
f"unknown method '{method}'. available methods are: {list(_BIN_RULES.keys())}"
|
||
)
|
||
return _BIN_RULES[method](data)
|
||
|
||
|
||
def plot_histogram(
|
||
data: np.ndarray,
|
||
name: str,
|
||
methods: list[str] | None = None,
|
||
) -> None:
|
||
if methods is None:
|
||
methods = list(_BIN_RULES.keys())
|
||
|
||
fig, axes = plt.subplots(
|
||
1, len(methods), figsize=(5 * len(methods), 4), sharey=False
|
||
)
|
||
if len(methods) == 1:
|
||
axes = [axes]
|
||
|
||
for ax, method in zip(axes, methods):
|
||
bins = get_number_of_bins(data, method)
|
||
print(f"{method} -- {name} -- {bins}")
|
||
ax.hist(data, bins=bins, color="steelblue", edgecolor="white", linewidth=0.5)
|
||
ax.set_title(f"{method}\n(k = {bins})", fontsize=10)
|
||
ax.set_xlabel("Значения")
|
||
ax.set_ylabel("Частота")
|
||
|
||
fig.suptitle(f"Гистограммы — {name}", fontsize=13, fontweight="bold")
|
||
plt.tight_layout()
|
||
|
||
plt.show()
|
||
|
||
def describe_distribution(data: np.ndarray, name: str) -> None:
|
||
n = len(data)
|
||
mean = np.mean(data)
|
||
median = np.median(data)
|
||
std = np.std(data, ddof=1)
|
||
skewness = sp_stats.skew(data)
|
||
kurtosis = sp_stats.kurtosis(data)
|
||
q25, q75 = np.quantile(data, [0.25, 0.75])
|
||
iqr = q75 - q25
|
||
|
||
if abs(skewness) < 0.2:
|
||
sym = "симметричное"
|
||
elif skewness > 0:
|
||
sym = f"правосторонняя асимметрия (skew = {skewness:.3f})"
|
||
else:
|
||
sym = f"левосторонняя асимметрия (skew = {skewness:.3f})"
|
||
|
||
if kurtosis > 1:
|
||
tails = f"тяжёлые хвосты (kurt = {kurtosis:.3f})"
|
||
elif kurtosis < -1:
|
||
tails = f"лёгкие хвосты (kurt = {kurtosis:.3f})"
|
||
else:
|
||
tails = f"хвосты близки к нормальным (kurt = {kurtosis:.3f})"
|
||
|
||
lo, hi = q25 - 1.5 * iqr, q75 + 1.5 * iqr
|
||
outliers = data[(data < lo) | (data > hi)]
|
||
if len(outliers) == 0:
|
||
out_str = "выбросы не обнаружены"
|
||
else:
|
||
pct = 100 * len(outliers) / n
|
||
out_str = f"{len(outliers)} выброс(ов) ({pct:.1f}%) за пределами [{lo:.3f}, {hi:.3f}]"
|
||
|
||
kde = sp_stats.gaussian_kde(data)
|
||
x_grid = np.linspace(np.min(data), np.max(data), 1000)
|
||
kde_vals = kde(x_grid)
|
||
peaks = [i for i in range(1, len(kde_vals) - 1)
|
||
if kde_vals[i] > kde_vals[i-1] and kde_vals[i] > kde_vals[i+1]]
|
||
threshold = 0.1 * np.max(kde_vals)
|
||
significant_peaks = [x_grid[i] for i in peaks if kde_vals[i] > threshold]
|
||
|
||
if len(significant_peaks) == 1:
|
||
modes_str = f"унимодальное (пик ≈ {significant_peaks[0]:.3f})"
|
||
elif len(significant_peaks) == 2:
|
||
modes_str = f"бимодальное (пики ≈ {significant_peaks[0]:.3f} и {significant_peaks[1]:.3f})"
|
||
else:
|
||
modes_str = f"многомодальное ({len(significant_peaks)} пика)"
|
||
|
||
print(f"\n{'=' * 40}")
|
||
print(f" Описание формы распределения — {name}")
|
||
print(f"{'=' * 40}")
|
||
print(f" Симметрия: {sym}")
|
||
print(f" Хвосты: {tails}")
|
||
print(f" Выбросы: {out_str}")
|
||
print(f" Модальность: {modes_str}")
|
||
print(f" Среднее vs медиана: {mean:.4f} vs {median:.4f}")
|
||
|
||
|
||
def estimate_exponential_parameters(data: np.ndarray, name: str) -> None:
|
||
n = len(data)
|
||
mean = np.mean(data)
|
||
var_ub = np.var(data, ddof=1)
|
||
|
||
lambda_mm = 1.0 / mean
|
||
|
||
lambda_mm2 = 1.0 / np.sqrt(var_ub)
|
||
|
||
lambda_mle = n / np.sum(data)
|
||
|
||
def log_likelihood(lam: float) -> float:
|
||
return n * np.log(lam) - lam * np.sum(data)
|
||
|
||
ll_mm = log_likelihood(lambda_mm)
|
||
ll_mle = log_likelihood(lambda_mle)
|
||
|
||
print(f"\n{'=' * 50}")
|
||
print(f" Оценка параметров — {name}")
|
||
print(f" Модель: Экспоненциальное распределение Exp(λ)")
|
||
print(f"{'=' * 50}")
|
||
|
||
print(f"\n [Метод моментов]")
|
||
print(f" EX = {lambda_mm:.6f}")
|
||
print(f" DX = {lambda_mm2:.6f} (через дисперсию)")
|
||
|
||
print(f"\n [Метод максимального правдоподобия]")
|
||
print(f" ММП = {lambda_mle:.6f}")
|
||
|
||
print(f"\n [Сравнение]")
|
||
print(f" ММ = {lambda_mm:.6f}")
|
||
print(f" ММП = {lambda_mle:.6f}")
|
||
print(f" Разница = {abs(lambda_mle - lambda_mm):.2e}")
|
||
|
||
if abs(lambda_mle - lambda_mm) < 1e-10:
|
||
comment = ()
|
||
else:
|
||
comment = ("Небольшое расхождение обусловлено численными погрешностями.")
|
||
print(f"\n Оценки совпадают: для Exp(λ) метод моментов (через EX) и МПП дают одинаковую формулу λ̂ = 1/x̄.")
|
||
|
||
lambdas = np.linspace(lambda_mle * 0.5, lambda_mle * 1.8, 500)
|
||
ll_vals = [log_likelihood(l) for l in lambdas]
|
||
|
||
plt.figure(figsize=(7, 4))
|
||
plt.plot(lambdas, ll_vals, color="steelblue", linewidth=1.8,
|
||
label="log L(λ)")
|
||
plt.axvline(lambda_mle, color="crimson", linestyle="--", linewidth=1.2,
|
||
label=f"λ̂_MLE = {lambda_mle:.4f}")
|
||
plt.axvline(lambda_mm, color="darkorange", linestyle=":", linewidth=1.2,
|
||
label=f"λ̂_MM = {lambda_mm:.4f}")
|
||
plt.scatter([lambda_mle, lambda_mm],
|
||
[log_likelihood(lambda_mle), log_likelihood(lambda_mm)],
|
||
color=["crimson", "darkorange"], zorder=5)
|
||
plt.title(f"Логарифм правдоподобия — {name}", fontweight="bold")
|
||
plt.xlabel("λ")
|
||
plt.ylabel("log L(λ)")
|
||
plt.legend()
|
||
plt.grid(True, alpha=0.3)
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
return lambda_mle
|
||
|
||
|
||
def estimate_probability(data: np.ndarray, name: str, lambda_mle: float) -> None:
|
||
n = len(data)
|
||
mean = np.mean(data)
|
||
std_ub = np.std(data, ddof=1)
|
||
c_hat = np.min(data)
|
||
|
||
x0 = mean + std_ub
|
||
|
||
p_empirical = np.sum(data > x0) / n
|
||
|
||
p_parametric = np.exp(-lambda_mle * (x0 - c_hat))
|
||
|
||
print(f"\n{'=' * 50}")
|
||
print(f" Оценка P(X > x0) — {name}")
|
||
print(f" Модель: Exp(λ={lambda_mle:.4f}, c={c_hat:.4f})")
|
||
print(f"{'=' * 50}")
|
||
print(f" Порог x0 = x̄ + σ̂ = {mean:.4f} + {std_ub:.4f} = {x0:.4f}")
|
||
print(f"\n Эмпирическая P(X > x0) = {p_empirical:.6f} ({int(p_empirical*n)}/{n} наблюдений)")
|
||
print(f" Параметрическая P(X > x0) = {p_parametric:.6f}")
|
||
print(f"\n Абсолютное расхождение: {abs(p_parametric - p_empirical):.2e}")
|
||
rel = abs(p_parametric - p_empirical) / p_parametric * 100
|
||
print(f" Относительное расхождение: {rel:.2f}%")
|
||
|
||
if rel < 5:
|
||
comment = "Расхождение незначительное — модель хорошо описывает данные."
|
||
elif rel < 15:
|
||
comment = "Умеренное расхождение — модель в целом подходит, но есть отклонения."
|
||
else:
|
||
comment = "Заметное расхождение — стоит перепроверить выбор модели или оценки параметров."
|
||
print(f" Вывод: {comment}")
|
||
|
||
x_grid = np.linspace(np.min(data) - 0.5, np.max(data) + 2, 500)
|
||
|
||
cdf_param = np.where(x_grid >= c_hat, 1 - np.exp(-lambda_mle * (x_grid - c_hat)), 0.0)
|
||
|
||
from statsmodels.distributions.empirical_distribution import ECDF
|
||
ecdf = ECDF(data)
|
||
|
||
fig, ax = plt.subplots(figsize=(8, 4))
|
||
ax.plot(x_grid, 1 - ecdf(x_grid), color="steelblue", linewidth=1.5, label="Эмпирическая P(X > x)")
|
||
ax.plot(x_grid, 1 - cdf_param, color="crimson", linewidth=1.5, linestyle="--", label=f"Параметрическая P(X > x)")
|
||
ax.axvline(x0, color="darkorange", linestyle=":", linewidth=1.5, label=f"x₀ = {x0:.3f}")
|
||
ax.scatter([x0], [p_empirical], color="steelblue", zorder=5, s=60)
|
||
ax.scatter([x0], [p_parametric], color="crimson", zorder=5, s=60)
|
||
ax.annotate(f" эмп. = {p_empirical:.4f}", (x0, p_empirical), fontsize=9, color="steelblue")
|
||
ax.annotate(f" пар. = {p_parametric:.4f}", (x0, p_parametric), fontsize=9, color="crimson",
|
||
xytext=(x0, p_parametric - 0.03))
|
||
ax.set_title(f"P(X > x) — эмпирическая vs параметрическая — {name}", fontweight="bold")
|
||
ax.set_xlabel("x")
|
||
ax.set_ylabel("P(X > x)")
|
||
ax.legend()
|
||
ax.grid(True, alpha=0.3)
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
def estimate_grouped_moments(data: np.ndarray, name: str, method: str = "scott") -> None:
|
||
n = len(data)
|
||
bins = get_number_of_bins(data, method)
|
||
|
||
counts, edges = np.histogram(data, bins=bins)
|
||
midpoints = (edges[:-1] + edges[1:]) / 2
|
||
|
||
mean_grouped = np.sum(counts * midpoints) / n
|
||
var_grouped = np.sum(counts * (midpoints - mean_grouped) ** 2) / (n - 1)
|
||
std_grouped = np.sqrt(var_grouped)
|
||
|
||
mean_raw = np.mean(data)
|
||
var_raw = np.var(data, ddof=1)
|
||
std_raw = np.std(data, ddof=1)
|
||
|
||
print(f"\n{'=' * 55}")
|
||
print(f" Моменты по сгруппированной выборке — {name}")
|
||
print(f" Правило разбиения: {method} (k = {bins} интервалов)")
|
||
print(f"{'=' * 55}")
|
||
print(f" {'Характеристика':<28} {'По исходным':>12} {'По группир.':>12} {'Δ':>10}")
|
||
print(f" {'-'*62}")
|
||
print(f" {'Среднее EX':<28} {mean_raw:>12.6f} {mean_grouped:>12.6f} {abs(mean_grouped - mean_raw):>10.2e}")
|
||
print(f" {'Дисперсия DX (несмещ.)':<28} {var_raw:>12.6f} {var_grouped:>12.6f} {abs(var_grouped - var_raw):>10.2e}")
|
||
print(f" {'Ст. отклонение σ̂':<28} {std_raw:>12.6f} {std_grouped:>12.6f} {abs(std_grouped - std_raw):>10.2e}")
|
||
|
||
rel_mean = abs(mean_grouped - mean_raw) / abs(mean_raw) * 100
|
||
rel_var = abs(var_grouped - var_raw) / abs(var_raw) * 100
|
||
print(f"\n Относительное отклонение среднего: {rel_mean:.3f}%")
|
||
print(f" Относительное отклонение дисперсии: {rel_var:.3f}%")
|
||
|
||
if max(rel_mean, rel_var) < 1:
|
||
comment = "Группировка практически не вносит погрешности — интервалов достаточно."
|
||
elif max(rel_mean, rel_var) < 5:
|
||
comment = "Незначительная потеря точности из-за группировки — результат приемлем."
|
||
else:
|
||
comment = "Заметная погрешность группировки — стоит увеличить число интервалов."
|
||
print(f" Вывод: {comment}")
|
||
|
||
fig, ax = plt.subplots(figsize=(8, 4))
|
||
ax.bar(edges[:-1], counts, width=np.diff(edges), align="edge",
|
||
color="steelblue", edgecolor="white", linewidth=0.5,
|
||
alpha=0.7, label="Частоты интервалов")
|
||
ax.scatter(midpoints, counts, color="darkorange", zorder=5, s=50, label="Середины интервалов")
|
||
ax.axvline(mean_raw, color="crimson", linestyle="--", linewidth=1.5,
|
||
label=f"x̄ исходное = {mean_raw:.4f}")
|
||
ax.axvline(mean_grouped, color="darkgreen", linestyle=":", linewidth=1.5,
|
||
label=f"x̄ группир. = {mean_grouped:.4f}")
|
||
ax.set_title(f"Гистограмма с оценками среднего — {name}", fontweight="bold")
|
||
ax.set_xlabel("Значения")
|
||
ax.set_ylabel("Частота")
|
||
ax.legend(fontsize=9)
|
||
ax.grid(True, alpha=0.3)
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
|
||
def analyze_column(name: str, data: np.ndarray) -> None:
|
||
stats = compute_statistics(data)
|
||
print_statistics(name, stats)
|
||
describe_distribution(data, name)
|
||
plot_variation_series(data, name)
|
||
plot_ecdf(data, name)
|
||
plot_histogram(data, name)
|
||
|
||
lambda_mle = estimate_exponential_parameters(data, name)
|
||
estimate_probability(data, name, lambda_mle)
|
||
estimate_grouped_moments(data, name)
|
||
|
||
|
||
uploaded = files.upload()
|
||
columns = load_data(uploaded)
|
||
analyze_column("X3", columns["X3"])
|