Άσκηση 1

Γράψτε μια συνάρτηση που παίρνει ως όρισμα ένα πίνακα και επιστρέφει τις συντεταγμένες (αριθμός γραμμής και αριθμός στήλης) του μέγιστου στοιχείου και του ελάχιστου στοιχείου του πίνακα, Αν υπάρχουν περισσότερα από ένα μέγιστα ή ελάχιστα στοιχεία θα επιστρέφει τις συντεταγμένες όλων. (Υπόδειξη: Δείτε τη συνάρτηση which)

Απάντηση

Η συνάρτηση which παίρνει όρισμα μια λογική συνθήκη που περιλαμβάνει ένα διάνυσμα ή πίνακα, και επιστρέφει τις συντεταγμένες των στοιχείων που ικανοποιούν τη συνθήκη. Για παράδειγμα αν \(x\) είναι το διάνυσμα των 10 πρώτων άρτιων θετικών ακεραίων

x=seq(2,20,by=2)
x
##  [1]  2  4  6  8 10 12 14 16 18 20

η εντολή \(x>3\) δίνει για κάθε στοιχείο TRUE ή FALSE ανάλογα με το αν το αντίστοιχο στοιχείο ικανοποιεί τη συνθήκη \(>3\)

x>8
##  [1] FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

ενώ η εντολή

which(x>8)
## [1]  5  6  7  8  9 10

δίνει τις συντεταγμένες των στοιχείων που είναι μεγαλύτερα του 8, δηλαδή από το 5ο μέχρι και το 10ο.

Στο πρόβλημα που εξετάζουμε η λογική συνθήκη για να βρεθούν τα στοιχεία του \(x\) που είναι ίσα με το μέγιστο στοιχείο είναι \(x==max(x)\). Επομένως ορίζουμε τη συνάρτηση

wminmax=function(x)
{
  xmax=max(x)
  xmin=min(x)
  wmin=which(x==xmin,arr.ind=T)
  wmax=which(x==xmax,arr.ind=T)
  return(list(wmin, wmax))
}

Η εντολή which καλείται με την παράμετρο arr.ind=T, έτσι ώστε αν το x είναι πίνακας, οι συντεταγμένες να δίνονται με τη μορφή ζεύγους (αρ. γραμμής, αρ. στήλης). Διαφορετικά δίνονται ως ακέραιοι αριθμοί (π.χ. σε ένα πίνακα 4Χ5 το \(x[6]\) αντιστοιχεί στο 6ο στοιχείο του πίνακα πηγαίνοντας κατά στήλη, δηλαδή στο στοιχείο \(x[2,2]\)). Το αποτέλεσμα της συνάρτησης επιστρέφει με τη μορφή λίστας, επειδή τα wmin και wmax μπορεί να είναι διανύσματα ή πίνακες που δεν έχουν τον ίδιο αριθμό στοιχείων. Καλώντας τη συνάρτηση για το διάνυσμα \(x\) που ορίσαμε παραπάνω

l=wminmax(x)
l
## [[1]]
## [1] 1
## 
## [[2]]
## [1] 10

ενώ για τον πίνακα

y=matrix(c(3,2,2,3,5,4,3,2,3,5,6,3,2,6,4,3,2,2,2,3),4,5)
y
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    3    5    3    2    2
## [2,]    2    4    5    6    2
## [3,]    2    3    6    4    2
## [4,]    3    2    3    3    3
wminmax(y)
## [[1]]
##      row col
## [1,]   2   1
## [2,]   3   1
## [3,]   4   2
## [4,]   1   4
## [5,]   1   5
## [6,]   2   5
## [7,]   3   5
## 
## [[2]]
##      row col
## [1,]   3   3
## [2,]   2   4

το αποτέλεσμα είναι μια λίστα της οποίας το πρώτο στοιχείο είναι ένας πίνακας με τις συντεταγμένες των στοιχείων του \(y\) που είναι ίσα με 2 (ελάχιστο), και το δεύτερο των στοιχείων που είναι ίσα με 6 (μέγιστο).

Άσκηση 2

Γράψτε μια συνάρτηση που παίρνει ως όρισμα ένα data.frame \(A\), και σε κάθε στήλη τύπου numeric βρίσκει τις ελλείπουσες τιμές (NA) και τις αντικαθιστά με τη μέση τιμή των υπόλοιπων στοιχείων της στήλης. Επιστρέφει ως αποτέλεσμα ένα νέο data.frame με τις συμπληρωμένες τιμές ενώ το αρχικό παραμένει ως έχει.

Απάντηση

Πρώτα βλέπουμε κάποια θέματα σχετικά με NA και NaN στοιχεία σε ένα διάνυσμα, πίνακα ή dataframe. Το ΝΑ σημαίνει ότι το αντίστοιχο στοιχείο δεν έχει οριστεί, δηλαδή είναι ελλείπουσα τιμή (missing value). Το NaN προκύπτει όταν το αποτέλεσμα μιας μαθηματικής πράξης είναι απροσδιόριστο, π.χ. 0/0:

x=0/0
x
## [1] NaN

Η συνάρτηση is.na(x) είναι true αν το στοιχείο είναι είτε NA είτε NaN, ενώ η is.nan(x) είναι true αν το στοιχείο είναι NaN.

Έστω ένα διάνυσμα που περιέχει στοιχεία NA ή NaN. Αν υπολογίσουμε τη μέση τιμή, θα πάρουμε αποτέλεσμα NA ή ΝaN, γιατί η συνάρτηση προσπαθεί να υπολογίσει τη μέση τιμή όλων των στοιχείων. Για να υπολογίσουμε τη μέση τιμή παραλείποντας τα NA/NaN στοιχεία, πρέπει να περιλάβουμε την παράμετρο na.rm=TRUE μέσα στη mean:

x=c(1,1,2,NA,3,5)
mean(x)
## [1] NA
mean(x,na.rm=TRUE)
## [1] 2.4

Με την na.rm=TRUE, η μέση τιμή υπολογίζεται μόνο για τα στοιχεία που έχουν τιμή, που στο παραπάνω παράδειγμα είναι 5 στο πλήθος και όχι 6 που είναι το μήκος του \(x\).

Επομένως για να αντικαταστήσουμε τα στοιχεία NA/NaN ενός διανύσματος με τη μέση τιμή των υπόλοιπων, η εντολή είναι

x[is.na(x)]=mean(x,na.rm=TRUE)
x
## [1] 1.0 1.0 2.0 2.4 3.0 5.0

Με βάση τα παραπάνω για να αντικαταστήσουμε τα NA/NaN στοιχεία ενός dataframe με τους μέσους των αντίστοιχων στηλών, πρέπει πρώτα να ελέγξουμε ποιες στήλες έχουν αριθμητικά στοιχεία με την εντολή is.numeric και μόνο για αυτές να γίνει η αντικατάσταση, διαφορετικά για τις στήλες που περιέχουν μη αριθμητικά στοιχεία θα επιστρέψει σφάλμα. Επομένως ορίζουμε την παρακάτω συνάρτηση:

fillna=function(A)
{
  e=A
  for (i in 1:ncol(e))
  {
    if (is.numeric(e[,i]))
      {
      e[is.na(e[,i]),i]=mean(e[,i],na.rm = TRUE)
      }
  }
  return(e)
}

Ορίζουμε ένα ενδεικτικό data.frame και παίρνουμε:

x=data.frame(a=c(2,3,NA,1,4), b=c(NA,NA,2,NA,NA), c=c("A","B","C",NA,NA))
x
##    a  b    c
## 1  2 NA    A
## 2  3 NA    B
## 3 NA  2    C
## 4  1 NA <NA>
## 5  4 NA <NA>
y=fillna(x)
y
##     a b    c
## 1 2.0 2    A
## 2 3.0 2    B
## 3 2.5 2    C
## 4 1.0 2 <NA>
## 5 4.0 2 <NA>

Άσκηση 3

Δημιουργήστε μια συνάρτηση που παίρνει ως όρισμα ένα διάνυσμα αριθμών \(x=(x_1, x_2, \ldots, x_n)\) και έναν αριθμό \(A\) και επιστρέφει ένα διάνυσμα \((a,b,c)\), όπου:

(α) \(a=\sum_{j=1}^n \frac{e^{-x_j}+1}{j+2}\),

(β) \(b\) είναι ίσο με το πλήθος των στοιχείων του \(x\) που είναι μεγαλύτερα ή ίσα από το 80% του μέγιστου στοιχείου του \(x\).

(γ) \(c\) είναι το ελάχιστο απαιτούμενο πλήθος στοιχείων του \(x\) ξεκινώντας από το πρώτο, έτσι ώστε το άθροισμα τους να είναι μεγαλύτερο του \(Α\). Αν το άθροισμα όλων είναι μικρότερο ή ίσο του \(A\), τότε \(c=n+1\).

Απάντηση

Ορίζουμε τη συνάρτηση

ask3=function(x,A)
{
  n=length(x)
  e=1:n
  # Compute a
  a=sum(exp(-x)/(e+2))
  xmax=max(x)
  # Compute b
  b=sum(x>=0.8*xmax)
  # Compute c
  c=0
  s=0
  while((s<=A) & (c<=n))
  {
    c=c+1
    s=s+x[c]
  }
  #Return result
  return(c(a,b,c))
}

Για τον υπολογισμό του \(b\), η έκφραση x>=0.8xmax δίνει ένα λογικό διάνυσμα με στοιχεία TRUE/FALSE ανάλογα με το αν ισχύει η x[i]>=0.8xmax. Επειδή τα στοιχεία TRUE αντιστοιχούν σε τιμή 1 και τα FALSE σε τιμή 0, το άθροισμα στοιχείων αυτού του διανύσματος δίνει τον αριθμό στοιχείων του x που ικανοποιούν τη συνθήκη.

Εφαρμόζοντας τη συνάρτηση ask3 στο διάνυσμα

x=1:5
x
## [1] 1 2 3 4 5

παίρνουμε

ask3(x,7)
## [1] 0.1704329 2.0000000 4.0000000
ask3(x,16)
## [1] 0.1704329 2.0000000 6.0000000

Άσκηση 4

(α) Εισαγάγετε τον πίνακα x (10Χ8) που περιέχεται στο αρχείο x.Rdata.

(β) Δημιουργήστε ένα πίνακα f διαστάσεων 2Χ8, ο οποίος στην πρώτη γραμμή περιέχει τον αριθμό των στοιχείων κάθε στήλης του x που είναι μικρότερα ή ίσα από τη μέση τιμή της αντίστοιχης στήλης και στη δεύτερη γραμμή τον αριθμό των στοιχείων κάθε στήλης του x που είναι μεγαλύτερα από τη μέση τιμή της αντίστοιχης στήλης. (Για παράδειγμα, στην πρώτη στήλη του x η μέση τιμή είναι 9.6 και πέντε στοιχεία της στήλης είναι μικρότερα της μέσης τιμής ενώ 5 είναι μεγαλύτερα. Επομένως η πρώτη στήλη του πίνακα f θα πρέπει να είναι το διάνυσμα (5,5).)

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

(δ) Εφαρμόστε τη συνάρτηση abovebelowmean στον πίνακα x. Το αποτέλεσμα πρέπει να είναι το ίδιο με αυτό της ερώτησης (β).

Απάντηση

(α)

load("x.Rdata")

(β) Πρώτα υπολογίζουμε τους μέσους όρους των στηλών του x. Αυτό μπορεί να γίνει με διάφορους τρόπους. Ο πιο άμεσος είναι με τη χρήση της συνάρτησης colMeans που υπολογίζει τους μέσους όρους στηλών ενός πίνακα

xcolmeans=colMeans(x)
xcolmeans
## [1]  9.6  9.2 11.4 10.8 11.7  6.1 10.9  7.7

Εναλλακτικά μπορούμε να εφαρμόσουμε τη συνάρτηση mean στις στήλες του x μέσω της apply:

xcolmeans2=apply(x,2,mean)
xcolmeans2
## [1]  9.6  9.2 11.4 10.8 11.7  6.1 10.9  7.7

Για την έυρεση των συχνοτήτων που ζητούνται εφαρμόζουμε ένα for loop στις στήλες του x. Για τη στήλη j υπολογίζουμε τις δύο συχνότητες και αποθηκεύουμε το διάνυσμα των δύο τιμών στην j γραμμή του πίνακα f:

n=ncol(x)
f=matrix(rep(0, 2*n), n, 2)
for (j in 1:n)
{
  nbelow=sum(x[,j]<=xcolmeans[j])
  nabove=sum(x[,j]>xcolmeans[j])
  f[j,]=c(nbelow,nabove)
}
f
##      [,1] [,2]
## [1,]    5    5
## [2,]    6    4
## [3,]    3    7
## [4,]    4    6
## [5,]    5    5
## [6,]    5    5
## [7,]    3    7
## [8,]    5    5

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

Η συνάρτηση abovebelowmean μπορεί να διαρθρωθεί χρησιμοποιώντας τις εντολές μέσα στο for loop του προηγούμενου ερωτήματος:

abovebelowmean=function(v)
{
  m=mean(v)
  nbelow=sum(v<=m)
  nabove=sum(v>m)
  return(c(nbelow, nabove))
}

Για παράδειγμα εφαρμόζοντας τη συνάρτηση στην πρώτη στήλη του πίνακα παίρνουμε

abovebelowmean(x[,1])
## [1] 5 5

(δ) Τώρα μπορούμε να εφαρμόσουμε την abovebelowmean στον x μέσω της apply σε στήλες:

f1=apply(x,2,abovebelowmean)
f1
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    5    6    3    4    5    5    3    5
## [2,]    5    4    7    6    5    5    7    5

Το αποτέλεσμα f1 είναι ένας πίνακας 2χ8, όπου κάθε στήλη περιέχει το διάνυσμα συχνοτήτων της αντίστοιχης στήλης του x. Για να το φέρουμε σε μορφή συγκρίσιμη με τον f παίρνουμε τον ανάστροφο πίνακα και συγκρίνουμε τα δύο αποτελέσματα

f2=t(f1)
f2
##      [,1] [,2]
## [1,]    5    5
## [2,]    6    4
## [3,]    3    7
## [4,]    4    6
## [5,]    5    5
## [6,]    5    5
## [7,]    3    7
## [8,]    5    5
f2==f
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
## [3,] TRUE TRUE
## [4,] TRUE TRUE
## [5,] TRUE TRUE
## [6,] TRUE TRUE
## [7,] TRUE TRUE
## [8,] TRUE TRUE

Άσκηση 5

Δημιουργήστε μια συνάρτηση bincount που παίρνει ως όρισμα ένα διάνυσμα μη αρνητικών αριθμών \(x=(x_1, x_a2, \ldots, x_n)\), και ένα διάνυσμα αριθμών \(b=(b_1, b_2, \ldots, b_k)\), και επιστρέφει διάνυσμα συχνοτήτων στο οποίο η πρώτη συνιστώσα είναι ο αριθμός των στοιχείων του \(x\) που είναι μικρότερα από \(b_1\), η δεύτερη είναι ο αριθμός στοιχείων του \(x\) που είναι μεταξύ \(b_1\) και \(b_2\), κλπ μέχρι την τελευταία συνιστώσα που δίνει τον αριθμό στοιχείων του \(x\) που είναι μεγαλύτερα του \(b_k\).

Απάντηση Η συνάρτηση πρέπει να επιστρέφει ένα διάνυσμα \(v\) διάστασης \(k+1\), όπου \(v_1=\sum_{j=1}^n 1(x_j < b_1)\), \(v_2=\sum_{j=1}^n 1(b_1 \leq x_j < b_2), \ldots, v_k=\sum_{j=1}^n 1(b_{k-1}\leq x_j < b_k), v_{k+1}=\sum_{j=1}^n 1(x_j \geq b_k)\). Επειδή τα \(v_1\) και \(v_{k+1}\) έχουν διαφορετικό τύπο από τα υπόλοια, αυτά τα υπολογίζουμε χωριστά, ενώ τα \(v_2, \ldots, v_k\) τα υπολογίζουμε με ένα for loop. Ορίζουμε τη συνάρτηση

bincount = function(x,b)
{
  k=length(b)
  v=rep(0,k+1)
  v[1]=sum(x<b[1])
  for (j in 2:k)
  {
    v[j]=sum((x >= b[j-1])&(x < b[j]))
  }
  v[k+1]=sum(x >= b[k])
  return(v)
}

Εφαρμόζουμε τη συνάρτηση bincount σε ένα διάνσυμα 1000 τυχαίων αριθμών από την ομοιόμορφη κατανομή μεταξυ 0 και 1, με break points \(0.1, 0.2, \ldots, 0.9\):

x=runif(1000,0,1)
b=seq(0.1, 0.9, by=0.1)
freq=bincount(x,b)
freq
##  [1]  98  89 112 102 120 100  94  99 100  86

Γο άθροισμα των στοιχείων του freq είναι 1000 όπως αναμένεται,

sum(freq)
## [1] 1000

και επειδή το μήκος του \(x\) είναι αρκετά μεγάλο, οι σχετικές συχνότητες των 10 διαστημάτων είναι αρκετά κοντά στη θεωρητική τιμή των αντίστοιχων πιθανοτήτων που είναι όλες ίσες με 0.1.

Άσκηση 6

Δημιουργήστε ένα data.frame με 100 γραμμές και 3 στήλες και ονόματα στηλών a,b,c. Η στήλη 1 περιέχει τους αριθμούς 2,4,..200, η στήλη 2 100 τυχαίους αριθμούς από την ομοιόμορφη κατανομή μεταξύ 10 και 20 και η τρίτη στήλη τυχαίους λογικούς χαρακτήρες (TRUE/FALSE). Σώστε το data.frame σε αρχείο excel.

Απάντηση Δημιουργούμε πρώτα τις στήλες μία μία. Για την \(a\) χρησιμοποιούμε την εντολή seq:

a=seq(2,200,by=2)
a
##   [1]   2   4   6   8  10  12  14  16  18  20  22  24  26  28  30  32  34  36
##  [19]  38  40  42  44  46  48  50  52  54  56  58  60  62  64  66  68  70  72
##  [37]  74  76  78  80  82  84  86  88  90  92  94  96  98 100 102 104 106 108
##  [55] 110 112 114 116 118 120 122 124 126 128 130 132 134 136 138 140 142 144
##  [73] 146 148 150 152 154 156 158 160 162 164 166 168 170 172 174 176 178 180
##  [91] 182 184 186 188 190 192 194 196 198 200

Για την \(b\) χρησιμοποιούμε τη συνάρτηση runif(n,a,b) που επιστρέφει ένα διάνυσμα \(n\) τυχαίων αριθμών από την ομοιόμορφη κατανομή μεταξύ \(a\) και \(b\):

b=runif(100,10,20)
b
##   [1] 16.67734 19.79229 19.01414 15.99433 18.75384 16.29073 15.83449 19.50262
##   [9] 18.89757 15.21535 13.79724 16.41669 16.37337 11.07398 10.41962 11.23912
##  [17] 17.09689 11.79289 10.76668 13.71130 18.26806 15.10134 16.90762 11.35586
##  [25] 19.81611 18.33908 17.00127 14.82511 19.56485 12.79860 17.96724 10.66896
##  [33] 18.12507 16.95616 14.31442 10.44795 16.75685 17.80779 12.23390 14.01606
##  [41] 11.45490 16.03235 13.49448 16.19527 19.45271 15.25154 12.25195 15.40252
##  [49] 15.13080 15.26687 15.86154 19.99663 10.00757 14.88905 15.83804 11.69791
##  [57] 17.65338 14.93315 16.95797 19.79120 19.73575 15.20441 11.41364 19.65142
##  [65] 18.34496 14.39214 12.15193 15.81267 13.09884 14.37358 17.59721 10.55633
##  [73] 15.98517 10.92114 10.07945 19.63482 11.28546 10.49644 18.31529 13.91114
##  [81] 17.75401 18.25147 13.79867 18.38516 13.41122 18.43285 18.62420 11.49516
##  [89] 10.80017 19.86131 14.44442 14.26737 12.20773 16.63682 11.65793 17.61107
##  [97] 12.63906 11.63414 13.53653 11.87014

Για τη στήλη \(c\) πρέπει να δημιουργήσουμε ένα διάνυσμα λογικών μεταβλητών κάθε μια από τις οποίες είναι TRUE/FALSE με πιθανότητες 1/2, 1/2. Ένας τρόπος για να γίνει αυτό είναι να δημιουργήσουμε ένα διάνυσμα από 0 ή 1 ισοπίθανα, και να το μετατρέψουμε σε logical με την as.logical. Για να δημιουργήσουμε έναν τυχαίο αριθμό 0 ή 1 με πιθανότητες 1/2, 1/2 παρατηρούμε ότι αντιστοιχεί σε κατανομή Bernoulli(1/2), δηλαδή Binomial(1, 1/2). Επομένως μπορούμε να χρησιμοποιήσουμε τη συνάρτηση rbinom(k,n,p) που επιστρέφει \(n\) τυχαίους αριθμούς από την κατανομή Binomial(n,p):

c=as.logical(rbinom(100,1,1/2))
c
##   [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
##  [13]  TRUE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE
##  [25]  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE
##  [37]  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
##  [49]  TRUE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE  TRUE
##  [61] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE
##  [73] FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE
##  [85]  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
##  [97] FALSE  TRUE  TRUE  TRUE

Τώρα μπορούμε να ορίσουμε το ζητούμενο dataframe:

A=data.frame(a,b,c)
A[1:10,]
##     a        b     c
## 1   2 16.67734  TRUE
## 2   4 19.79229  TRUE
## 3   6 19.01414  TRUE
## 4   8 15.99433  TRUE
## 5  10 18.75384  TRUE
## 6  12 16.29073  TRUE
## 7  14 15.83449  TRUE
## 8  16 19.50262 FALSE
## 9  18 18.89757  TRUE
## 10 20 15.21535  TRUE

Για να το σώσουμε, χρησιμοποιούμε την εντολή write.table που σώζει το dataframe σε μορφή csv, η οποία μπορεί να εισαχθεί στο excel.

write.table(A,file="A.csv", sep=';', col.names = T)

Άσκηση 7

Το αρχείο water.xlsx περιέχει δεδομένα σχετικά με την κατανάλωση νερού ενός αριθμού ατόμων. Για κάθε άτομο δίνονται η κατανάλωση νερού σε ένα έτος (σε κ.μ.), το εισόδημα του ατόμου (σε ευρω/μήνα) και το εργασιακό καθεστώς (εργαζόμενος/συνταξιούχος).

(α) Εισαγάγετε το αρχείο σε ένα data.frame

(β) Υπολογίστε τη μέση κατανάλωση νερού των εργαζόμενων και των συνταξιούχων.

(γ) Υπολογίστε το συντελεστή συσχέτισης ανάμεσα στο εισόδημα και την κατανάλωση νερού για τους εργαζόμενους, τους συνταξιούχους και όλα τα άτομα του δείγματος.

(δ) Επαναλάβετε το (γ) έχοντας αποκλείσει τα άτομα με εισόδημα μικρότερο από 500 ευρώ το μήνα ή κατανάλωση νερού πάνω από 400 κμ.

(ε) Δημιουργήστε ένα γράφημα της κατανάλωσης νερού ως προς το εισόδημα.

(στ) Επαναλάβετε το (ε) όπου στο γράφημα διακρίνονται με κάποιο τρόπο οι εργαζόμενοι από τους συνταξιούχους (διαφορετικό χρώμα ή σύμβολο ή και τα δύο)

Απάντηση

(α) Για να εισαγάγουμε το αρχείο, πρώτα το σώζουμε σε μορφή csv με διαχωριστικό χαρακτήρα ανάμεσα στα πεδία το “;”.

w=read.csv("water.csv",sep = ';')

(β)

waterconsworking=mean(w[w$retire=='no',1])
waterconsworking
## [1] 256.0286
waterconsretired=mean(w[w$retire=='yes',1])
waterconsretired
## [1] 167.0548

(γ)

corworking=cor(w[w$retire=='no',1],w[w$retire=='no',2])
corworking
## [1] 0.3533838
corretired=cor(w[w$retire=='yes',1],w[w$retire=='yes',2])
corretired
## [1] 0.3518266
corall=cor(w[,1],w[,2])
corall
## [1] 0.4177872

Άσκηση 8

Το αρχείο heart.csv περιέχει δεδομένα σχετικά με ασθενείς που ήρθαν στα επείγοντα ενός νοσοκομείου με πόνο στο στήθος. Η στήλη AHD έχει το τελικό αποτέλεσμα των καρδιολογικών εξετάσεων σχετικά με το αν ο ασθενής διαγνώστηκε με καρδιακή νόσο ή όχι. Οι υπόλοιπες στήλες περιέχουν την ηλικία, το φύλο, την κατηγοριοποίηση του πόνου στήθους (ChestPain), την αρτηριακή πίεση σε ηρεμία (RestBP), και την χοληστερόλη (Chol).

(α) Εισαγάγετε το αρχείο σε ένα data.frame με το όνομα heart.

(β) Προσθέστε μια νέα στήλη στο heart που περιέχει τη χοληστερόλη κατηγοριοποιημένη σε 4 κατηγορίες (low, medium, high, veryhigh) ανάλογα με το αν η μέτρηση χοληστερόλης είναι μικρότερη από το 1ο τεταρτημόριο, μεταξύ 1ου και 2ου, μεταξύ 2ου και 3ου, και μεγαλύτερη από το 3ο τεταρτημόριο, αντίστοιχα.

Η απάντηση να δοθεί με δύο τρόπους: (1) με κατάλληλο for loop για όλες τις γραμμές του dataframe, (2) χρησιμοποιώντας συνδυασμό της συνάρτησης cut με τη συνάρτηση quantile (δείτε το help).

(γ) Δημιουργήστε δύο πίνακες συνάφειας με τις συχνότητες της καρδιακής νόσου ως προς την κατηγορία του πόνου και ως προς την κατηγορία χοληστερόλης (η εντολή table δημιουργεί πίνακες συνάφειας).

(δ) Δημιουργήστε ένα σχήμα που αποτελείται από 3 μικρότερα σχήματα σε μια σειρά. Το πρώτο σχήμα περιλαμβάνει δύο χωριστά boxplot της αρτηριακής πίεσης για άτομα χωρίς και με καρδιακή νόσο. Το δεύτερο σχήμα περιλαμβάνει δύο χωριστά scatterplot της χοληστερόλης ως προς την ηλικία, για άτομα χωρίς και με καρδιακή νόσο, με χρώματα μπλε και κόκκινο, αντίστοιχα. Το τρίτο σχήμα περιλαμβάνει ένα διάγραμμα όλων των μετρήσεων χοληστερόλης με τρεις οριζόντιες γραμμές στα επίπεδα της μέσης τιμής και δύο τυπικών αποκλίσεων πάνω και κάτω από τη μέση τιμή. (Σημείωση: Η εντολή plot(x), χωρίς δεύτερο διάνυσμα δημιουργεί διάγγραμα με τις τιμές του x στον κάθετο άξονα και τον αριθμό παρατήρησης στον οριζόντιο άξονα.)

(ε) Πραγματοποιήστε ένα t-test για το αν υπάρχει στατιστικά σημαντική διαφορά στη μέση αρτηριακή πίεση μεταξύ ατόμων με και χωρίς καρδιακή νόσο. Περάστε την τιμή του p-value στη μεταβλητή p. 

Άσκηση 9

Το αρχείο graduates.Rdata περιέχει ένα data.frame με ένα τυχαίο δείγμα από φοιτητές/απόφοιτους ενός ακαδημαϊκού τμήματος. Οι στήλες είναι ο Αριθμός Μητρώου, το φύλο, ο χρόνος πρώτης εγγραφής, ο χρόνος αποφοίτησης αν έχει αποφοιτήσει ή ο χρόνος διαγραφής αν έχει διαγραφεί, η ακαδημαϊκή κατάσταση (C=ενεργός φοιτητής, D=έχει διαγραφεί, G=έχει πάρει πτυχίο) και ο βαθμός πτυχίου αν έχει αποφοιτήσει.

Το αρχείο μπορεί να εισαχθεί με την εντολή load(“graduates.Rdata”).

(α) Προσθέστε μια μεταβλητή στο data.frame, με όνομα Duration που περιέχει τη διάρκεια σπουδών για όσους έχουν αποφοιτήσει (για τους διαγραφέντες και όσους συνεχίζουν θα είναι NA).

(β) Δημιουργήστε ένα πίνακα συχνοτήτων που περιέχει στοιχεία για τη διάρκεια λήψης πτυχίου ανά έτος εισαγωγής ως εξής: κάθε γραμμή του πίνακα περιέχει απόφοιτους που έχουν εισαχθεί κατά το έτος 2010, 2011, κλπ. Κάθε στήλη περιέχει τον αριθμό που μπήκαν στο τμήμα αυτή τη χρονιά και πήραν πτυχίο σε 4, 5, 6 ,… χρόνια. Επίσης ο πίνακας περιέχει τα αθροίσματα γραμμών στην τελευταία στήλη και τα αθροίσματα στηλών στην τελευταία γραμμή:

YearIn 4 5 6 7 Total
2010 25 30 35 15 105
2010 30 20 15 20 85
————- ————- —– —– —– —–
Total 55 50 50 35 190

(γ) Επαναλάβετε το ερώτημα (β) ανά έτος αποφοίτησης.

(δ) Προσθέστε μια μεταβλητή στο data.frame που περιέχει τον χαρακτηρισμό του πτυχίου σε Άριστα (>=8.5), Λίαν Καλώς (7-8.49) και Καλώς (5-6.99). Δημιουργήστε ένα πίνακα αντίστοιχο του ερωτήματος (γ) που περιέχει συχνότητες αποφοίτων με τον χαρακτηρισμό του πτυχίου, ανά έτος αποφοίτησης (μαζί με τα σύνολα γραμμών και στηλών).

(ε) Δημιουργήστε ένα scatter plot του μέσου βαθμού πτυχίου των αποφοίτων ανά έτος αποφοίτησης στο οποίο διακρίνονται οι απόφοιτοι ανάλογα με το φύλο.

(στ) Ελέγξτε σε επίπεδο σημαντικότητας 5% αν υπάρχει στατιστικά σημαντική διαφορά στο μέσο βαθμό πτυχίου ανάμεσα σε άντρες και γυναίκες αποφοίτους.

Απάντηση

Το αρχείο εισάγεται με την εντολή

load("graduates.Rdata")
dim(graduates)
## [1] 800   6
graduates[1:10,]
##      AM Sex YearIn YearOut Status  GPA
## 1  1001   F   2016    2020      G 8.49
## 2  1002   F   2011    2015      G 6.27
## 3  1003   M   2015    2019      G 6.44
## 4  1004   M   2014    2015      D   NA
## 5  1005   F   2011    2018      G 7.14
## 6  1006   M   2011    2016      G 6.12
## 7  1007   M   2019      NA      C   NA
## 8  1008   M   2017      NA      C   NA
## 9  1009   M   2015    2021      G 5.80
## 10 1010   F   2017      NA      C   NA

και δημιουργείται το dataframe graduates. Το ότι το αρχείο και το dataframe έχουν το ίδιο όνομα δεν είναι αναγκαίο, απλώς ορίστηκαν έτσι. Επίσης ένα αρχείο .Rdata μπορεί να περιέχει περισσότερα από ένα αντικείμενα, τύπου όχι μόνο dataframes, αλλά γενικά.

(α)

graduates$Duration=graduates$YearOut-graduates$YearIn
graduates[1:10,]
##      AM Sex YearIn YearOut Status  GPA Duration
## 1  1001   F   2016    2020      G 8.49        4
## 2  1002   F   2011    2015      G 6.27        4
## 3  1003   M   2015    2019      G 6.44        4
## 4  1004   M   2014    2015      D   NA        1
## 5  1005   F   2011    2018      G 7.14        7
## 6  1006   M   2011    2016      G 6.12        5
## 7  1007   M   2019      NA      C   NA       NA
## 8  1008   M   2017      NA      C   NA       NA
## 9  1009   M   2015    2021      G 5.80        6
## 10 1010   F   2017      NA      C   NA       NA

(β) Για την πινακοποίηση του Duration μπορούμε να χρησιμοποιήσουμε την εντολή table που δημιουργεί 2Χ2 πίνακες συχνοτήτων:

d=table(graduates$YearIn,graduates$Duration)
d
##       
##         0  1  2  3  4  5  6  7  8  9 10 11
##   2010  1  0  0  0  6 11 11 18 13  5  4  9
##   2011  1  2  0  0 18 10  4 17 17  8 11  0
##   2012  1  2  0  0 15  7 12 16 13 12  0  0
##   2013  1  0  2  1 16 14  9  8 15  0  0  0
##   2014  1  2  0  3 16 30 21 15  0  0  0  0
##   2015  0  1  0  1 29 26 37  0  0  0  0  0
##   2016  0  1  1  3 33 28  0  0  0  0  0  0
##   2017  1  0  1  3  1  0  0  0  0  0  0  0
##   2018  1  2  0  2  0  0  0  0  0  0  0  0
##   2019  3  1  1  0  0  0  0  0  0  0  0  0
sum(d)
## [1] 574

Παρατηρούμε ότι ο συνολικός αριθμός των φοιτητών που έχουν καταμετρηθεί στον πίνακα d είναι 574, επειδή στη συνάρτηση table δεν χρησιμοποιήθηκαν οι γραμμές με Duration=NA. Όμως υπάρχουν γραμμές με Duration=0, 1, 2 ή 3, που αντιστοιχούν σε φοιτητές που έχουν διαγραφεί και δεν πρέπει να καταμετρηθούν στον πίνακα. Ένας τρόπος να το αντιμετωπίσουμε είναι να δημιουργήσουμε ένα νέο dataframe μόνο με τους πτυχιούχους:

ptyx=graduates[graduates$Status=="G",]
dim(ptyx)
## [1] 503   7

Η συνάρτηση table στο νέο dataframe δίνει

duryearin=table(ptyx$YearIn, ptyx$Duration)
duryearin
##       
##         4  5  6  7  8  9 10 11
##   2010  6  7 10 16 13  5  4  9
##   2011 14 10  4 17 16  8 10  0
##   2012 14  6 12 15 13 11  0  0
##   2013 14 13  9  8 15  0  0  0
##   2014 15 29 20 14  0  0  0  0
##   2015 27 25 35  0  0  0  0  0
##   2016 33 26  0  0  0  0  0  0
sum(duryearin)
## [1] 503

και βλέπουμε ότι όλοι οι φοιτητές έχουν χρόνο λήψης πτυχίου τουλάχιστον 4. Το αντικείμενο d είναι τύπου πίνακα:

is.matrix(d)
## [1] TRUE

και η πρώτη στήλη με τα έτη εισαγωγής που εμφανίζεται στην οθόνη δεν είναι στήλη του πίνακα αλλά τα ονόματα των γραμμών όπως και η πρώτη γραμμή τα ονόματα των στηλών:

dimnames(duryearin)
## [[1]]
## [1] "2010" "2011" "2012" "2013" "2014" "2015" "2016"
## 
## [[2]]
## [1] "4"  "5"  "6"  "7"  "8"  "9"  "10" "11"

Για να προσθέσουμε τη γραμμή και τη στήλη με τα αθροίσματα στηλών και γραμμών αντίστοιχα μπορούμε να χρησιμοποιήσουμε την εντολή apply:

rowsums=apply(duryearin,1,sum)
rowsums
## 2010 2011 2012 2013 2014 2015 2016 
##   70   79   71   59   78   87   59

Μπορούμε να προσθέσουμε το διάνυσμα rowsums στην τελευταία στήλη του duryearin με την εντολή cbind:

duryearin=cbind(duryearin, rowsums)
duryearin
##       4  5  6  7  8  9 10 11 rowsums
## 2010  6  7 10 16 13  5  4  9      70
## 2011 14 10  4 17 16  8 10  0      79
## 2012 14  6 12 15 13 11  0  0      71
## 2013 14 13  9  8 15  0  0  0      59
## 2014 15 29 20 14  0  0  0  0      78
## 2015 27 25 35  0  0  0  0  0      87
## 2016 33 26  0  0  0  0  0  0      59

Αντίστοιχα για τα αθροίσματα των στηλών

colsums=apply(duryearin,2,sum)
duryearin=rbind(duryearin,colsums)
duryearin
##           4   5  6  7  8  9 10 11 rowsums
## 2010      6   7 10 16 13  5  4  9      70
## 2011     14  10  4 17 16  8 10  0      79
## 2012     14   6 12 15 13 11  0  0      71
## 2013     14  13  9  8 15  0  0  0      59
## 2014     15  29 20 14  0  0  0  0      78
## 2015     27  25 35  0  0  0  0  0      87
## 2016     33  26  0  0  0  0  0  0      59
## colsums 123 116 90 70 57 24 14  9     503

(γ)

duryearout=table(ptyx$YearOut, ptyx$Duration)
duryearout
##       
##         4  5  6  7  8  9 10 11
##   2014  6  0  0  0  0  0  0  0
##   2015 14  7  0  0  0  0  0  0
##   2016 14 10 10  0  0  0  0  0
##   2017 14  6  4 16  0  0  0  0
##   2018 15 13 12 17 13  0  0  0
##   2019 27 29  9 15 16  5  0  0
##   2020 33 25 20  8 13  8  4  0
##   2021  0 26 35 14 15 11 10  9
sum(duryearout)
## [1] 503
rowsums=apply(duryearout,1,sum)
duryearout=cbind(duryearout, rowsums)
colsums=apply(duryearout,2,sum)
duryearout=rbind(duryearout,colsums)
duryearout
##           4   5  6  7  8  9 10 11 rowsums
## 2014      6   0  0  0  0  0  0  0       6
## 2015     14   7  0  0  0  0  0  0      21
## 2016     14  10 10  0  0  0  0  0      34
## 2017     14   6  4 16  0  0  0  0      40
## 2018     15  13 12 17 13  0  0  0      70
## 2019     27  29  9 15 16  5  0  0     101
## 2020     33  25 20  8 13  8  4  0     111
## 2021      0  26 35 14 15 11 10  9     120
## colsums 123 116 90 70 57 24 14  9     503

(δ) Δημιουργούμε μια συνάρτηση ptyxcat που διακρίνει τις περιπτώσεις βαθμού πτυχίου

ptyxcat=function(x)
{
  if (x>=8.5) p="Arista" else if (x>=7) p="LianKalws" else p="Kalws"
  return(p)
}

Χρησιμποποιούμε την εντολή sapply για να εφαρμόσουμε τη συνάρτηση ptyxcat στο διάνυσμα των βαθμών πτυχίου του dataframe ptyx και προσθέτουμε το διάνυσμα των κατηγοριών στο dataframe με την εντολή cbind:

pcat=sapply(ptyx$GPA, ptyxcat)
ptyx=cbind(ptyx, pcat)
ptyx[1:10,]
##      AM Sex YearIn YearOut Status  GPA Duration      pcat
## 1  1001   F   2016    2020      G 8.49        4 LianKalws
## 2  1002   F   2011    2015      G 6.27        4     Kalws
## 3  1003   M   2015    2019      G 6.44        4     Kalws
## 5  1005   F   2011    2018      G 7.14        7 LianKalws
## 6  1006   M   2011    2016      G 6.12        5     Kalws
## 9  1009   M   2015    2021      G 5.80        6     Kalws
## 11 1011   F   2011    2018      G 6.50        7     Kalws
## 14 1014   M   2011    2018      G 7.47        7 LianKalws
## 16 1016   F   2016    2020      G 5.90        4     Kalws
## 17 1017   F   2014    2020      G 7.18        6 LianKalws

Τώρα πινακοποιούμε την κατηγορία πτυχίου ανά έτος αποφοίτησης, ανάλογα με το ερώτημα (γ):

catyearout=table(ptyx$YearOut, ptyx$pcat)
rowsums=apply(catyearout,1,sum)
catyearout=cbind(catyearout, rowsums)
colsums=apply(catyearout,2,sum)
catyearout=rbind(catyearout,colsums)
catyearout
##         Arista Kalws LianKalws rowsums
## 2014         0     4         2       6
## 2015         3    12         6      21
## 2016         5    16        13      34
## 2017         9    20        11      40
## 2018        10    40        20      70
## 2019        13    52        36     101
## 2020        18    52        41     111
## 2021        20    63        37     120
## colsums     78   259       166     503

Επειδή η πινακοποίηση έγινε αλφαβητικά ως προς το pcat, για να εμφανιστούν οι κατηγορίες με τη σειρά Καλώς, Λίαν Καλώς και Άριστα, μπορούμε να κάνουμε αναδιάταξη των στηλών του catyearout:

catyearout=catyearout[,c(2,3,1,4)]
catyearout
##         Kalws LianKalws Arista rowsums
## 2014        4         2      0       6
## 2015       12         6      3      21
## 2016       16        13      5      34
## 2017       20        11      9      40
## 2018       40        20     10      70
## 2019       52        36     13     101
## 2020       52        41     18     111
## 2021       63        37     20     120
## colsums   259       166     78     503

(ε) Θα διακρίνουμε το φύλο με διαφορετικό χρώμα, κόκκινο για F και μπλε για M. Η στήλη Sex είναι τύπου factor:

is.factor(ptyx$Sex)
## [1] TRUE
ptyx[1:5,2]
## [1] F F M F M
## Levels: F M
as.numeric(ptyx[1:5,2])
## [1] 1 1 2 1 2

Τα επίπεδα F,M αντιστοιχούν σε 1 και 2 αντίστοιχα. Τώρα δημιουργούμε ένα διάνυσμα δύο χρωμάτων που θέλουμε να αντιστοιχούν στα επίπεδα F, M:

pcolors=c("red","blue")

Τα ονόματα των χρωμάτων είναι από τη λίστα των χρωμάτων που αναγνωρίζει η εντολή plot. Τώρα μπορούμε να καλέσουμε την εντολή plot για να δημιουργήσουμε ένα scatter plot μεταξύ GPA και YearOut, χρησιμοποιώντας για χρώμα την επιλογή col=pcolors[ptyx$Sex]. Αυτό σημαίνει ότι για κάθε σημείο το χρώμα θα είναι αυτό που αντιστοιχεί στο pcolors[1] αν πρόκειται για επίπεδο 1 (F) και στο pcolors[2] αν πρόκειται για επίπεδο 2 (M). H επιλογή pch ορίζει το σύμβολο κάθε σημείου και η επιλογή cex το μέγεθος σε κλάσμα του default:

plot(ptyx$GPA~ptyx$YearOut, col=pcolors[ptyx$Sex],cex=0.7,pch="+")

Με αντίστοιχο τρόπο θα μπορούσαμε να διαφοροποιήσουμε τα σημεία του διαγράμματος και ως προς το σύμβολο:

psymbols=c(16,17)
plot(ptyx$GPA~ptyx$YearOut, col=pcolors[ptyx$Sex],cex=0.7,pch=psymbols[ptyx$Sex])

Για τους κωδικούς των συμβόλων δείτε το help(points).

(στ) Πρόκειται για έλεγχο διαφοράς μέσων τιμών δύο πληθυσμών με άγνωστη διασπορά. Εφαρμόζουμε τον δίπλευρο έλεγχο \(t\) με υποθέσεις \[H_0 : \mu_F = \mu_M, H_1 : \mu_F \neq \mu_M\]:

gpatest=t.test(ptyx[ptyx$Sex=="F",6], ptyx[ptyx$Sex=="M",6],conf.level = 0.95, alternative = "two.sided")
gpatest
## 
##  Welch Two Sample t-test
## 
## data:  ptyx[ptyx$Sex == "F", 6] and ptyx[ptyx$Sex == "M", 6]
## t = -1.7945, df = 497.95, p-value = 0.07334
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4464595  0.0202206
## sample estimates:
## mean of x mean of y 
##  6.908945  7.122065

Από το p-value=0.0733429 προκύπτει ότι δεν υπάρχει στατιστικά σημαντική διαφορά στο μέσο βαθμό πτυχίου μεταξύ των δύο φύλων.

Άσκηση 10

Σκοπός της άσκησης είναι να ελέγξει μέσω προσομοίωσης ότι ο τύπος για το διάστημα εμπιστοσύνης του μέσου μιας κανονικής κατανομής με γνωστή διασπορά \[ \overline{x}_n \pm z_{1-\alpha/2} \frac{\sigma}{\sqrt{n}} \] όντως έχει πιθανότητα κάλυψης ίση με το επίπεδο εμπιστοσύνης \(1-\alpha\), όπου \(\overline{x}_n\) είναι ο δειγματικός μέσος ενός τυχαίου δείγματος μεγέθους \(n\) από την κανονική κατανομή με άγνωστη μέση τιμή \(\mu\) και γνωστή διασπορά \(\sigma^2\).

Δημιουργήστε μια συνάρτηση confintverify που παίρνει ως ορίσματα \(\mu, \sigma, \alpha, n, N\) και

  1. Δημιουργεί ένα πίνακα \(Ν\times n\) από \(N\) ανεξάρτητα μεταξύ τους τυχαία δείγματα μεγέθους \(n\) από την κανονική κατανομή \({\cal N}(\mu, \sigma^2)\) (αρκεί να δημιουργηθεί ένας πίνακας \(N\times n\) από αυτή την κανονική κατανομή).

  2. Για κάθε γραμμή του πίνακα (δηλαδή για καθένα από τα τυχαία δείγματα) υπολογίζει τα άκρα του \((1-\alpha)\cdot 100\%\) διαστήματος εμπιστοσύνης για το \(\mu\) χρησιμποιώντας τον παραπάνω τύπο.

  3. Δημιουργεί ένα dataframe με όνομα results 3 στηλών, του οποίου οι δύο πρώτες στήλες περιέχουν τα άκρα των \(N\) διαστημάτων εμπιστοσύνης και η τρίτη στήλη τις λογικές μεταβλητές που δείχνουν αν κάθε ένα από τα \(N\) διαστήματα εμπιστοσύνης που έχουν υπολογιστεί περιέχει την πραγματική μέση τιμή \(\mu\).

  4. Υπολογίζει το ποσοστό \(p\) των διαστημάτων εμπιστοσύνης που καλύπτουν την πραγματική μέση τιμή \(\mu\).

Η συνάρτηση θα επιστρέφει ως αποτέλεσμα το dataframe results και το ποσοστό κάλυψης \(p\).

Απάντηση

Ορίζουμε τη συνάρτηση:

confintverify = function(m, s, a, n, N)
{
  # Create N samples of size n each from N(m,s)
  samplematrix=matrix(rnorm(N*n, m, s), N, n)
  # Initialize the dataframe of results
  cilower=rep(0,N)
  ciupper=rep(0,N)
  covermu=rep(F,N)
  results=data.frame(cilower, ciupper,covermu)
  # Compute z_(1-a/2) and the half length of the confidence interval (which is constant and independent of the sample
  # since sigma is known)
  z=qnorm(1-a/2, 0, 1)
  halflength = z *s/sqrt(n)
  # Create a for loop to compute the sample mean for each row of the samplematrix and the corresponding C.I. 
  for (i in 1:N)
  {
    xbar=mean(samplematrix[i,])
    cl=xbar-halflength
    cu=xbar+halflength
    cv=(cl <= m)&(m <= cu)
    results[i, ]=c(cl, cu, cv)
  }
  # Compute the number of correct intervals and their proportion. Note that the sum of a vector of logical variables is the number of them that are true (why?)
  p=sum(results$covermu)/N
  # Create a list to return results and p
  y=list(results, p)
  names(y)=c("results", "p")
  return(y)
}

Καλούμε τη συνάρτηση με 1000 δείγματα μεγέθους 30 το καθένα, και διαστήματα εμπιστοσύνης 99%:

cint=confintverify(10, 4, 0.01, 30, 1000)
cint$results[1:10,]
##     cilower  ciupper covermu
## 1  7.262364 11.02460       1
## 2  7.515943 11.27818       1
## 3  6.704079 10.46632       1
## 4  8.597019 12.35926       1
## 5  7.287609 11.04985       1
## 6  7.805713 11.56795       1
## 7  9.399272 13.16151       1
## 8  8.184339 11.94658       1
## 9  7.122869 10.88511       1
## 10 9.286242 13.04848       1
cint$p
## [1] 0.987

Παρατηρούμε ότι το ποσοστό κάλυψης είναι πολύ κοντά στο 0.99.