Διαστήματα εμπιστοσύνης για μέσο κανονικής με άγνωστη διασπορά

Ας υποθέσουμε ότι η \(X\) ακολουθεί κανονική \(N(\mu,\sigma^2)\) και θέλουμε να παράγουμε διαστήματα εμπιστοσύνης για το \(\mu\) χωρίς να ξέρουμε το \(\sigma^2\), όπως το ξέραμε στην προηγούμενη διάλεξη.

Μια λύση σε αυτό το πρόβλημα, που είναι εφικτή όταν το δείγμα μας είναι μεγάλο, είναι να κάνουμε πρώτα μια εκτίμηση για το \(\sigma^2\) και να χρησιμοποιήσουμε αυτή την εκτίμηση στη θέση του πραγματικού \(\sigma^2\) στην ανάλυση που κάναμε στην προηγούμενη διάλεξη, να προσποιηθούμε με άλλα λόγια ότι γνωρίζουμε το ποιο είναι το \(\sigma^2\) και ότι η εκτίμησή μας είναι ακριβής. Κάτι τέτοιο έχει πιθανότητες να δουλέψει μόνο όταν η εκτίμησή μας για το \(\sigma\) είναι εγγυημένα ακριβής (στα πλαίσια πάντα της θεωρίας πιθανοτήτων όπου εγγυημένο σημαίνει ότι με μεγάλη πιθανότητα είναι σωστό), και όταν το δείγμα μας είναι μικρό αυτό δε μπορεί να συμβεί.

Ας συμβολίσουμε και πάλι με \(\overline{X}_N = \frac{X_1+\cdots+X_N}{N}\) το δειγματικό μέσο του δείγματος και με \(S^2 = \frac{1}{N-1} \sum_{i=1}^N (X_i-\overline{X}_N)^2\) τη δειγματική διασπορά.

Θυμόμαστε από την προηγούμενη διάλεξη ότι και η \(\overline{X}_N\) είναι κανονική και μάλιστα με μέσο \(\mu\) και διασπορά \(\sigma^2/N\).

Το διάστημα εμπιστοσύνης που θα κατασκευάσουμε θα είναι της μορφής \[ (\overline{X}_N - L, \overline{X}_N+L), \] όπου με \(L\) συμβολίζουμε το μισό μήκος του διαστήματος. Το να μην ανήκει ο μέσος \(\mu\) στο διάστημα αυτό είναι ισοδύναμο με την ανισότητα \[ {{\left|{\overline{X}_N - \mu}\right|}} > L, \] και αυτό θέλουμε συμβαίνει με πιθανότητα το πολύ \(p\), όπου \(1-p\) είναι το επίπεδο εμπιστοσύνης του διαστήματος που θα φτιάξουμε. Η ΤΜ \[ \frac{{{\overline{X}_N}}-\mu}{\sigma/\sqrt{N}} \] ακολουθεί τη \(N(0,1)\), όπως είδαμε στην προηγούμενη διάλεξη, πλην όμως τώρα δε γνωρίζουμε το \(\sigma\) για να μπορέσουμε να χρησιμοποιήσουμε τη μεταβλητή. Αυτό που μπορούμε όμως να κάνουμε είναι να την αντικαταστήσουμε στην ανάλυσή μας με την ΤΜ \[ Y = \frac{{{\overline{X}_N}}-\mu}{S/\sqrt{N}}, \] όπου \(S = \sqrt{S^2}\) είναι η τετραγωνική ρίζα της δειγματικής διασποράς. Ας προσποιηθούμε προς το παρόν ότι γνωρίζουμε την κατανομή της \(Y\) και ας ονομάσουμε \(F(\cdot)\) τη συνάρτηση κατανομής της. Ας είναι επίσης \(y_0 > 0\) μια τιμή για την οποία \[ F(-y_0) \le \frac{p}{2}, \ \ \ F(y_0) \ge 1- \frac{p}{2} \] και μάλιστα επιλέγουμε το \(y_0\) όσο πιο μικρό γίνεται κατά τα άλλα.

Υποθέτουμε επίσης ότι το \(L\) και το \(N\) είναι τόσο μεγάλα ώστε να ισχύει \[ \frac{L}{S/\sqrt{N}} \ge y_0. \] Τότε λόγω αυτής της ανισότητας θα έχουμε το επιθυμητό αποτέλεσμα: \[ \begin{align*} {{\mathbb P}\left[{{{\left|{{{\overline{X}_N}}-\mu}\right|}} > L}\right]} &= {{\mathbb P}\left[{{{\overline{X}_N}}-\mu < -L}\right]}+{{\mathbb P}\left[{{{\overline{X}_N}}-\mu > L}\right]}\\ &= {{\mathbb P}\left[{\frac{{{\overline{X}_N}}-\mu}{S/\sqrt{N}} < -\frac{L}{S/\sqrt{N}}}\right]} + {{\mathbb P}\left[{\frac{{{\overline{X}_N}}-\mu}{S/\sqrt{N}} > \frac{L}{S/\sqrt{N}}}\right]}\\ &\le {{\mathbb P}\left[{\frac{{{\overline{X}_N}}-\mu}{S/\sqrt{N}} < -y_0}\right]} + {{\mathbb P}\left[{\frac{{{\overline{X}_N}}-\mu}{S/\sqrt{N}} > y_0}\right]}\\ &= F(-y_0) + 1-F(y_0)\\ &\le \frac{p}{2}+\frac{p}{2}\\ &= p. \end{align*} \] Συνοψίζοντας, η διαδικασία για να παράγουμε ένα διάστημα εμπιστοσύνης για το \(\mu\) είναι:

  1. Να βρούμε τη μικρότερη τιμή \(y_0\) για την οποία ισχύουν οι ανισότητες \[ F(-y_0) \le \frac{p}{2}, \ \ \ F(y_0) \ge 1- \frac{p}{2}. \]

  2. Να βεβαιωθούμε ότι το \(L\) και το \(N\) είναι τέτοια ώστε να ισχύει η ανισότητα \[ \frac{L}{S/\sqrt{N}} \ge y_0. \]

Η κατανομή \(F\), της ΤΜ \(\frac{{{\overline{X}_N}}-\mu}{S/\sqrt{N}}\), υπολογίστηκε από τον W. Gosset (δείτε εδώ την ενδιαφέρουσα ιστορία) και ονομάζεται κατανομή \(t\), ή κατανομή \(t\) του Student (Student’s t-distribution).

Πρόκειται για μια κατανομή που μοιάζει με την τυπική κανονική και μάλιστα συγκλίνει σε αυτήν για \(N \to \infty\) (η κατανομή \(t\) έχει το \(N\) ως παράμετρο). Η πυκνότητά της δίνεται από τον τύπο \[ f_N(t) = \frac{\Gamma(N/2)}{\sqrt{\pi(N-1)}\Gamma((N-1)/2)} (1+\frac{t^2}{N-1})^{-N/2}. \] Εδώ η συνάρτηση \(\Gamma(\cdot)\) ορίζεται (για θετικά \(t\)) από το ολοκλήρωμα \[ \Gamma(t) = \int_0^\infty x^{t-1}e^{-x}\,dx, \] και ικανοποιεί την \(\Gamma(k) = (k-1)!\) για κάθε θετικό ακέραιο \(k\), ιδιότητα από την οποία μπορεί κανείς να υπολογίσει τη σταθερά που εμφανίζεται στην πυκνότητα \(f_N\) παραπάνω χωρίς χρήση ολοκληρωμάτων (αλλά προκύπτει διαφορετικός τύπος για άρτια και για περιττά \(N\)).

Από τον παραπάνω τύπο για την πυκνότητα \(t\) γίνεται φανερό ότι πρόκειται για άρτια συνάρτηση. Η ποσότητα \(N-1\) ονομάζεται “βαθμοί ελευθερίας” της κατανομής \(t\) και πρέπει να είναι τουλάχιστον 1.

Δείχνουμε πώς μπορούμε να σχεδιάσουμε μερικές από αυτές τις πυκνότητες, για διάφορες τιμές βαθμών ελευθερίας, στη γλώσσα R. Παρατηρείστε πως για μεγάλους βαθμούς ελευθερίας η κατανομή \(t\) σχεδόν συμπίπτει με την τυπική κανονική. Στη γλώσσα R η πυκνότητα της t δίνεται από τη συνάρτηση dt.

x = seq(-4, 4, length=100) # τιμές του x για τα οποία υπολογίζεται η πυκνότητα
hx <- dnorm(x) # η τυπική κανονική πυκνότητα στα σημεία x

degf <- c(1, 3, 8, 30) # λίστα με τους βαθμούς ελευθερίας για τους οποίους υπολογίζεται η πυκνότητα t
colors <- c("red", "blue", "darkgreen", "gold", "black") # διαφορετικά χρώματα για διαφορετικούς βαθμούς ελευθερίας
labels <- c("df=1", "df=3", "df=8", "df=30", "N(0,1)") # κείμενο στο γράφημα για διαφορετικούς βαθμούς ελευθερίας

plot(x, hx, type="l", lty=2, xlab="x", # σχεδιάζουμε την τυπική κανονική
     ylab="Density", main="t-distributions")

for (i in 1:4){ # σχεδιάζουμε την t για κάθε βαθμό ελευθερίας
  lines(x, dt(x,degf[i]), lwd=2, col=colors[i])
}

legend("topright", inset=.05, title="Distributions", # τοποθετούμενο επεξηγηματικό κείμενο στο γράφημα
       labels, lwd=2, lty=c(1, 1, 1, 1, 2), col=colors)

plot of chunk unnamed-chunk-2

Εφαρμογή Δίνουμε τώρα μια εφαρμογή της παραπάνω μεθόδου πάνω σε κάποιες μετρήσεις της ταχύτητας του φωτός (δείτε εδώ για την ιστορία).Τα δεδομένα έχουν τοποθετηθεί στο αρχείο sl.txt και τα εισάγουμε με την εντολή

sl = read.table("sl.txt", header=TRUE)
N = length(sl$time)
##    time
## 1    28
## 2    22
## 3    36
## 4    26
## 5    28
## 6    28
## 7    26
## 8    24
## 9    32
## 10   30
## 11   27
## 12   24
## 13   33
## 14   21
## 15   36
## 16   32
## 17   31
## 18   25
## 19   24
## 20   25
## 21   28
## 22   36
## 23   27
## 24   32
## 25   34
## 26   30
## 27   25
## 28   26
## 29   26
## 30   25
## 31   23
## 32   21
## 33   30
## 34   33
## 35   29
## 36   27
## 37   29
## 38   28
## 39   22
## 40   26
## 41   27
## 42   16
## 43   31
## 44   29
## 45   36
## 46   32
## 47   28
## 48   40
## 49   19
## 50   37
## 51   23
## 52   32
## 53   29
## 54   24
## 55   25
## 56   27
## 57   24
## 58   16
## 59   29
## 60   20
## 61   28
## 62   27
## 63   39
## 64   23

Πρόκειται για N=64 μετρήσεις του χρόνου που έκανε το φως να ταξιδέψει μια σταθερή απόσταση 7.44373 km, μετρημένες σε μsec. Υποθέτουμε ότι οι μετρήσεις αυτές προέρχονται από μια κανονική κατανομή με μέσο \(\mu\). Σκοπός μας είναι να δώσουμε ένα διάστημα εμπιστοσύνης 95% για το μέσο χρόνο αυτό \(\mu\). Η παράμετρος \(p\) λοιπόν είναι δεδομένη και ίση με \(0.05\), όπως δεδομένη είναι και η παράμετρος \(N\):

p = 0.05
N = length(sl$time)

Απομένει να προσδιορίσουμε την παράμετρο \(L\), το μισό μήκος του διαστήματος εμπιστοσύνης του οποίου το μέσο είναι ο δειγματικός μέσος:

samplemean = mean(sl$time)
samplemean
## [1] 27.75

Υπολογίζουμε επίσης τη δειγματική τυπική απόκλιση \(S\):

S = sqrt(var(sl$time))
S
## [1] 5.083

Αρχίζουμε την ανάλυσή μας βρίσκοντας την τιμή \(y_0\) χρησιμοποιώντας τη συνάρτηση qt της R, η οποία είναι η αντίστροφη συνάρτηση της συνάρτησης κατανομής \(t\). Η qt παίρνει μόνο μια παράμετρο, την df (βαθμοί ελευθερίας), την οποία τη θέτουμε ίση με \(N-1\) = 63.

y0 = -qt(p/2, df=N-1)
y0
## [1] 1.998

Έχοντας βρει το \(y_0\) ορίζουμε το \(L\) από την εξίσωση \[ L = \frac{y_0 S}{\sqrt{N}} \]

L = y0*S/sqrt(N)
L
## [1] 1.27

Προκύπτει η τιμή (@) \(L=1.2698\), άρα το διάστημα εμπιστοσύνης 95% που φτιάξαμε είναι το \[ (27.75-1.2698, 27.75+1.2698) = (26.4802, 29.0198). \]