13.4-dars: Taxminlar va diagnostika

1. CHIZIQLILIK     — y va belgilar orasidagi bog'liqlik chiziqli
      buzilsa: BASHORAT ham, talqin ham noto'g'ri
2. MUSTAQILLIK     — kuzatuvlar (qoldiqlar) o'zaro bog'liq emas
      buzilsa: standart xatolar noto'g'ri (vaqt qatori, klasterlar)
3. GOMOSKEDASTIKLIK — qoldiq dispersiyasi barcha bashoratlarda bir xil
      buzilsa: koeffitsiyent baholari to'g'ri, CI va p NOTO'G'RI
4. NORMAL QOLDIQLAR — qoldiqlar taxminan normal taqsimlangan
      buzilsa: kichik namunada CI noto'g'ri; katta namunada muhim emas 4.7-bob
5. MULTIKOLLINEARLIK YO'Q — belgilar o'zaro kuchli bog'liq emas 13.3-bob
      buzilsa: koeffitsiyentlar beqaror, bashorat zarar ko'rmaydi

Asosiy grafiklar (10.6, 12.6):
  1. qoldiq ~ bashorat    — egrilik (chiziqlilik), konus (geteroskedastiklik)
  2. Q-Q grafigi          — normallik
  3. qoldiq ~ har belgi   — qaysi belgi shaklni buzyapti
  4. qoldiq ~ tartib/vaqt — mustaqillik (avtokorrelyatsiya)
  5. leverage ~ qoldiq    — ta'sirli nuqtalar

Sonli testlar: Breusch-Pagan (geteroskedastiklik), Durbin-Watson (avtokorrelyatsiya),
               Jarque-Bera (normallik), RESET (shakl)

Belgisi: qoldiq tarqoqligi bashorat bilan birga o'sadi (konus shakli)
Sabablari: ko'paytiruvchi xato, y ning o'ng qiyshiqligi, birlashtirilgan ma'lumot

Oqibati: koeffitsiyentlar BIAS SIZ, lekin standart xatolar noto'g'ri
         (odatda juda kichik → CI tor, p-qiymat kichik → soxta ishonch)

Yechimlar:
  1. log(y) yoki sqrt(y) transformatsiyasi
  2. ROBUST (HC3) standart xatolar — cov_type="HC3"
  3. Vaznli eng kichik kvadratlar (WLS)
  4. GLM (Gamma, Poisson — 13.10)

NORMALLIK KERAK:
  · kichik namunada (n < 30-50) aniq CI va p-qiymat uchun
  · bashorat intervallari (prediction interval) uchun

NORMALLIK DEYARLI KERAK EMAS:
  · katta namunada — MLT tufayli koeffitsiyent taqsimoti normal 4.7-bob
  · nuqtaviy bashorat uchun

Tekshirish: Q-Q grafigi (eng yaxshi), gistogramma; testlar — katta n da haddan tashqari sezgir

leverage (h_ii) — kuzatuv X fazosida qanchalik chetda
  o'rtacha leverage = p/n;  2p/n dan katta — yuqori

standartlashtirilgan qoldiq — y bo'yicha qanchalik chetda

Cook masofasi — ikkalasini birlashtiradi: "shu nuqta olib tashlansa,
                barcha bashoratlar qancha o'zgaradi"
  D > 4/n  — e'tibor bering;  D > 1 — juda ta'sirli

⚠️ Yuqori leverage o'zi yomon emas; ta'sirli nuqta = yuqori leverage + katta qoldiq

BUZILISH                TA'SIR              YECHIM
chiziqlilik             hamma narsa         nochiziqli belgilar 13.5-bob, boshqa model
mustaqillik             CI, p               klaster/HAC standart xatolar, aralash model
geteroskedastiklik      CI, p               HC3, log(y), WLS, GLM
normallik (kichik n)    CI, p               bootstrap 11.8-bob, transformatsiya
multikollinearlik       koeffitsiyent        birlashtirish, Ridge 13.7-bob
ta'sirli nuqtalar       hammasi             tekshirish, robust regressiya (13.11)

import numpy as np
import statsmodels.api as sm
from scipy import stats
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import durbin_watson

nat = sm.OLS(y, sm.add_constant(X)).fit()
qoldiq, bashorat = nat.resid, nat.fittedvalues

het_breuschpagan(qoldiq, nat.model.exog)[1]        # geteroskedastiklik p
durbin_watson(qoldiq)                              # ~2 — avtokorrelyatsiya yo'q
stats.probplot(qoldiq, dist="norm")                # Q-Q ma'lumoti

ta = nat.get_influence()
ta.hat_matrix_diag                                 # leverage
ta.cooks_distance[0]                               # Cook masofasi

robust = sm.OLS(y, sm.add_constant(X)).fit(cov_type="HC3")   # robust SE
QOIDA: avval chiziqlilik · qoldiq grafigi · HC3 arzon · outlierni o'chirma

Chiziqlilik — eng muhim · Gomoskedastiklik va mustaqillik — CI/p
Normallik — katta n da muhim emas · Multikollinearlik — koeffitsiyent
Qoldiq~bashorat grafigi — markaziy vosita · HC3 — eng arzon yechim
Cook > 4/n — tekshiring (o'chirmang)

"""Taxminlarning turli buzilishlari (real numpy/statsmodels)."""

import numpy as np
import statsmodels.api as sm
from scipy import stats
from statsmodels.stats.diagnostic import het_breuschpagan


def main() -> None:
    rng = np.random.default_rng(12)
    n = 600
    x = rng.uniform(1, 20, n)

    holatlar = {
        "yaxshi": 5 + 2.0 * x + rng.normal(0, 4, n),
        "egrilik": 5 + 0.22 * x ** 2 + rng.normal(0, 4, n),
        "konus": 5 + 2.0 * x + rng.normal(0, 0.55 * x, n),
        "og'ir dum": 5 + 2.0 * x + rng.standard_t(2, n) * 3,
    }

    print("=== 1. Umumiy ko'rsatkichlar ===")
    print(f"  {'holat':<10} {'R^2':>7} {'qoldiq SD':>10}")
    natijalar = {}
    for nom, y in holatlar.items():
        nat = sm.OLS(y, sm.add_constant(x)).fit()
        natijalar[nom] = nat
        print(f"  {nom:<10} {nat.rsquared:>7.3f} {np.sqrt(nat.mse_resid):>10.2f}")

    print("\n=== 2. Qoldiq ~ bashorat: choraklar bo'yicha o'rtacha ===")
    for nom, nat in natijalar.items():
        b, e = nat.fittedvalues, nat.resid
        chegara = np.quantile(b, [0, 0.25, 0.5, 0.75, 1.0])
        ort = [e[(b >= chegara[i]) & (b <= chegara[i + 1])].mean() for i in range(4)]
        print(f"  {nom:<10}: [{', '.join(f'{v:6.2f}' for v in ort)}]")
    print("  egrilik holatida naqsh: + - - + ")

    print("\n=== 3. Geteroskedastiklik (Breusch-Pagan) ===")
    for nom, nat in natijalar.items():
        p = het_breuschpagan(nat.resid, nat.model.exog)[1]
        b, e = nat.fittedvalues, nat.resid
        nisbat = e[b >= np.median(b)].std() / e[b < np.median(b)].std()
        print(f"  {nom:<10}: p = {p:.2e}, SD nisbati = {nisbat:.2f}")

    print("\n=== 4. Normallik (Jarque-Bera) ===")
    for nom, nat in natijalar.items():
        jb = stats.jarque_bera(nat.resid)
        print(f"  {nom:<10}: p = {jb.pvalue:.2e}, "
              f"basiqlik = {stats.kurtosis(nat.resid):6.2f}")
    print("  ⭐ R^2 o'xshash, lekin muammolar butunlay boshqa")


if __name__ == "__main__":
    main()

=== 1. Umumiy ko'rsatkichlar ===
  holat          R^2  qoldiq SD
  yaxshi       0.884       4.09
  egrilik      0.937       6.96
  konus        0.764       6.51
  og'ir dum    0.343      15.74

=== 2. Qoldiq ~ bashorat: choraklar bo'yicha o'rtacha ===
  yaxshi    : [  0.03,   0.04,  -0.05,  -0.02]
  egrilik   : [  4.99,  -5.22,  -4.62,   4.84]
  konus     : [ -0.12,   0.09,   0.29,  -0.26]
  og'ir dum : [  0.83,  -1.17,   0.11,   0.23]
  egrilik holatida naqsh: + - - +

=== 3. Geteroskedastiklik (Breusch-Pagan) ===
  yaxshi    : p = 6.35e-01, SD nisbati = 1.02
  egrilik   : p = 1.86e-02, SD nisbati = 0.92
  konus     : p = 3.26e-30, SD nisbati = 2.71
  og'ir dum : p = 5.61e-01, SD nisbati = 0.27

=== 4. Normallik (Jarque-Bera) ===
  yaxshi    : p = 2.03e-01, basiqlik =  -0.01
  egrilik   : p = 9.18e-03, basiqlik =  -0.50
  konus     : p = 1.78e-11, basiqlik =   1.41
  og'ir dum : p = 0.00e+00, basiqlik = 380.66
  ⭐ R^2 o'xshash, lekin muammolar butunlay boshqa

"""Oddiy va robust standart xatolar (real numpy/statsmodels)."""

import numpy as np
import statsmodels.api as sm


def main() -> None:
    rng = np.random.default_rng(21)
    HAQIQIY = 2.0

    print("=== 1. Bitta namunada ikki xil standart xato ===")
    n = 500
    x = rng.uniform(1, 20, n)
    y = 5 + HAQIQIY * x + rng.normal(0, 0.55 * x, n)
    oddiy = sm.OLS(y, sm.add_constant(x)).fit()
    robust = sm.OLS(y, sm.add_constant(x)).fit(cov_type="HC3")
    print(f"  koeffitsiyent (ikkalasida bir xil): {oddiy.params[1]:.4f} / "
          f"{robust.params[1]:.4f}")
    print(f"  oddiy SE  = {oddiy.bse[1]:.4f}, CI = "
          f"[{oddiy.conf_int()[1, 0]:.3f}, {oddiy.conf_int()[1, 1]:.3f}]")
    print(f"  robust SE = {robust.bse[1]:.4f}, CI = "
          f"[{robust.conf_int()[1, 0]:.3f}, {robust.conf_int()[1, 1]:.3f}]")

    print("\n=== 2. Qamrov tajribasi (500 takrorlash) ===")
    qamrov_oddiy = qamrov_robust = 0
    N = 500
    for _ in range(N):
        xs = rng.uniform(1, 20, n)
        ys = 5 + HAQIQIY * xs + rng.normal(0, 0.55 * xs, n)
        Xc = sm.add_constant(xs)
        o = sm.OLS(ys, Xc).fit()
        r = sm.OLS(ys, Xc).fit(cov_type="HC3")
        lo, hi = o.conf_int()[1]
        qamrov_oddiy += lo <= HAQIQIY <= hi
        lo, hi = r.conf_int()[1]
        qamrov_robust += lo <= HAQIQIY <= hi
    print(f"  oddiy CI qamrovi : {qamrov_oddiy / N:.1%}  (kerak 95%)")
    print(f"  robust CI qamrovi: {qamrov_robust / N:.1%}")

    print("\n=== 3. Gomoskedastik holatda farq yo'q ===")
    y2 = 5 + HAQIQIY * x + rng.normal(0, 6, n)
    o2 = sm.OLS(y2, sm.add_constant(x)).fit()
    r2 = sm.OLS(y2, sm.add_constant(x)).fit(cov_type="HC3")
    print(f"  oddiy SE = {o2.bse[1]:.4f}, robust SE = {r2.bse[1]:.4f}")
    print("  (robust SE deyarli bepul — shubha bo'lsa ishlating)")

    print("\n=== 4. log(y) transformatsiyasi ===")
    y3 = np.exp(1.5 + 0.12 * x) * rng.lognormal(0, 0.3, n)
    xom = sm.OLS(y3, sm.add_constant(x)).fit()
    log_m = sm.OLS(np.log(y3), sm.add_constant(x)).fit()
    e1, b1 = xom.resid, xom.fittedvalues
    e2, b2 = log_m.resid, log_m.fittedvalues
    print(f"  xom y:    SD nisbati = "
          f"{e1[b1 >= np.median(b1)].std() / e1[b1 < np.median(b1)].std():.2f}")
    print(f"  log(y):   SD nisbati = "
          f"{e2[b2 >= np.median(b2)].std() / e2[b2 < np.median(b2)].std():.2f}")
    print("  ⭐ Ko'paytiruvchi xatoda log(y) muammoni yo'qotadi")


if __name__ == "__main__":
    main()

=== 1. Bitta namunada ikki xil standart xato ===
  koeffitsiyent (ikkalasida bir xil): 1.9931 / 1.9931
  oddiy SE  = 0.0541, CI = [1.887, 2.099]
  robust SE = 0.0601, CI = [1.875, 2.111]

=== 2. Qamrov tajribasi (500 takrorlash) ===
  oddiy CI qamrovi : 92.6%  (kerak 95%)
  robust CI qamrovi: 94.4%

=== 3. Gomoskedastik holatda farq yo'q ===
  oddiy SE = 0.0465, robust SE = 0.0456
  (robust SE deyarli bepul — shubha bo'lsa ishlating)

=== 4. log(y) transformatsiyasi ===
  xom y:    SD nisbati = 2.50
  log(y):   SD nisbati = 1.06
  ⭐ Ko'paytiruvchi xatoda log(y) muammoni yo'qotadi

"""Leverage, qoldiq va Cook masofasi (real numpy/statsmodels)."""

import numpy as np
import statsmodels.api as sm


def main() -> None:
    rng = np.random.default_rng(31)
    n = 100
    x = rng.uniform(0, 10, n)
    y = 4 + 1.5 * x + rng.normal(0, 1.2, n)

    print("=== 1. Asl model ===")
    nat = sm.OLS(y, sm.add_constant(x)).fit()
    print(f"  nishab = {nat.params[1]:.3f}")

    holatlar = {
        "markazda katta qoldiq": (5.0, 4 + 1.5 * 5.0 + 9.0),
        "chetda, chiziq ustida": (25.0, 4 + 1.5 * 25.0),
        "chetda, chiziqdan uzoq": (25.0, 4 + 1.5 * 25.0 - 18.0),
    }

    print("\n=== 2. Bitta nuqta qo'shamiz ===")
    for nom, (xn, yn) in holatlar.items():
        xs, ys = np.append(x, xn), np.append(y, yn)
        m = sm.OLS(ys, sm.add_constant(xs)).fit()
        ta = m.get_influence()
        print(f"  {nom:<24}: nishab {m.params[1]:6.3f} "
              f"({m.params[1] / nat.params[1] - 1:+6.1%}), "
              f"leverage {ta.hat_matrix_diag[-1]:.3f}, "
              f"Cook {ta.cooks_distance[0][-1]:7.3f}")

    print("\n=== 3. Chegaralar ===")
    print(f"  o'rtacha leverage = p/n = {2 / (n + 1):.4f}, "
          f"ogohlantirish chegarasi 2p/n = {4 / (n + 1):.4f}")
    print(f"  Cook chegarasi 4/n = {4 / (n + 1):.4f}")

    print("\n=== 4. Asl ma'lumotda ta'sirli nuqtalar ===")
    ta = nat.get_influence()
    cook = ta.cooks_distance[0]
    yuqori = np.argsort(cook)[::-1][:3]
    for i in yuqori:
        print(f"  i={i:>3}: x={x[i]:5.2f}, qoldiq={nat.resid[i]:+6.2f}, "
              f"leverage={ta.hat_matrix_diag[i]:.3f}, Cook={cook[i]:.4f}")
    print(f"  4/n dan katta: {(cook > 4 / n).sum()} ta")
    print("  ⭐ Yuqori leverage o'zi muammo emas — qoldiq bilan birga muammo")


if __name__ == "__main__":
    main()

=== 1. Asl model ===
  nishab = 1.476

=== 2. Bitta nuqta qo'shamiz ===
  markazda katta qoldiq   : nishab  1.473 ( -0.2%), leverage 0.010, Cook   0.179
  chetda, chiziq ustida   : nishab  1.483 ( +0.5%), leverage 0.333, Cook   0.025
  chetda, chiziqdan uzoq  : nishab  1.185 (-19.7%), leverage 0.333, Cook  14.637

=== 3. Chegaralar ===
  o'rtacha leverage = p/n = 0.0198, ogohlantirish chegarasi 2p/n = 0.0396
  Cook chegarasi 4/n = 0.0396

=== 4. Asl ma'lumotda ta'sirli nuqtalar ===
  i= 46: x= 9.70, qoldiq= +2.68, leverage=0.035, Cook=0.0932
  i= 18: x= 0.19, qoldiq= -2.40, leverage=0.042, Cook=0.0914
  i= 96: x= 0.64, qoldiq= +2.14, leverage=0.037, Cook=0.0626
  4/n dan katta: 5 ta
  ⭐ Yuqori leverage o'zi muammo emas — qoldiq bilan birga muammo

"""Bitta funksiyada barcha tekshiruvlar (real numpy/scipy/statsmodels)."""

import numpy as np
import statsmodels.api as sm
from scipy import stats
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.stattools import durbin_watson


def tekshir(y, X, nomlar) -> dict:
    """Regressiya diagnostikasi — qisqa hisobot."""
    Xc = sm.add_constant(X)
    nat = sm.OLS(y, Xc).fit()
    e, b = nat.resid, nat.fittedvalues
    chegara = np.quantile(b, [0, 0.25, 0.5, 0.75, 1.0])
    egrilik = [e[(b >= chegara[i]) & (b <= chegara[i + 1])].mean() for i in range(4)]
    return {
        "R2": nat.rsquared,
        "egrilik": max(abs(v) for v in egrilik) / np.std(e),
        "het_p": het_breuschpagan(e, Xc)[1],
        "SD_nisbat": e[b >= np.median(b)].std() / e[b < np.median(b)].std(),
        "dw": durbin_watson(e),
        "jb_p": stats.jarque_bera(e).pvalue,
        "max_VIF": max(variance_inflation_factor(Xc, i)
                       for i in range(1, Xc.shape[1])),
        "cook_yuqori": int((nat.get_influence().cooks_distance[0]
                            > 4 / len(y)).sum()),
    }


def main() -> None:
    rng = np.random.default_rng(44)
    n = 800
    maydon = rng.lognormal(np.log(70), 0.33, n).clip(28, 260)
    xona = np.clip((maydon / 24).round(), 1, 7)
    markaz = (rng.random(n) < 0.3).astype(float)
    nomlar = ["maydon", "xona", "markaz"]
    X = np.column_stack([maydon, xona, markaz])

    print("=== 1. Yaxshi model ===")
    y1 = 18 + 1.52 * maydon + 4 * xona + 38 * markaz + rng.normal(0, 13, n)
    for k, v in tekshir(y1, X, nomlar).items():
        print(f"  {k:<12}: {v:.4f}" if isinstance(v, float) else f"  {k:<12}: {v}")

    print("\n=== 2. Geteroskedastik model ===")
    y2 = 18 + 1.52 * maydon + 4 * xona + 38 * markaz + rng.normal(0, 0.12 * maydon, n)
    n2 = tekshir(y2, X, nomlar)
    for k, v in n2.items():
        print(f"  {k:<12}: {v:.4f}" if isinstance(v, float) else f"  {k:<12}: {v}")

    print("\n=== 3. Nochiziqli model ===")
    y3 = 18 + 0.011 * maydon ** 2 + 38 * markaz + rng.normal(0, 13, n)
    n3 = tekshir(y3, X, nomlar)
    for k, v in n3.items():
        print(f"  {k:<12}: {v:.4f}" if isinstance(v, float) else f"  {k:<12}: {v}")

    n1 = tekshir(y1, X, nomlar)
    print("\n=== 4. Xulosa va tavsiyalar ===")
    print(f"  yaxshi      : egrilik {n1['egrilik']:.3f}, het_p {n1['het_p']:.2f} "
          f"→ chiziqlilik va gomoskedastiklik joyida")
    print(f"  geteroskedastik: SD nisbati {n2['SD_nisbat']:.2f}, het_p "
          f"{n2['het_p']:.1e} → HC3 yoki log(y)")
    print(f"  nochiziqli  : egrilik {n3['egrilik']:.3f} (yaxshida {n1['egrilik']:.3f}) "
          f"→ x^2 belgisi 13.5-bob")
    print(f"  VIF uchala modelda ham {n1['max_VIF']:.1f} — maydon va xona bog'liq 13.3-bob;")
    print("    u faqat koeffitsiyentlarga ta'sir qiladi, bashoratga emas")
    print(f"  Cook > 4/n: {n1['cook_yuqori']} ta ({n1['cook_yuqori'] / len(y1):.1%}) — "
          f"bu chegara odatda ~5% qatorni belgilaydi, hammasi muammo emas")
    print("  ⭐ Diagnostika — bitta funksiya, har loyihada takrorlanadi")


if __name__ == "__main__":
    main()

=== 1. Yaxshi model ===
  R2          : 0.9286
  egrilik     : 0.1249
  het_p       : 0.4520
  SD_nisbat   : 1.0804
  dw          : 2.0243
  jb_p        : 0.3020
  max_VIF     : 13.4558
  cook_yuqori : 39

=== 2. Geteroskedastik model ===
  R2          : 0.9592
  egrilik     : 0.0993
  het_p       : 0.0000
  SD_nisbat   : 1.6189
  dw          : 2.0246
  jb_p        : 0.0000
  max_VIF     : 13.4558
  cook_yuqori : 42

=== 3. Nochiziqli model ===
  R2          : 0.8951
  egrilik     : 0.4414
  het_p       : 0.0000
  SD_nisbat   : 1.3134
  dw          : 2.0474
  jb_p        : 0.0000
  max_VIF     : 13.4558
  cook_yuqori : 29

=== 4. Xulosa va tavsiyalar ===
  yaxshi      : egrilik 0.125, het_p 0.45 → chiziqlilik va gomoskedastiklik joyida
  geteroskedastik: SD nisbati 1.62, het_p 6.0e-24 → HC3 yoki log(y)
  nochiziqli  : egrilik 0.441 (yaxshida 0.125) → x^2 belgisi 13.5-bob
  VIF uchala modelda ham 13.5 — maydon va xona bog'liq 13.3-bob;
    u faqat koeffitsiyentlarga ta'sir qiladi, bashoratga emas
  Cook > 4/n: 39 ta (4.9%) — bu chegara odatda ~5% qatorni belgilaydi, hammasi muammo emas
  ⭐ Diagnostika — bitta funksiya, har loyihada takrorlanadi

stats.shapiro(qoldiq)                                             # ⚠️
# avval qoldiq ~ bashorat grafigi (chiziqlilik)                   # ✅

nat = sm.OLS(y, X).fit()               # tor CI                   # ⚠️
nat = sm.OLS(y, X).fit(cov_type="HC3")                            # ✅

df = df[np.abs(qoldiq) < 3 * qoldiq.std()]                        # ⚠️
# tekshiring: xatomi yoki haqiqiy kuzatuvmi                       # ✅

if het_p < 0.05: print("muammo")       # n = 100000               # ⚠️
print(f"SD nisbati {nisbat:.2f}")      # effekt kattaligi         # ✅

sm.OLS(y, X).fit()                     # kunlik ma'lumot          # ⚠️
durbin_watson(nat.resid); cov_type="HAC"                          # ✅

# test qoldiqlarini tahlil qilish                                 # ⚠️
# o'quvda diagnostika, testda faqat yakuniy baho 12.3-bob           # ✅

# log(y) modeli: "koeffitsiyent 0.12 ming so'm"                   # ⚠️
# "+12.7%" 13.3-bob                                                 # ✅

1.  # beshta taxmin?
2.  # eng muhimi?
3.  # geteroskedastiklik nimani buzadi?
4.  # eng oson yechimi?
5.  # normallik qachon muhim?
6.  # leverage nima?
7.  # Cook masofasi nima?
8.  # Cook chegarasi?
9.  # Durbin-Watson qiymati?
10. # markaziy diagnostika grafigi?
11. # katta namunada testlar nega aldaydi?
12. # diagnostika qaysi to'plamda?

1.  stats.shapiro(resid)   # yagona tekshiruv

2.  nat = sm.OLS(y, X).fit()   # konus shaklli qoldiqlar

3.  df = df[np.abs(z) < 3]   # tekshirmasdan

4.  if bp_p < 0.05: print("muammo")   # n = 200000

5.  sm.OLS(y, X).fit()   # kunlik vaqt qatori