https://pastein.ru/t/p5a
Загрузка данных
"""
Скелет вычислительного эксперимента к ВКР:
"Бюджетное распределение превентивного контроля по ожидаемому ущербу
от поздних дефектов в проектах мобильной разработки".
Что делает скрипт
-----------------
1. Генерирует полусинтетический «мир» НЕЗАВИСИМЫМ генератором: истинная
вероятность позднего дефекта зависит от скрытого размера задачи, компонента,
опыта исполнителя и новизны через cloglog-связь (а НЕ через логистическую
функцию от Story Points). Story Points — зашумлённая оценка размера.
2. Обучает модель вероятности (логистическая регрессия) только на признаках,
доступных при планировании (SP; опционально — компонент).
3. Сравнивает политики распределения усиленного контроля при ОДНОМ И ТОМ ЖЕ
бюджете B и в ОДНОМ И ТОМ ЖЕ мире (истинные P, C, g), то есть различается
только правило решения.
4. Считает метрики: истинные ожидаемые затраты, предотвращённый ущерб,
recall@budget, регрет к оракулу, Brier/ECE/ROC-AUC/PR-AUC модели.
5. Сохраняет таблицы в CSV. Числа, которые выдаёт скрипт, — результаты ВАШЕГО
запуска при ВАШИХ параметрах; параметры по умолчанию — допущения-заглушки.
Запуск
------
python experiment_skeleton.py --reps 50 --sprints 30 --out results
python experiment_skeleton.py --legacy # проверочный расчёт примера §2.4 ВКР
python experiment_skeleton.py --sweep budget # свип по бюджету
python experiment_skeleton.py --sweep k_mismatch
python experiment_skeleton.py --sweep cost_cv
python experiment_skeleton.py --sweep threshold
python experiment_skeleton.py --sweep miscalibration --budget-fixed 5000
python experiment_skeleton.py --sweep shift
python experiment_skeleton.py --control-check # частный случай = пороговое правило
Зависимости: numpy, scipy, pandas, scikit-learn (Python 3.10+).
"""
from __future__ import annotations
import argparse
import dataclasses
import os
from dataclasses import dataclass, replace
import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import average_precision_score, brier_score_loss, roc_auc_score
FIB = np.array([1, 2, 3, 5, 8, 13, 21])
COMPONENTS = ("ui", "core", "network", "payments", "auth")
# ----------------------------------------------------------------------------
# 1. Конфигурация мира (все значения — допущения, подлежат свипу)
# ----------------------------------------------------------------------------
@dataclass(frozen=True)
class Config:
n_hist: int = 6000 # исторические задачи для обучения модели
test_share: float = 0.30 # доля отложенной выборки для метрик модели
n_sprint: int = 100 # задач в планируемом спринте
base_rate: float = 0.15 # целевая доля поздних дефектов в мире
comp_share: tuple = (0.40, 0.30, 0.15, 0.10, 0.05)
comp_kappa: tuple = (1.0, 2.0, 2.0, 8.0, 6.0) # критичность компонента
comp_prod_share: tuple = (0.25, 0.30, 0.35, 0.40, 0.40) # доля production среди поздних
comp_shift: tuple = (0.0, 0.2, 0.3, 0.4, 0.3) # вклад компонента в P
C_qa: float = 1500.0 # стоимость позднего дефекта, найденного в QA
C_prod: float = 6000.0 # стоимость дефекта в production
cost_cv: float = 0.3 # индивидуальный разброс стоимости (лог-норм. сигма)
cost_est_noise: float = 0.0 # ошибка оценки C в решении (лог-норм. сигма)
c0: float = 500.0 # стоимость полного усиленного контроля
c_gamma: float = 0.0 # c_i = c0*(SP/5)^gamma; 0 -> константа как в ВКР
k_true: float = 1.6 # истинная эффективность: g(x)=exp(-k x)
k_mismatch: float = 1.0 # k_assumed = k_true * k_mismatch
sp_noise: float = 0.35 # шум оценки SP (лог-шкала)
use_component_feature: bool = False # расширенная модель: SP + компонент
threshold_T: float = 0.25 # порог базовой линии (2)
budget_mult: float = 1.0 # B = budget_mult * (расход пороговой политики)
budget_fixed: float | None = None # если задан, B не зависит от пороговой политики
# (нужно для свипов по T и по калибровке)
robust_z: float = 1.2816 # квантиль для верхней границы P (90%)
miscal_temp: float = 1.0 # искажение калибровки: logit -> t*logit
miscal_shift: float = 0.0 # и сдвиг: + b
shift_eta: float = 0.0 # сдвиг распределения: рост истинного риска в спринтах
shift_logsize: float = 0.0 # сдвиг распределения: задачи в спринтах крупнее
# ----------------------------------------------------------------------------
# 2. Генератор мира (независимый от модели)
# ----------------------------------------------------------------------------
def _nearest_fib(v: np.ndarray) -> np.ndarray:
idx = np.abs(np.log(v)[:, None] - np.log(FIB)[None, :]).argmin(axis=1)
return FIB[idx]
def _latent(cfg: Config, n: int, rng: np.random.Generator, size_shift: float = 0.0) -> dict:
size = rng.lognormal(mean=np.log(3.0) + size_shift, sigma=0.8, size=n)
sp = _nearest_fib(size * np.exp(rng.normal(0.0, cfg.sp_noise, n)))
comp = rng.choice(len(COMPONENTS), size=n, p=np.array(cfg.comp_share))
senior = rng.random(n) < 0.6
novelty = rng.random(n) < 0.2
eta_wo_a0 = (0.9 * np.log(size) + np.array(cfg.comp_shift)[comp]
+ 0.6 * novelty - 0.5 * senior)
return dict(size=size, sp=sp, comp=comp, eta=eta_wo_a0)
def calibrate_intercept(cfg: Config, seed: int = 12345) -> float:
"""Подбирает a0, чтобы средняя истинная P в мире равнялась base_rate."""
lat = _latent(cfg, 200_000, np.random.default_rng(seed))
lo, hi = -10.0, 5.0
for _ in range(80):
mid = 0.5 * (lo + hi)
p = 1.0 - np.exp(-np.exp(mid + lat["eta"]))
lo, hi = (mid, hi) if p.mean() < cfg.base_rate else (lo, mid)
return 0.5 * (lo + hi)
def make_tasks(cfg: Config, n: int, a0: float, rng: np.random.Generator,
size_shift: float = 0.0) -> pd.DataFrame:
lat = _latent(cfg, n, rng, size_shift)
p_true = 1.0 - np.exp(-np.exp(a0 + lat["eta"])) # cloglog-связь
comp = lat["comp"]
q = np.array(cfg.comp_prod_share)[comp]
kappa = np.array(cfg.comp_kappa)[comp]
c_expected = kappa * ((1 - q) * cfg.C_qa + q * cfg.C_prod) # известна по компоненту
s = cfg.cost_cv
c_true = c_expected * rng.lognormal(-s * s / 2, s, n) # индивидуальный разброс
c_prev = cfg.c0 * (lat["sp"] / 5.0) ** cfg.c_gamma
y = (rng.random(n) < p_true).astype(int)
return pd.DataFrame(dict(sp=lat["sp"], comp=comp, p_true=p_true, y=y,
C_exp=c_expected, C_true=c_true, c_prev=c_prev))
# ----------------------------------------------------------------------------
# 3. Модель вероятности и её качество
# ----------------------------------------------------------------------------
def design(df: pd.DataFrame, cfg: Config) -> np.ndarray:
cols = [df["sp"].to_numpy(float)]
if cfg.use_component_feature:
for j in range(1, len(COMPONENTS)):
cols.append((df["comp"].to_numpy() == j).astype(float))
return np.column_stack(cols)
class ProbModel:
"""Логистическая регрессия без SMOTE; дельта-метод для верхней границы P."""
def __init__(self, cfg: Config):
self.cfg = cfg
self.lr = LogisticRegression(C=1e6, max_iter=2000)
def fit(self, df: pd.DataFrame) -> "ProbModel":
X = design(df, self.cfg)
self.lr.fit(X, df["y"].to_numpy())
Xa = np.column_stack([np.ones(len(X)), X])
p = self.lr.predict_proba(X)[:, 1]
info = Xa.T @ (Xa * (p * (1 - p))[:, None])
self.cov = np.linalg.pinv(info)
return self
def _logit(self, df: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
X = design(df, self.cfg)
Xa = np.column_stack([np.ones(len(X)), X])
z = Xa @ np.r_[self.lr.intercept_, self.lr.coef_.ravel()]
se = np.sqrt(np.einsum("ij,jk,ik->i", Xa, self.cov, Xa))
return z, se
def predict(self, df: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
z, se = self._logit(df)
z = self.cfg.miscal_temp * z + self.cfg.miscal_shift # искажение калибровки
p_hat = 1 / (1 + np.exp(-z))
p_up = 1 / (1 + np.exp(-(z + self.cfg.robust_z * se)))
return p_hat, p_up
def ece(y: np.ndarray, p: np.ndarray, bins: int = 10) -> float:
edges = np.linspace(0, 1, bins + 1)
idx = np.clip(np.digitize(p, edges) - 1, 0, bins - 1)
total = 0.0
for b in range(bins):
m = idx == b
if m.any():
total += m.mean() * abs(y[m].mean() - p[m].mean())
return float(total)
def model_quality(model: ProbModel, test: pd.DataFrame) -> dict:
p_hat, _ = model.predict(test)
y = test["y"].to_numpy()
return dict(roc_auc=roc_auc_score(y, p_hat), pr_auc=average_precision_score(y, p_hat),
brier=brier_score_loss(y, p_hat), ece=ece(y, p_hat),
mean_p_hat=float(p_hat.mean()), base_rate=float(y.mean()))
# ----------------------------------------------------------------------------
# 4. Политики распределения контроля: каждая возвращает x in [0,1]^N
# ----------------------------------------------------------------------------
def g_exp(x: np.ndarray, k: float) -> np.ndarray:
return np.exp(-k * x)
def alloc_budget_exp(L: np.ndarray, c: np.ndarray, k: float, B: float | None,
tol: float = 1e-9) -> tuple[np.ndarray, float]:
"""min sum c x + L exp(-k x) s.t. sum c x <= B, 0<=x<=1. Возвращает (x, theta)."""
L = np.maximum(L, 0.0)
def x_of(theta: float) -> np.ndarray:
ratio = L * k / (theta * c)
with np.errstate(divide="ignore"):
x = np.where(ratio > 1.0, np.log(ratio) / k, 0.0)
return np.clip(x, 0.0, 1.0)
x1 = x_of(1.0)
if B is None or (c * x1).sum() <= B:
return x1, 1.0
lo, hi = 1.0, max(1.0, float((L * k / c).max())) * 1.0001
while hi - lo > tol * lo:
mid = 0.5 * (lo + hi)
if (c * x_of(mid)).sum() > B:
lo = mid
else:
hi = mid
return x_of(hi), hi
def alloc_linear_greedy(L: np.ndarray, c: np.ndarray, r: float, B: float) -> np.ndarray:
"""Абляция «без убывающей отдачи»: g(x)=1-(1-r)x -> дробный рюкзак."""
gain = L * (1 - r) - c # чистый выигрыш полного контроля
order = np.argsort(-(L * (1 - r)) / c)
x, spent = np.zeros_like(L), 0.0
for i in order:
if gain[i] <= 0 or spent >= B:
continue
take = min(1.0, (B - spent) / c[i])
x[i], spent = take, spent + take * c[i]
return x
def alloc_binary_knapsack_greedy(L, c, r, B):
"""Абляция «бинарный контроль»: x in {0,1}, жадно по w/c (точно при равных c)."""
w = L * (1 - r) - c
order = np.argsort(-w / c)
x, spent = np.zeros_like(L), 0.0
for i in order:
if w[i] > 0 and spent + c[i] <= B + 1e-9:
x[i], spent = 1.0, spent + c[i]
return x
def alloc_proportional_sp(sp: np.ndarray, c: np.ndarray, B: float) -> np.ndarray:
lo, hi = 0.0, 1.0 / sp.min()
for _ in range(100):
a = 0.5 * (lo + hi)
if (c * np.minimum(1, a * sp)).sum() > B:
hi = a
else:
lo = a
return np.minimum(1, lo * sp)
def alloc_random(c: np.ndarray, B: float, rng: np.random.Generator) -> np.ndarray:
x, spent = np.zeros_like(c), 0.0
for i in rng.permutation(len(c)):
if spent + c[i] <= B + 1e-9:
x[i], spent = 1.0, spent + c[i]
return x
def policies(df: pd.DataFrame, p_hat, p_up, cfg: Config, rng) -> dict:
c = df["c_prev"].to_numpy(float)
k_a = cfg.k_true * cfg.k_mismatch
r_a = np.exp(-k_a)
C_hat = df["C_exp"].to_numpy(float) * rng.lognormal(
-cfg.cost_est_noise ** 2 / 2, cfg.cost_est_noise, len(df)) if cfg.cost_est_noise > 0 \
else df["C_exp"].to_numpy(float)
x_thr = (p_hat > cfg.threshold_T).astype(float) # (2) текущий подход ВКР
B = cfg.budget_fixed if cfg.budget_fixed is not None \
else cfg.budget_mult * float((c * x_thr).sum())
L_hat, L_up = p_hat * C_hat, p_up * C_hat
out = {
"1a_no_control": np.zeros(len(df)),
"1b_uniform_same_B": np.full(len(df), min(1.0, B / c.sum())),
"2_threshold_T": x_thr,
"3_random_same_B": alloc_random(c, B, rng),
"4_proportional_SP": alloc_proportional_sp(df["sp"].to_numpy(float), c, B),
"5_proposed": alloc_budget_exp(L_hat, c, k_a, B)[0],
"5r_proposed_robust": alloc_budget_exp(L_up, c, k_a, B)[0],
"A1_no_cost_heterog": alloc_budget_exp(p_hat * C_hat.mean(), c, k_a, B)[0],
"A2_no_budget": alloc_budget_exp(L_hat, c, k_a, None)[0],
"A3_no_dimin_returns": alloc_linear_greedy(L_hat, c, r_a, B),
"A5_binary_levels": alloc_binary_knapsack_greedy(L_hat, c, r_a, B),
# оракул при том же бюджете B — верхняя граница для политик 1b–5, A1, A3, A5
"oracle": alloc_budget_exp(df["p_true"].to_numpy() * df["C_true"].to_numpy(),
c, cfg.k_true, B)[0],
# оракул без бюджета — сравнивать только с A2 (у A2 расход не ограничен B)
"oracle_no_budget": alloc_budget_exp(df["p_true"].to_numpy() * df["C_true"].to_numpy(),
c, cfg.k_true, None)[0],
}
return out, B
# ----------------------------------------------------------------------------
# 5. Оценка политик в ИСТИННОМ мире (общие случайные числа)
# ----------------------------------------------------------------------------
def evaluate(df: pd.DataFrame, x: np.ndarray, cfg: Config, u: np.ndarray) -> dict:
c = df["c_prev"].to_numpy(float)
p, C = df["p_true"].to_numpy(), df["C_true"].to_numpy()
p_res = p * g_exp(x, cfg.k_true)
L_true = p * C
exp_cost = float((c * x).sum() + (p_res * C).sum())
realized_defects = u < p_res # CRN: один u на задачу для всех политик
return dict(
spend=float((c * x).sum()),
exp_late_defects=float(p_res.sum()),
exp_total_cost=exp_cost,
avoided_loss=float((L_true * (1 - g_exp(x, cfg.k_true))).sum()),
recall_at_budget_lossw=float(L_true[x > 1e-6].sum() / L_true.sum()),
realized_cost=float((c * x).sum() + C[realized_defects].sum()),
)
def run(cfg: Config, reps: int, sprints: int, seed: int = 42) -> tuple[pd.DataFrame, pd.DataFrame]:
a0 = calibrate_intercept(cfg)
rows, qrows = [], []
for rep in range(reps):
rng = np.random.default_rng([seed, rep])
hist = make_tasks(cfg, cfg.n_hist, a0, rng)
idx = rng.permutation(len(hist))
n_te = int(cfg.test_share * len(hist))
test, train = hist.iloc[idx[:n_te]], hist.iloc[idx[n_te:]]
model = ProbModel(cfg).fit(train)
qrows.append(dict(rep=rep, **model_quality(model, test)))
for s in range(sprints):
# спринты могут генерироваться со сдвигом относительно истории (устойчивость к сдвигу)
sprint = make_tasks(cfg, cfg.n_sprint, a0 + cfg.shift_eta, rng, cfg.shift_logsize)
p_hat, p_up = model.predict(sprint)
pol, B = policies(sprint, p_hat, p_up, cfg, rng)
u = rng.random(len(sprint))
for name, x in pol.items():
rows.append(dict(rep=rep, sprint=s, policy=name, budget=B,
**evaluate(sprint, x, cfg, u)))
return pd.DataFrame(rows), pd.DataFrame(qrows)
def summarize(res: pd.DataFrame, ref: str = "2_threshold_T") -> pd.DataFrame:
"""Средние по повторам и парные разности с базовой линией (2): 95% bootstrap ДИ + Wilcoxon."""
per_rep = res.groupby(["rep", "policy"]).mean(numeric_only=True).reset_index()
base = per_rep[per_rep.policy == ref].set_index("rep")["exp_total_cost"]
oracle = per_rep[per_rep.policy == "oracle"].set_index("rep")["exp_total_cost"]
out = []
rng = np.random.default_rng(0)
for name, g in per_rep.groupby("policy"):
g = g.set_index("rep")
d = (g["exp_total_cost"] - base).to_numpy()
boot = [rng.choice(d, len(d)).mean() for _ in range(2000)]
p_w = stats.wilcoxon(d).pvalue if name != ref and np.any(d != 0) else np.nan
out.append(dict(
policy=name,
exp_total_cost=g["exp_total_cost"].mean(),
spend=g["spend"].mean(),
exp_late_defects=g["exp_late_defects"].mean(),
avoided_loss=g["avoided_loss"].mean(),
recall_at_budget_lossw=g["recall_at_budget_lossw"].mean(),
regret_vs_oracle=(g["exp_total_cost"] - oracle).mean(),
diff_vs_threshold=d.mean(),
ci95_low=np.percentile(boot, 2.5), ci95_high=np.percentile(boot, 97.5),
wilcoxon_p=p_w,
realized_cost_p95=res[res.policy == name]["realized_cost"].quantile(0.95),
))
return pd.DataFrame(out).sort_values("exp_total_cost")
# ----------------------------------------------------------------------------
# 6. Свипы чувствительности
# ----------------------------------------------------------------------------
SWEEPS = {
"budget": ("budget_mult", [0.25, 0.5, 1.0, 1.5, 2.0]),
"threshold": ("threshold_T", [0.10, 0.15, 0.20, 0.25, 0.30, 0.40]),
"k_mismatch": ("k_mismatch", [0.5, 0.75, 1.0, 1.5, 2.0]),
"cost_cv": ("cost_cv", [0.0, 0.3, 0.6, 1.0]),
"cost_est_noise": ("cost_est_noise", [0.0, 0.25, 0.5, 1.0]),
"miscalibration": ("miscal_temp", [0.5, 0.75, 1.0, 1.5, 2.0]),
"underestimation": ("miscal_shift", [-1.0, -0.5, 0.0, 0.5]),
"k_true": ("k_true", [0.5, 1.0, 1.6, 2.5]),
"shift": ("shift_eta", [0.0, 0.25, 0.5]),
"shift_size": ("shift_logsize", [0.0, 0.3, 0.6]),
}
def run_sweep(base: Config, key: str, reps: int, sprints: int, out: str) -> pd.DataFrame:
field, values = SWEEPS[key]
frames = []
for v in values:
cfg = replace(base, **{field: v})
res, _ = run(cfg, reps, sprints)
s = summarize(res)
s.insert(0, field, v)
frames.append(s)
df = pd.concat(frames, ignore_index=True)
df.to_csv(os.path.join(out, f"sweep_{key}.csv"), index=False)
return df
# ----------------------------------------------------------------------------
# 7. Проверочный расчёт примера §2.4 текущей редакции ВКР (аудит)
# ----------------------------------------------------------------------------
def legacy_check() -> None:
R = lambda c: 0.5 / (1 + np.exp(2.8 - 0.38 * c))
counts = {1: 25, 2: 28, 3: 20, 5: 15, 8: 9, 13: 3} # рис. 2.2, N=100
sp = np.repeat(list(counts), list(counts.values())).astype(float)
r = R(sp)
hi = r > 0.25
no_ctrl = r.sum()
ml = r[~hi].sum() + 0.10 * hi.sum()
print(f"Средний R по рис. 2.2: {r.mean():.4f} (в тексте база 0,15)")
print(f"Ожид. поздних дефектов без контроля в ТОМ ЖЕ мире: {no_ctrl:.2f} -> {3000*no_ctrl:,.0f} у.е.")
print(f"ML-правило T=0,25: {ml:.2f} дефектов, затраты {3000*ml + 500*hi.sum():,.0f} у.е.")
print(f"Экономия к согласованной базе: {3000*no_ctrl - 3000*ml - 500*hi.sum():,.0f} у.е.")
def control_check(seed: int = 7) -> None:
"""Контрольное предсказание H1_0 (Утверждение 2б): при однородной цене дефекта,
равных c_i, бинарном контроле и неактивном бюджете оптимальное правило
совпадает с порогом T* = c / ((1 - r) C). Расхождение = ошибка реализации."""
cfg = Config(cost_cv=0.0, comp_kappa=(1.0,) * 5, comp_prod_share=(0.3,) * 5)
a0 = calibrate_intercept(cfg)
rng = np.random.default_rng(seed)
model = ProbModel(cfg).fit(make_tasks(cfg, cfg.n_hist, a0, rng))
sprint = make_tasks(cfg, 500, a0, rng)
p_hat, _ = model.predict(sprint)
C = float(sprint["C_exp"].iloc[0]); c = sprint["c_prev"].to_numpy(float)
r = np.exp(-cfg.k_true)
x_thr = (p_hat > cfg.c0 / ((1 - r) * C)).astype(float)
x_opt = alloc_binary_knapsack_greedy(p_hat * C, c, r, B=np.inf)
print("T* =", round(cfg.c0 / ((1 - r) * C), 4), "| совпадение выбора:", bool(np.all(x_thr == x_opt)))
# ----------------------------------------------------------------------------
def main() -> None:
ap = argparse.ArgumentParser()
ap.add_argument("--reps", type=int, default=50)
ap.add_argument("--sprints", type=int, default=30)
ap.add_argument("--out", default="results")
ap.add_argument("--sweep", choices=list(SWEEPS), default=None)
ap.add_argument("--component-feature", action="store_true")
ap.add_argument("--legacy", action="store_true")
ap.add_argument("--budget-fixed", type=float, default=None,
help="фиксированный бюджет B на спринт, у.е. (для свипов threshold/miscalibration)")
ap.add_argument("--control-check", action="store_true",
help="контрольное предсказание: частный случай должен совпасть с порогом")
a = ap.parse_args()
if a.legacy:
legacy_check()
return
if a.control_check:
control_check()
return
os.makedirs(a.out, exist_ok=True)
cfg = Config(use_component_feature=a.component_feature, budget_fixed=a.budget_fixed)
pd.Series(dataclasses.asdict(cfg)).to_csv(os.path.join(a.out, "config.csv"))
if a.sweep:
print(run_sweep(cfg, a.sweep, a.reps, a.sprints, a.out).to_string(index=False))
return
res, q = run(cfg, a.reps, a.sprints)
res.to_csv(os.path.join(a.out, "raw.csv"), index=False)
q.to_csv(os.path.join(a.out, "model_quality.csv"), index=False)
s = summarize(res)
s.to_csv(os.path.join(a.out, "summary.csv"), index=False)
pd.set_option("display.width", 200)
print(q.describe().loc[["mean", "std"]].to_string())
print(s.to_string(index=False))
if __name__ == "__main__":
main()