Tutorial – Session 7¶
Intervalestimation og konfidensintervaller¶
Denne tutorial følger Ross, kapitel 7, afsnit 7.1, 7.3, 7.4 og 7.5 (jf. sessionbeskrivelsen). For hvert emne gives først de centrale formler og fortolkninger, derefter Python-eksempler med numpy, scipy.stats og matplotlib.
Notebook-version: Tutorial_7_notebook.ipynb
Ekstra plot-notebook: plot_normal_distribution.ipynb
1. Introduktion og opsætning¶
Vi arbejder med konfidensintervaller for middelværdi, varians, forskel mellem middelværdier og andele, samt prædiktionsintervaller for en ny observation. Biblioteker:
numpy– tal og stikprøvestatistikkerscipy.stats–norm,t,chi2(og sandsynligheder)matplotlib.pyplot– plots (fx normalfordelingsplot)
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, t, chi2, probplot
np.random.seed(42)
plt.rcParams["figure.figsize"] = (10, 6)
plt.style.use("ggplot")
2. Normal probability plot (QQ-plot) og normalitetscheck¶
Flere af intervallerne i denne session (særligt klassiske \(t\)-intervaller, prædiktionsintervaller i normalmodellen og \(\chi^2\)-intervaller for varians) bygger på antagelser om normalfordeling. Derfor er det god praksis at starte med et normal probability plot (QQ-plot).
Hvorfor tjekke normalitet?¶
- For små og mellemstore stikprøver kan ikke-normalitet påvirke intervalkvaliteten markant.
- Et QQ-plot giver et hurtigt visuelt check af, om data følger en omtrent lineær normalstruktur.
- Tydelige afvigelser i halerne eller enkeltstående outliers er et signal om, at man bør være forsigtig med standardformlerne.
Hvordan fortolker man QQ-plot?¶
- Punkter tæt på en ret linje: data er omtrent normalfordelte.
- Systematisk krumning (S-form): skævhed eller tunge/lette haler.
- Få punkter langt fra linjen: mulige outliers.
Du kan også se en separat og mere visuelt orienteret notebook her: plot_normal_distribution.ipynb.
# Eksempel: normal probability plot (QQ-plot)
sample_data = np.array([3.2, 3.5, 3.6, 3.8, 4.0, 4.1, 4.2, 4.5, 4.8, 5.4])
fig, ax = plt.subplots()
probplot(sample_data, dist="norm", plot=ax)
ax.set_title("QQ-plot: hurtig normalitetsvurdering")
plt.tight_layout()
plt.show()
3. Indledning til konfidensintervaller for middelværdi¶
Et \(100(1-\alpha)\,\%\) konfidensinterval for en parameter er konstrueret sådan, at hvis vi gentager forsøget (ny stikprøve, samme procedure), så vil cirka \(100(1-\alpha)\,\%\) af de således beregnede intervaller indeholde den sande parameter.
Vigtigt: Det er ikke korrekt at sige, at “sandsynligheden for at \(\mu\) ligger i dette interval er \(95\,\%\)” efter data er observeret — intervallet er fast, og \(\mu\) er ikke en tilfældig variabel i den frekventistiske fortolkning.
I de næste to underafsnit gennemgår vi CI for middelværdi i normalmodellen for hhv. kendt og ukendt populationsstandardafvigelse.
3.1 Konfidensinterval for \(\mu\) når \(\sigma\) er kendt (normalpopulation)¶
Antag \(X_1,\ldots,X_n\) uafhængige \(N(\mu,\sigma^2)\) med kendt \(\sigma\). Så er $\(Z = \frac{\bar{X}-\mu}{\sigma/\sqrt{n}} \sim N(0,1).\)$
Tosidet \(100(1-\alpha)\,\%\)-CI: $\(\bar{x} \pm z_{\alpha/2}\,\frac{\sigma}{\sqrt{n}}.\)$
Ensidet nedre grænse (\(\mu\) er mindst … med konfidens \(1-\alpha\)): $\(\bar{x} - z_{\alpha}\,\frac{\sigma}{\sqrt{n}} \leq \mu.\)$
Ensidet øvre analogt med \(+\) og \(z_\alpha\).
# 1-til-1 implementering af formlen (uden helper-funktion)
xbar = 100.0
sigma = 15.0
n = 36
confidence = 0.95
alpha = 1 - confidence
se = sigma / np.sqrt(n)
z_crit = norm.ppf(1 - alpha / 2)
lo = xbar - z_crit * se
hi = xbar + z_crit * se
print(f"95% tosidet z-CI: [{lo:.4f}, {hi:.4f}]")
# Ensidet nedre grænse:
z_alpha = norm.ppf(1 - alpha)
lower_bound = xbar - z_alpha * se
print(f"95% ensidet nedre grænse: mu >= {lower_bound:.4f}")
# (valgfrit) helper-funktion til genbrug:
def ci_mean_sigma_known(xbar, sigma, n, confidence=0.95, side="two"):
alpha = 1 - confidence
se = sigma / np.sqrt(n)
if side == "two":
z = norm.ppf(1 - alpha / 2)
return xbar - z * se, xbar + z * se
if side == "lower":
z = norm.ppf(1 - alpha)
return xbar - z * se, np.inf
if side == "upper":
z = norm.ppf(1 - alpha)
return -np.inf, xbar + z * se
raise ValueError("ukendt side")
3.2 Konfidensinterval for \(\mu\) når \(\sigma\) er ukendt (\(t\)-fordeling)¶
Erstat \(\sigma\) med \(S\) (stikprøvestandardafvigelse). Pivot: $\(T = \frac{\bar{X}-\mu}{S/\sqrt{n}} \sim t_{n-1}.\)$
Tosidet CI: $\(\bar{x} \pm t_{\alpha/2,\,n-1}\,\frac{s}{\sqrt{n}}.\)$
# 1-til-1 implementering af t-intervallet (uden helper-funktion)
x = np.array([9.8, 10.2, 10.1, 9.9, 10.5, 9.7, 10.0, 10.3, 9.8, 10.1])
confidence = 0.95
n = len(x)
xbar = np.mean(x)
s = np.std(x, ddof=1)
df = n - 1
alpha = 1 - confidence
se = s / np.sqrt(n)
t_crit = t.ppf(1 - alpha / 2, df)
lo = xbar - t_crit * se
hi = xbar + t_crit * se
print(f"95% tosidet t-CI: [{lo:.4f}, {hi:.4f}]")
# Ensidet nedre grænse:
t_alpha = t.ppf(1 - alpha, df)
lower_bound = xbar - t_alpha * se
print(f"95% ensidet nedre grænse: mu >= {lower_bound:.4f}")
# (valgfrit) helper-funktion til genbrug:
def ci_mean_sigma_unknown(x, confidence=0.95, side="two"):
x = np.asarray(x, dtype=float)
n = len(x)
xbar = np.mean(x)
s = np.std(x, ddof=1)
df = n - 1
alpha = 1 - confidence
se = s / np.sqrt(n)
if side == "two":
tcrit = t.ppf(1 - alpha / 2, df)
return xbar - tcrit * se, xbar + tcrit * se
if side == "lower":
tcrit = t.ppf(1 - alpha, df)
return xbar - tcrit * se, np.inf
if side == "upper":
tcrit = t.ppf(1 - alpha, df)
return -np.inf, xbar + tcrit * se
raise ValueError("side skal være 'two', 'lower' eller 'upper'")
4. Prædiktionsinterval for en ny observation (Ross 7.3.2)¶
Målet er én ny observation \(X_{n+1}\), ikke \(\mu\).
- \(\sigma\) kendt: med pivot der udnytter uafhængighed mellem \(\bar{X}\) og den nye observation fås typisk $\(\bar{x} \pm z_{\alpha/2}\,\sigma\sqrt{1+\frac{1}{n}}.\)$
- \(\sigma\) ukendt (erstattes med \(s\)): $\(\bar{x} \pm t_{\alpha/2,\,n-1}\,s\sqrt{1+\frac{1}{n}}.\)$
Bemærk faktoren \(\sqrt{1+1/n}\) i stedet for kun \(\sigma/\sqrt{n}\) — der er både usikkerhed på \(\mu\) og spredning på den nye måling.
def prediction_interval(x, confidence=0.95, sigma=None, side="two"):
x = np.asarray(x, dtype=float)
n = len(x)
xbar = np.mean(x)
s = np.std(x, ddof=1)
df = n - 1
alpha = 1 - confidence
if sigma is not None:
if side == "two":
z = norm.ppf(1 - alpha / 2)
half = z * sigma * np.sqrt(1 + 1 / n)
else:
raise NotImplementedError("brug tosidet eller udvid med z_alpha")
return xbar - half, xbar + half
if side == "two":
tcrit = t.ppf(1 - alpha / 2, df)
half = tcrit * s * np.sqrt(1 + 1 / n)
return xbar - half, xbar + half
raise NotImplementedError
5. Konfidensintervaller for \(\sigma^2\) og \(\sigma\) (\(\chi^2\), normalpopulation)¶
For normalfordelte data gælder $\(\frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1}.\)$
Tosidet \(100(1-\alpha)\,\%\) CI for \(\sigma^2\): $\(\left(\frac{(n-1)s^2}{\chi^2_{\alpha/2,\,n-1}},\;\frac{(n-1)s^2}{\chi^2_{1-\alpha/2,\,n-1}}\right).\)$
Ensidet øvre grænse for \(\sigma^2\) (ønsker “\(\sigma^2\) er højst …”): $\(\sigma^2 \leq \frac{(n-1)s^2}{\chi^2_{1-\alpha,\,n-1}}.\)$ Tag kvadratrod for et øvre grænse-udtryk for \(\sigma\).
def ci_var_sigma_sq(x, confidence=0.95, bound="two-sided"):
x = np.asarray(x, dtype=float)
n = len(x)
s2 = np.var(x, ddof=1)
nu = n - 1
alpha = 1 - confidence
if bound == "two-sided":
chi_lo = chi2.ppf(alpha / 2, nu)
chi_hi = chi2.ppf(1 - alpha / 2, nu)
return nu * s2 / chi_hi, nu * s2 / chi_lo
if bound == "upper-var": # øvre grænse for sigma^2 (ensidet)
chi_c = chi2.ppf(1 - alpha, nu)
return 0.0, nu * s2 / chi_c
raise ValueError
6. To uafhængige stikprøver: \(\mu_1 - \mu_2\) (Ross 7.4)¶
Typisk setup i lærebøger: to uafhængige stikprøver fra normalfordelinger; ofte ukendte men ens varianser \(\sigma^2\). Så bruges pooled spredning $\(S_p^2 = \frac{(n_1-1)S_1^2 + (n_2-1)S_2^2}{n_1+n_2-2}\)$ og $\((\bar{X}_1-\bar{X}_2) \pm t_{\alpha/2,\,n_1+n_2-2}\; S_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}.\)$
Hvis \(\sigma_1^2\) og \(\sigma_2^2\) er kendte, bruges \(z\) med standardfejl \(\sqrt{\sigma_1^2/n_1+\sigma_2^2/n_2}\).
Tjek altid de præcise antagelser i din Ross-udgave (ens varianser m.m.).
def ci_diff_means_pooled(x1, x2, confidence=0.95):
x1, x2 = np.asarray(x1, float), np.asarray(x2, float)
n1, n2 = len(x1), len(x2)
m1, m2 = np.mean(x1), np.mean(x2)
v1, v2 = np.var(x1, ddof=1), np.var(x2, ddof=1)
df = n1 + n2 - 2
sp2 = ((n1 - 1) * v1 + (n2 - 1) * v2) / df
sp = np.sqrt(sp2)
se = sp * np.sqrt(1 / n1 + 1 / n2)
alpha = 1 - confidence
tcrit = t.ppf(1 - alpha / 2, df)
diff = m1 - m2
return diff - tcrit * se, diff + tcrit * se
# Illustration med simulerede data
a = np.random.normal(5.0, 1.0, 12)
b = np.random.normal(6.5, 1.0, 15)
print("95% CI for mu1 - mu2:", ci_diff_means_pooled(a, b))
7. Konfidensinterval for populationsandel \(p\) (Ross 7.5)¶
Med \(\hat{p}=X/n\) (andel succeser) og stor \(n\) bruges ofte Wald-intervallet $\(\hat{p} \pm z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}.\)$ Ensidet nedre grænse: \(\hat{p} - z_{\alpha}\sqrt{\hat{p}(1-\hat{p})/n} \leq p\).
For små \(n\) eller ekstremer \(\hat{p}\) kan normaltilnærmelsen være dårlig; et alternativ er Wilson-intervallet (se evt. statsmodels.stats.proportion.proportion_confint hvis du har pakken installeret).
def ci_proportion_wald(x_successes, n, confidence=0.95, side="two"):
p_hat = x_successes / n
alpha = 1 - confidence
se = np.sqrt(p_hat * (1 - p_hat) / n)
if side == "two":
z = norm.ppf(1 - alpha / 2)
return p_hat - z * se, p_hat + z * se
if side == "lower":
z = norm.ppf(1 - alpha)
return p_hat - z * se, 1.0
if side == "upper":
z = norm.ppf(1 - alpha)
return 0.0, p_hat + z * se
raise ValueError
8. Gennemgang af øvelserne (samme tal som i README)¶
Øvelse 1 – stempelringe (\(\sigma\) kendt, \(99\,\%\))
xbar, sigma, n = 74.036, 0.001, 15
lo, hi = ci_mean_sigma_known(xbar, sigma, n, confidence=0.99, side="two")
print(f"99% CI (tosidet): [{lo:.4f}, {hi:.4f}]")
lo1, _ = ci_mean_sigma_known(xbar, sigma, n, confidence=0.99, side="lower")
print(f"99% nedre grænse (ensidet): mu >= {lo1:.4f}")
Øvelse 2 – beton (\(\sigma^2=1000\), \(n=12\))
xbar, n = 3250, 12
sigma = np.sqrt(1000)
for conf in (0.95, 0.99):
lo, hi = ci_mean_sigma_known(xbar, sigma, n, confidence=conf, side="two")
print(f"{int(conf*100)}% CI: [{lo:.2f}, {hi:.2f}] (bredde {hi-lo:.2f})")
Øvelse 3 – CNN speedup (\(t\)-interval)
cnn = np.array([
3.775302, 3.350679, 4.217981, 4.030324, 4.639692,
4.139665, 4.395575, 4.824257, 4.268119, 4.584193,
4.930027, 4.315973, 4.600101,
])
n = len(cnn)
xbar = np.mean(cnn)
s = np.std(cnn, ddof=1)
print(f"n={n}, xbar={xbar:.3f}, s={s:.4f}")
fig, ax = plt.subplots()
probplot(cnn, dist="norm", plot=ax)
ax.set_title("Normalfordelingsplot (CNN speedup)")
plt.tight_layout()
plt.show()
print("95% CI:", ci_mean_sigma_unknown(cnn, 0.95, "two"))
print("95% nedre grænse:", ci_mean_sigma_unknown(cnn, 0.95, "lower"))
Øvelse 4 - papirvægt (\(\chi^2\), ensidet øvre for \(\sigma\))
paper = np.array([
3.481, 3.448, 3.485, 3.475, 3.472,
3.477, 3.472, 3.464, 3.472, 3.470,
3.470, 3.470, 3.477, 3.473, 3.474,
])
n = len(paper)
s = np.std(paper, ddof=1)
nu = n - 1
alpha = 0.05
chi_crit = chi2.ppf(1 - alpha, nu) # = chi^2_{0.95,14} i facit-noteringen
upper_sigma_sq = nu * s**2 / chi_crit
print(f"s={s:.5f}, øvre 95% grænse for sigma: {np.sqrt(upper_sigma_sq):.4f}")
Øvelse 5 - Ohio-andel
x, n = 412, 768
print("95% CI:", ci_proportion_wald(x, n, 0.95, "two"))
print("95% nedre grænse:", ci_proportion_wald(x, n, 0.95, "lower"))
9. Opsummering¶
| Situation | Fordeling / metode |
|---|---|
| \(\mu\), \(\sigma\) kendt | \(z\) |
| \(\mu\), \(\sigma\) ukendt | \(t_{n-1}\) |
| Ny observation | Prædiktion: \(\sqrt{1+1/n}\)-faktor |
| \(\sigma^2\) / \(\sigma\) | \(\chi^2_{n-1}\) |
| \(\mu_1-\mu_2\), ens ukendt varians | Pooled \(t_{n_1+n_2-2}\) |
| Andel \(p\), stor \(n\) | Wald (\(z\)) |
Du skal kunne: vælge rigtig fordeling, skelne CI for \(\mu\) fra prædiktionsinterval, og implementere beregningerne i Python.