import math
import numpy as np
import pandas as pd
y = [1.55, 3.80, 5.40, 2.60, 5.60, 8.30, 7.30, 8.40, 4.20, 5.60]
x1 = [0.30, 0.60, 1.10, 0.20, 0.80, 1.60, 1.50, 0.90, 0.90, 0.80]
x2 = [0.80, 1.90, 2.70, 1.30, 2.80, 4.10, 3.60, 4.20, 2.10, 2.80]
def _col(v):
v = np.asarray(v).reshape(-1)
return v[:, None]
def design_matrix(x1, x2, intercept=True):
x1c = _col(x1)
x2c = _col(x2)
if intercept:
ones = np.ones_like(x1c)
X = np.hstack([ones, x1c, x2c])
else:
X = np.hstack([x1c, x2c])
return X
def safe_inv(M):
try:
return np.linalg.inv(M)
except np.linalg.LinAlgError:
return np.linalg.pinv(M)
def _log_beta(a, b):
return math.lgamma(a) + math.lgamma(b) - math.lgamma(a + b)
def _betacf(a, b, x, max_iter=200, eps=3e-14):
am, bm = 1.0, 1.0
az = 1.0
qab = a + b
qap = a + 1.0
qam = a - 1.0
bz = 1.0 - qab * x / qap
if abs(bz) < 1e-300:
bz = 1e-300
em, tem = 0.0, 0.0
d = 0.0
ap, bp = 0.0, 0.0
for m in range(1, max_iter + 1):
em = float(m)
tem = em + em
d = em * (b - em) * x / ((qam + tem) * (a + tem))
ap = az + d * am
bp = bz + d * bm
if abs(bp) < 1e-300:
bp = 1e-300
d = -(a + em) * (qab + em) * x / ((a + tem) * (qap + tem))
app = ap + d * az
bpp = bp + d * bz
if abs(bpp) < 1e-300:
bpp = 1e-300
am, bm = ap / bpp, bp / bpp
az, bz = app / bpp, 1.0
if abs(app - ap) < eps * abs(app):
return az
return az
def betainc_regularized(a, b, x):
if x <= 0.0:
return 0.0
if x >= 1.0:
return 1.0
ln_bt = a * math.log(x) + b * math.log(1.0 - x) - _log_beta(a, b)
if x < (a + 1.0) / (a + b + 2.0):
cf = _betacf(a, b, x)
return math.exp(ln_bt) * cf / a
else:
cf = _betacf(b, a, 1.0 - x)
return 1.0 - math.exp((b * math.log(x) + a * math.log(1.0 - x)) - _log_beta(b, a)) * cf / b
def t_cdf(t, v):
if v <= 0:
return float("nan")
if t == 0.0:
return 0.5
x = v / (v + t * t)
ib = betainc_regularized(v / 2.0, 0.5, x)
if t > 0:
return 1.0 - 0.5 * ib
else:
return 0.5 * ib
def t_ppf(p, v, tol=1e-10, max_iter=200):
if not (0.0 < p < 1.0):
return float("nan")
if p == 0.5:
return 0.0
lo, hi = -50.0, 50.0
for _ in range(max_iter):
mid = (lo + hi) / 2.0
c = t_cdf(mid, v)
if abs(c - p) < tol:
return mid
if c < p:
lo = mid
else:
hi = mid
return (lo + hi) / 2.0
def f_cdf(x, df1, df2):
if x <= 0.0:
return 0.0
z = (df1 * x) / (df1 * x + df2)
return betainc_regularized(df1 / 2.0, df2 / 2.0, z)
def f_ppf(p, df1, df2, tol=1e-10, max_iter=200):
if not (0.0 < p < 1.0):
return float("nan")
lo, hi = 0.0, 1.0
while f_cdf(hi, df1, df2) < p and hi < 1e8:
hi *= 2.0
for _ in range(max_iter):
mid = (lo + hi) / 2.0
c = f_cdf(mid, df1, df2)
if abs(c - p) < tol:
return mid
if c < p:
lo = mid
else:
hi = mid
return (lo + hi) / 2.0
def ols_fit(X, y):
y = _col(y)
n, p = X.shape
XtX = X.T @ X
XtX_inv = safe_inv(XtX)
Xty = X.T @ y
beta = XtX_inv @ Xty
yhat = X @ beta
e = y - yhat
RSS = float((e.T @ e).ravel()[0])
ybar = float(np.mean(y))
TSS = float(((y - ybar) ** 2).sum())
df = n - p
if df <= 0:
raise ValueError("Недостатньо ступенів вільності (n - p має бути > 0).")
s2 = RSS / df
cov_beta = s2 * XtX_inv
se = np.sqrt(np.diag(cov_beta)).reshape(-1, 1)
with np.errstate(divide="ignore", invalid="ignore"):
tvals = beta / se
tvals = np.where(np.isfinite(tvals), tvals, np.nan)
pvals = np.empty_like(tvals)
for i in range(tvals.size):
ti = float(abs(tvals[i, 0]))
pvals[i, 0] = 2.0 * (1.0 - t_cdf(ti, df))
tcrit = t_ppf(0.975, df)
ci_low = beta - tcrit * se
ci_high = beta + tcrit * se
R2 = 1.0 - RSS / TSS if TSS > 0 else float("nan")
k = p - 1
denom = (1.0 - R2) / (n - p) if (n - p) > 0 else float("nan")
numer = R2 / k if k > 0 else float("nan")
F = numer / denom if denom not in (None, 0.0) else float("nan")
F_pvalue = 1.0 - f_cdf(F, k, n - p) if math.isfinite(F) else float("nan")
F_crit = f_ppf(0.95, k, n - p) if k > 0 else float("nan")
return dict(
beta=beta, yhat=yhat, e=e, RSS=RSS, TSS=TSS, s2=s2, se=se,
tvals=tvals, pvals=pvals, ci_low=ci_low, ci_high=ci_high,
XtX=XtX, XtX_inv=XtX_inv, cov_beta=cov_beta,
R2=R2, F=F, F_pvalue=F_pvalue, F_crit=F_crit,
df_resid=df, k=k, p=p, n=n, tcrit=tcrit
)
def standardized_betas(y, x1, x2, beta):
y = np.asarray(y).reshape(-1)
x1 = np.asarray(x1).reshape(-1)
x2 = np.asarray(x2).reshape(-1)
sy = np.std(y, ddof=1)
sx1 = np.std(x1, ddof=1)
sx2 = np.std(x2, ddof=1)
b1, b2 = float(beta[1, 0]), float(beta[2, 0])
beta1_std = b1 * (sx1 / sy) if sy > 0 else float("nan")
beta2_std = b2 * (sx2 / sy) if sy > 0 else float("nan")
return beta1_std, beta2_std
def elasticities(y, x1, x2, beta):
y = np.asarray(y).reshape(-1)
x1 = np.asarray(x1).reshape(-1)
x2 = np.asarray(x2).reshape(-1)
my = float(np.mean(y))
mx1 = float(np.mean(x1))
mx2 = float(np.mean(x2))
b1, b2 = float(beta[1, 0]), float(beta[2, 0])
E1 = b1 * (mx1 / my) if my != 0 else float("nan")
E2 = b2 * (mx2 / my) if my != 0 else float("nan")
return E1, E2
def prediction_intervals(x_new, fit):
x_new = np.asarray(x_new).reshape(-1, 1)
beta = fit["beta"]
XtX_inv = fit["XtX_inv"]
s2 = float(fit["s2"])
s = math.sqrt(s2)
tcrit = fit["tcrit"]
df = int(fit["df_resid"])
yhat_new = float((x_new.T @ beta).ravel()[0])
leverage = float((x_new.T @ XtX_inv @ x_new).ravel()[0])
stderr_mean = float(s * math.sqrt(leverage))
CI_low = yhat_new - tcrit * stderr_mean
CI_high = yhat_new + tcrit * stderr_mean
stderr_pred = float(s * math.sqrt(1.0 + leverage))
PI_low = yhat_new - tcrit * stderr_pred
PI_high = yhat_new + tcrit * stderr_pred
return dict(
yhat=yhat_new,
CI95_low=CI_low, CI95_high=CI_high,
PI95_low=PI_low, PI95_high=PI_high,
df=df, tcrit=float(tcrit)
)
def fit_and_report(y, x1, x2):
src_df = pd.DataFrame({"y": y, "x1": x1, "x2": x2})
print("=== Вихідні дані (Варіант 2) ===")
print(src_df.round(4).to_string(index=True))
print(f"n = {len(y)} спостережень")
X = design_matrix(x1, x2, intercept=True)
fit = ols_fit(X, y)
print("\n--- Матриця X'X ---")
print(pd.DataFrame(fit["XtX"], index=["1","x1","x2"], columns=["1","x1","x2"]).round(6).to_string())
print("\n--- (X'X)^(-1) ---")
print(pd.DataFrame(fit["XtX_inv"], index=["a0","a1","a2"], columns=["a0","a1","a2"]).round(6).to_string())
a0, a1, a2 = (float(fit["beta"][0,0]), float(fit["beta"][1,0]), float(fit["beta"][2,0]))
print("\n=== (1) Вибіркове рівняння множинної регресії ===")
print(f"ŷ = {a0:.6f} + {a1:.6f}·x1 + {a2:.6f}·x2")
print("\n=== (2) Дисперсія збурень та коваріаційна матриця ===")
print(f"s_e^2 = {fit['s2']:.6f} (df = {fit['df_resid']})")
print("Cov(a) = s_e^2 · (X'X)^(-1):")
print(pd.DataFrame(fit["cov_beta"], index=["a0","a1","a2"], columns=["a0","a1","a2"]).round(6).to_string())
print("\n=== (3) t-тести для коефіцієнтів (α=0.05) та 95% ДІ ===")
coef_tbl = pd.DataFrame({
"coef":[a0, a1, a2],
"SE": fit["se"].reshape(-1),
"t": fit["tvals"].reshape(-1),
"p": fit["pvals"].reshape(-1),
"CI95_low": fit["ci_low"].reshape(-1),
"CI95_high": fit["ci_high"].reshape(-1),
}, index=["a0","a1","a2"])
print(coef_tbl.round(6).to_string())
R2 = float(fit["R2"])
R = float(math.sqrt(R2)) if R2 >= 0 else float("nan")
print("\n=== (4) Дисперсія та значущість моделі ===")
print(f"R^2 = {R2:.6f}, R = {R:.6f}")
print(f"F = {fit['F']:.6f} (df1 = {fit['k']}, df2 = {fit['df_resid']}), p-value = {fit['F_pvalue']:.6g}")
print(f"F_crit (α=0.05) = {fit['F_crit']:.6f}")
b1_std, b2_std = standardized_betas(y, x1, x2, fit["beta"])
E1, E2 = elasticities(y, x1, x2, fit["beta"])
comp_df = pd.DataFrame({"std_beta":[b1_std, b2_std], "elasticity":[E1, E2]},
index=["x1 (тис. м²)", "x2 (тис. люд/день)"])
print("\n=== (5) Порівняння роздільного впливу факторів ===")
print(comp_df.round(6).to_string())
stronger_std = "x1" if abs(b1_std) > abs(b2_std) else "x2"
stronger_el = "x1" if abs(E1) > abs(E2) else "x2"
print(f"→ За |std_beta| сильніший вплив: {stronger_std}; за |elasticity|: {stronger_el}.")
x_new = np.array([1.0, 1.2, 15.0])
pred = prediction_intervals(x_new, fit)
print("\n=== (6) Прогноз у точці x1=1.2, x2=15 з 95% інтервалами ===")
print(f"ŷ_new = {pred['yhat']:.6f}")
print(f"95% CI (середній відгук): [{pred['CI95_low']:.6f}; {pred['CI95_high']:.6f}]")
print(f"95% PI (нове спостереження): [{pred['PI95_low']:.6f}; {pred['PI95_high']:.6f}]")
print(f"(t_crit = {pred['tcrit']:.4f}, df = {pred['df']})")
print("\n=== ПІДСУМОК ===")
print(f"Модель: ŷ = {a0:.4f} + {a1:.4f}·x1 + {a2:.4f}·x2")
print(f"R^2 = {R2:.4f}, F({fit['k']},{fit['df_resid']}) = {fit['F']:.4f} (F_crit={fit['F_crit']:.4f})")
print(f"Прогноз у (1.2; 15): ŷ = {pred['yhat']:.4f}, 95% PI = [{pred['PI95_low']:.4f}; {pred['PI95_high']:.4f}]")
print("x1 — тис. м², x2 — тис. людей/день; y — млн грн.")
if __name__ == "__main__":
fit_and_report(y, x1, x2)
