Cheatsheet

Случайни експерименти

Избор на случайни елементи от вектор с връщане:

    sample(x=1:5, size=3, replace=T)

[1] 4 2 3

Избор на случайни елементи от вектор без връщане:

    sample(x=1:5, size=3, replace=F)

[1] 1 2 5

Генериране на пермутация:

    x <- 1:5
    sample(x, length(x), replace=F)

[1] 4 5 1 3 2

Намиране на честота на събитие (приблизителна вероятност):

В отдел на фирма работят 20 човека. За Коледа те решават да си разменят подаръци. В кутия слагат 20 листчета, на всяко от които има едно име. Всеки тегли листче (без да го връща) и подарява на този, чието име е изтеглил. Каква е вероятността поне един да изтегли своето име?

sim <- function(num_people) {
    x <- sample( 1:num_people, num_people, replace=F ) # permute
    d <- x - 1:num_people # diff between personal and received number
    any(d==0) # check if for any there is no diff between the two numbers
}

prob <- function(times, num_people) {
    result <- replicate(times, sim(num_people))
    sum(result)/length(result)
}

prob(100000, 20)

[1] 0.63365

Условна вероятност:

Имаме 3 карти: първата е бяла от двете страни, втората е черна от двете страни, а третата е бяла от едната и черна от другата страна. Всяка карта е поставена в затворена кутия. Избираме произволна кутия, отваряме я и виждаме, че горната страна на картата в нея е бяла. Каква е вероятността другата страна на картата също да е бяла?

sim <- function() {
    cards <- sample(c("bb", "ww", "bw"), 1)
    side <- sample(1:2, 1 )
    up <- substr(cards, start=side, stop=side)
    c(up, cards)
}

prob <- function(times) {
    result <- replicate(times, sim())
    sum(result[2,]=="ww") / sum(result[1,]=="w")
}

prob(100000)

[1] 0.665913

Случайни величини

Дискретни случайни величини

Бернулиево разпределение

X - брой успехи

successes <- 10
tries <- 20
success_probability <- 0.5
number_observations <- 10

# P(X = successes)
dbinom(successes, tries, success_probability)

[1] 0.1761971

# P(X <= successes)
pbinom(successes, tries, success_probability)

[1] 0.5880985

# generate number_observations numbers from the Binomial distribution
rbinom(number_observations, tries, success_probability)

 [1] 12 12  9  8 13 10 12 11 10 10

Геометрично разпределение

X - брой опити до успех

tries <- 2
success_probability <- 0.5
number_observations <- 10

# P(X = tries)
dgeom(tries - 1, success_probability)

[1] 0.25

# P(X <= tries)
pgeom(tries - 1, success_probability)

[1] 0.75

# generate number_observations numbers from the Geometric distribution
rgeom(number_observations, success_probability) + 1

 [1] 2 1 3 1 1 1 2 4 6 1

X - брой неуспехи преди първи успех

tries <- 2
success_probability <- 0.5
number_observations <- 10

# P(X = tries)
dgeom(tries, success_probability)

[1] 0.125

# P(X <= tries)
pgeom(tries, success_probability)

[1] 0.875

# generate number_observations numbers from the Geometric distribution
rgeom(number_observations, success_probability)

 [1] 1 2 0 2 1 1 0 0 1 0

Отрицателно биномно разпределение

X - брой опити преди r-ти успех включително

successes <- 10
tries <- 20
success_probability <- 0.5
number_observations <- 10

# P(X = tries)
dnbinom(tries - successes, successes, success_probability)

[1] 0.08809853

# P(X <= tries)
pnbinom(tries - successes, successes, success_probability)

[1] 0.5880985

# generate number_observations numbers from the Negative Binomial distribution
rnbinom(number_observations, successes, success_probability) + successes

 [1] 19 14 17 17 20 18 17 15 14 18

X - брой неуспехи преди r-ти успех

successes <- 10
failures <- 10
success_probability <- 0.5
number_observations <- 10

# P(X = failures)
dnbinom(failures, successes, success_probability)

[1] 0.08809853

# P(X <= tries)
pnbinom(failures, successes, success_probability)

[1] 0.5880985

# generate number_observations numbers from the Negative Binomial distribution
rnbinom(number_observations, successes, success_probability)

 [1] 13 15 21  5  7 12  4 21 12  8

Хипергеометрично резпределение

N - общ брой топки

M - брой бели топки

X - брой извадени бели топки (при извадени L без връщане)

chosen_white <- 5
white <- 10
all <- 20
tries <- 10
number_observations <- 10

# P(X = chosen_white)
dhyper(chosen_white, white, all - white, tries)

[1] 0.3437182

# P(X <= chosen_white)
phyper(chosen_white, white, all - white, tries)

[1] 0.6718591

# generate number_observations numbers from the Hypergeometric distribution
rhyper(number_observations, white, all - white, tries)

 [1] 8 6 5 6 6 5 7 5 7 5

Поасоново разпределение

X - брой случвания на събитието

occurences <- 1
rate <- 0.6
number_observations <- 10

# P(X = occurences)
dpois(occurences, rate)

[1] 0.329287

# P(X <= occurences)
ppois(occurences, rate)

[1] 0.8780986

# generate number_observations numbers from the Poisson distribution
rpois(number_observations, rate)

 [1] 2 1 2 0 2 0 0 1 1 1

Непрекъснати случайни величини

Равномерно разпределение

min <- 0
max <- 10
x <- 4
p <- 0.5
number_observations <- 5

# f(x)
dunif(x, min, max)

[1] 0.1

# P(X <= x)
punif(x, min, max)

[1] 0.4

# Q(x)
qunif(p, min, max)

[1] 5

# generate number_observations numbers from the Uniform distribution
runif(number_observations, min, max)

[1] 6.234611 9.345693 5.801829 7.479652 0.631416

Експоненциално разпределение

В теорията на вероятностите и статистиката, експоненциалното разпределение е разпределение на времето между събитията в процес на точка на Поасон, т.е. процес, при който събитията се случват непрекъснато и независимо с постоянна средна скорост.

rate <- 0.5
x <- 4
p <- 0.5
number_observations <- 5

# f(x)
dexp(x, rate)

[1] 0.06766764

# P(X <= x)
pexp(x, rate)

[1] 0.8646647

# Q(x)
qexp(p, rate)

[1] 1.386294

# generate number_observations numbers from the Exponential distribution
rexp(number_observations, rate)

[1] 1.464470 1.279188 2.916724 1.401556 3.430103

Нормално разпределение

mean <- 5
standard_deviation <- 2
x <- 6
p <- 0.5
number_observations <- 5

# f(x)
dnorm(x, mean, standard_deviation)

[1] 0.1760327

# P(X <= x)
pnorm(x, mean, standard_deviation)

[1] 0.6914625

# Q(x)
qnorm(p, mean, standard_deviation)

[1] 5

# generate number_observations numbers from the Normal distribution
rnorm(number_observations, mean, standard_deviation)

[1] 9.374021 7.834550 5.923814 9.338697 7.004951

Данни. Таблици и графики

Данни от тест

exam <- read.table("../datasets/examAB.txt", header=T)
exam

variant <- exam$variant

class(variant)

[1] "character"

character

table(variant)

variant
 A  B 
33 35 

table(variant) / length(variant)

variant
        A         B 
0.4852941 0.5147059 

sort(table(variant) / length(variant), decreasing=T)

variant
        B         A 
0.5147059 0.4852941 

barplot(table(variant))

barplot(table(variant) / length(variant))

barplot(sort(table(variant) / length(variant), decreasing=T))

pie(table(variant))

Други числови данни

wait <- c(2,3,3,5,5,2,7,10,4,3,1,7,11,10,5,6,3,8,5,12,5,3,8,5,7)

table(wait)

wait
 1  2  3  4  5  6  7  8 10 11 12 
 1  2  5  1  6  1  3  2  2  1  1 

wait_grouped <- cut(wait, breaks=seq(0,12,2))
table(wait_grouped)

wait_grouped
  (0,2]   (2,4]   (4,6]   (6,8]  (8,10] (10,12] 
      3       6       7       5       2       2 

table(wait_grouped)/length(wait)

wait_grouped
  (0,2]   (2,4]   (4,6]   (6,8]  (8,10] (10,12] 
   0.12    0.24    0.28    0.20    0.08    0.08 

barplot(table(wait_grouped))

hist(wait)

hist(wait, probability=T)

hist(wait, breaks=seq(0, 12, 3))

# хистограма с честотен полигон
h <- hist(wait, probability=T)
lines(
    x=c(
        min(h$breaks),
        h$mids,
        max(h$breaks)
    ),
    y=c(
      0,
      h$density,
      0
    ),
    type="l",
    lwd=2,
    col="darkorange1"
)

stripchart(
    wait,
    method="stack",
    pch=20,
    cex=1.5
)

# диаграма клон с листа (stem-and-leaf plot)
# интервали: [0, 4], (4, 9], (9, 12]
stem(wait)

  The decimal point is 1 digit(s) to the right of the |

  0 | 122333334
  0 | 555555677788
  1 | 0012

# увеличаване на броя подинтервали
# интервали: [0, 1] (1, 3] (3, 5], (5, 7], (7, 9], (9, 11] (11, 13]
stem(wait, scale=2)

  The decimal point is at the |

   0 | 0
   2 | 0000000
   4 | 0000000
   6 | 0000
   8 | 00
  10 | 000
  12 | 0

Числови характеристики на данните

x <- airquality$Temp

# средно
mean(x)

[1] 77.88235

# медиана
median(x)

[1] 79

# стандартно отклонение
sd(x)

[1] 9.46527

# p-квантил
quantile(x, p)

50% 
 79 

# Q1
quantile(x, 0.25)

25% 
 72 

# Q3
quantile(x, 0.75)

75% 
 85 

# IQR = Q3 - Q1
IQR(x)

[1] 13

summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  56.00   72.00   79.00   77.88   85.00   97.00 

boxplot(x, horizontal=T)

Многомерни данни

library(MASS)
?survey

str(survey)

'data.frame':   237 obs. of  12 variables:
 $ Sex   : Factor w/ 2 levels "Female","Male": 1 2 2 2 2 1 2 1 2 2 ...
 $ Wr.Hnd: num  18.5 19.5 18 18.8 20 18 17.7 17 20 18.5 ...
 $ NW.Hnd: num  18 20.5 13.3 18.9 20 17.7 17.7 17.3 19.5 18.5 ...
 $ W.Hnd : Factor w/ 2 levels "Left","Right": 2 1 2 2 2 2 2 2 2 2 ...
 $ Fold  : Factor w/ 3 levels "L on R","Neither",..: 3 3 1 3 2 1 1 3 3 3 ...
 $ Pulse : int  92 104 87 NA 35 64 83 74 72 90 ...
 $ Clap  : Factor w/ 3 levels "Left","Neither",..: 1 1 2 2 3 3 3 3 3 3 ...
 $ Exer  : Factor w/ 3 levels "Freq","None",..: 3 2 2 2 3 3 1 1 3 3 ...
 $ Smoke : Factor w/ 4 levels "Heavy","Never",..: 2 4 3 2 2 2 2 2 2 2 ...
 $ Height: num  173 178 NA 160 165 ...
 $ M.I   : Factor w/ 2 levels "Imperial","Metric": 2 1 NA 2 2 1 1 2 2 2 ...
 $ Age   : num  18.2 17.6 16.9 20.3 23.7 ...

# извежда данните като таблица
# fix(survey)

summary(survey)

     Sex          Wr.Hnd          NW.Hnd        W.Hnd          Fold    
 Female:118   Min.   :13.00   Min.   :12.50   Left : 18   L on R : 99  
 Male  :118   1st Qu.:17.50   1st Qu.:17.50   Right:218   Neither: 18  
 NA's  :  1   Median :18.50   Median :18.50   NA's :  1   R on L :120  
              Mean   :18.67   Mean   :18.58                            
              3rd Qu.:19.80   3rd Qu.:19.73                            
              Max.   :23.20   Max.   :23.50                            
              NA's   :1       NA's   :1                                
     Pulse             Clap       Exer       Smoke         Height     
 Min.   : 35.00   Left   : 39   Freq:115   Heavy: 11   Min.   :150.0  
 1st Qu.: 66.00   Neither: 50   None: 24   Never:189   1st Qu.:165.0  
 Median : 72.50   Right  :147   Some: 98   Occas: 19   Median :171.0  
 Mean   : 74.15   NA's   :  1              Regul: 17   Mean   :172.4  
 3rd Qu.: 80.00                            NA's :  1   3rd Qu.:180.0  
 Max.   :104.00                                        Max.   :200.0  
 NA's   :45                                            NA's   :28     
       M.I           Age       
 Imperial: 68   Min.   :16.75  
 Metric  :141   1st Qu.:17.67  
 NA's    : 28   Median :18.58  
                Mean   :20.37  
                3rd Qu.:20.17  
                Max.   :73.00  
                               

# equilavnt notations
survey[, 'Age']

  [1] 18.250 17.583 16.917 20.333 23.667 21.000 18.833 35.833 19.000 22.333
 [11] 28.500 18.250 18.750 17.500 17.167 17.167 19.333 18.333 19.750 17.917
 [21] 17.917 18.167 17.833 18.250 19.167 17.583 17.500 18.083 21.917 19.250
 [31] 41.583 17.500 39.750 17.167 17.750 18.000 19.000 17.917 35.500 19.917
 [41] 17.500 17.083 28.583 17.500 17.417 18.500 18.917 19.417 18.417 30.750
 [51] 18.500 17.500 18.333 17.417 20.000 18.333 17.167 17.417 17.667 18.417
 [61] 20.333 17.333 17.500 19.833 18.583 18.000 30.667 16.917 19.917 18.333
 [71] 17.583 17.833 17.667 17.417 17.750 20.667 23.583 17.167 17.083 18.750
 [81] 16.750 20.167 17.667 17.167 17.167 17.250 18.000 18.750 21.583 17.583
 [91] 19.667 18.000 19.667 17.083 22.833 17.083 19.417 23.250 18.083 19.083
[101] 18.917 17.750 20.833 20.167 17.667 18.250 17.000 18.500 18.583 17.750
[111] 24.167 18.167 21.167 17.917 17.417 20.500 22.917 18.917 18.917 20.083
[121] 17.500 18.250 17.500 17.417 21.000 19.833 17.667 18.083 18.000 18.333
[131] 20.000 18.750 19.083 18.500 18.417 19.167 21.500 19.333 21.417 18.667
[141] 17.500 21.083 17.250 19.000 19.167 19.000 23.000 32.667 20.000 20.167
[151] 25.500 18.167 23.500 70.417 43.833 23.583 21.083 44.250 19.667 17.917
[161] 18.417 21.167 17.500 29.083 19.917 18.500 18.167 32.750 17.417 17.333
[171] 73.000 18.667 18.500 18.667 17.750 17.250 36.583 23.083 19.250 17.167
[181] 23.417 17.083 17.250 23.833 18.750 21.167 24.667 18.500 20.333 20.083
[191] 18.917 27.333 18.917 17.250 18.167 26.500 17.000 17.167 19.167 17.500
[201] 19.250 21.333 18.583 20.167 18.667 17.083 17.417 18.583 19.500 18.500
[211] 17.167 17.250 17.500 20.417 17.083 21.250 19.250 19.333 19.167 18.917
[221] 20.917 17.333 18.167 20.750 19.917 18.667 18.417 17.417 20.333 19.333
[231] 18.167 20.750 17.667 16.917 18.583 17.167 17.750

survey$Age

  [1] 18.250 17.583 16.917 20.333 23.667 21.000 18.833 35.833 19.000 22.333
 [11] 28.500 18.250 18.750 17.500 17.167 17.167 19.333 18.333 19.750 17.917
 [21] 17.917 18.167 17.833 18.250 19.167 17.583 17.500 18.083 21.917 19.250
 [31] 41.583 17.500 39.750 17.167 17.750 18.000 19.000 17.917 35.500 19.917
 [41] 17.500 17.083 28.583 17.500 17.417 18.500 18.917 19.417 18.417 30.750
 [51] 18.500 17.500 18.333 17.417 20.000 18.333 17.167 17.417 17.667 18.417
 [61] 20.333 17.333 17.500 19.833 18.583 18.000 30.667 16.917 19.917 18.333
 [71] 17.583 17.833 17.667 17.417 17.750 20.667 23.583 17.167 17.083 18.750
 [81] 16.750 20.167 17.667 17.167 17.167 17.250 18.000 18.750 21.583 17.583
 [91] 19.667 18.000 19.667 17.083 22.833 17.083 19.417 23.250 18.083 19.083
[101] 18.917 17.750 20.833 20.167 17.667 18.250 17.000 18.500 18.583 17.750
[111] 24.167 18.167 21.167 17.917 17.417 20.500 22.917 18.917 18.917 20.083
[121] 17.500 18.250 17.500 17.417 21.000 19.833 17.667 18.083 18.000 18.333
[131] 20.000 18.750 19.083 18.500 18.417 19.167 21.500 19.333 21.417 18.667
[141] 17.500 21.083 17.250 19.000 19.167 19.000 23.000 32.667 20.000 20.167
[151] 25.500 18.167 23.500 70.417 43.833 23.583 21.083 44.250 19.667 17.917
[161] 18.417 21.167 17.500 29.083 19.917 18.500 18.167 32.750 17.417 17.333
[171] 73.000 18.667 18.500 18.667 17.750 17.250 36.583 23.083 19.250 17.167
[181] 23.417 17.083 17.250 23.833 18.750 21.167 24.667 18.500 20.333 20.083
[191] 18.917 27.333 18.917 17.250 18.167 26.500 17.000 17.167 19.167 17.500
[201] 19.250 21.333 18.583 20.167 18.667 17.083 17.417 18.583 19.500 18.500
[211] 17.167 17.250 17.500 20.417 17.083 21.250 19.250 19.333 19.167 18.917
[221] 20.917 17.333 18.167 20.750 19.917 18.667 18.417 17.417 20.333 19.333
[231] 18.167 20.750 17.667 16.917 18.583 17.167 17.750

# 12-та колона
survey[ ,12]

  [1] 18.250 17.583 16.917 20.333 23.667 21.000 18.833 35.833 19.000 22.333
 [11] 28.500 18.250 18.750 17.500 17.167 17.167 19.333 18.333 19.750 17.917
 [21] 17.917 18.167 17.833 18.250 19.167 17.583 17.500 18.083 21.917 19.250
 [31] 41.583 17.500 39.750 17.167 17.750 18.000 19.000 17.917 35.500 19.917
 [41] 17.500 17.083 28.583 17.500 17.417 18.500 18.917 19.417 18.417 30.750
 [51] 18.500 17.500 18.333 17.417 20.000 18.333 17.167 17.417 17.667 18.417
 [61] 20.333 17.333 17.500 19.833 18.583 18.000 30.667 16.917 19.917 18.333
 [71] 17.583 17.833 17.667 17.417 17.750 20.667 23.583 17.167 17.083 18.750
 [81] 16.750 20.167 17.667 17.167 17.167 17.250 18.000 18.750 21.583 17.583
 [91] 19.667 18.000 19.667 17.083 22.833 17.083 19.417 23.250 18.083 19.083
[101] 18.917 17.750 20.833 20.167 17.667 18.250 17.000 18.500 18.583 17.750
[111] 24.167 18.167 21.167 17.917 17.417 20.500 22.917 18.917 18.917 20.083
[121] 17.500 18.250 17.500 17.417 21.000 19.833 17.667 18.083 18.000 18.333
[131] 20.000 18.750 19.083 18.500 18.417 19.167 21.500 19.333 21.417 18.667
[141] 17.500 21.083 17.250 19.000 19.167 19.000 23.000 32.667 20.000 20.167
[151] 25.500 18.167 23.500 70.417 43.833 23.583 21.083 44.250 19.667 17.917
[161] 18.417 21.167 17.500 29.083 19.917 18.500 18.167 32.750 17.417 17.333
[171] 73.000 18.667 18.500 18.667 17.750 17.250 36.583 23.083 19.250 17.167
[181] 23.417 17.083 17.250 23.833 18.750 21.167 24.667 18.500 20.333 20.083
[191] 18.917 27.333 18.917 17.250 18.167 26.500 17.000 17.167 19.167 17.500
[201] 19.250 21.333 18.583 20.167 18.667 17.083 17.417 18.583 19.500 18.500
[211] 17.167 17.250 17.500 20.417 17.083 21.250 19.250 19.333 19.167 18.917
[221] 20.917 17.333 18.167 20.750 19.917 18.667 18.417 17.417 20.333 19.333
[231] 18.167 20.750 17.667 16.917 18.583 17.167 17.750

# 5-ти ред
survey[5, ]

Крос-таблица

attach(survey)

table(Smoke, W.Hnd)

       W.Hnd
Smoke   Left Right
  Heavy    1    10
  Never   13   175
  Occas    3    16
  Regul    1    16

# премахване на липсващите наблюдения
table(Smoke, W.Hnd, useNA="no")

       W.Hnd
Smoke   Left Right
  Heavy    1    10
  Never   13   175
  Occas    3    16
  Regul    1    16

Сравняване на таблици

# таблица с относителния дял
prop.table(table(Smoke, W.Hnd))

       W.Hnd
Smoke          Left       Right
  Heavy 0.004255319 0.042553191
  Never 0.055319149 0.744680851
  Occas 0.012765957 0.068085106
  Regul 0.004255319 0.068085106

# таблица с редови проценти
prop.table(table(Smoke, W.Hnd), 1)

       W.Hnd
Smoke         Left      Right
  Heavy 0.09090909 0.90909091
  Never 0.06914894 0.93085106
  Occas 0.15789474 0.84210526
  Regul 0.05882353 0.94117647

# таблица с колонни проценти
prop.table(table(Smoke, W.Hnd), 2)

       W.Hnd
Smoke         Left      Right
  Heavy 0.05555556 0.04608295
  Never 0.72222222 0.80645161
  Occas 0.16666667 0.07373272
  Regul 0.05555556 0.07373272

barplot(
    table(Smoke, W.Hnd),
    legend=T,
    args.legend=list(
        x="topleft",
        inset=0.05
    )
)

barplot(
    table(W.Hnd, Smoke),
    beside=T,
    legend=T,
    args.legend=list(
        x="topright",
        inset=0.05
    )
)

Сравняване на стойности на категорийни променливи и корелация

Ако данните са близо до “възходяща” права линия, то корелацията им е близко до 1. Това показва, че те “растат” заедно.

Ако данните са близо до “низходяща” права линия, то корелацията им е близко до -1. Това показва, че те “растат” противоположно.

Ако корелацията им е близо до 0, то не наблюдаваме тенденции на линейна зависимост между тях.

# кутуя с мустаци
boxplot(Pulse ~ Smoke)

# лентова диаграма
stripchart(Pulse ~ Smoke, vertical=T)

# точкова диаграма
plot(Height, Pulse)

cor(Height, Pulse, use="complete.obs")

[1] -0.08394058

plot(Wr.Hnd, Height)

cor(Wr.Hnd, Height, use="complete.obs")

[1] 0.6009909

plot(Wr.Hnd, NW.Hnd)

cor(Wr.Hnd, NW.Hnd, use="complete.obs")

[1] 0.9483103

# преобразуваме данните, за да получим отрицателна корелация
plot(Wr.Hnd, -1 * NW.Hnd)

cor(Wr.Hnd, -1 * NW.Hnd, use="complete.obs")

[1] -0.9483103

Централна гранична теорема

Централна гранична теорема

Централната гранична теорема гласи, че разпределението на средното на извадката е приблизително нормално, ако размерът на извадката е достатъчно голям, дори ако разпределението на популацията не е нормално.

Пример 1

Времето на живот на електрическа крушка от даден тип има експоненциално разпределение със средно 900 часа. Измерено е времето на живот на 100 случайно избрани крушки. Каква е вероятността полученото средно време да е над 980 часа?

Използваме тази формула:

a <- (980-900)/(900/sqrt(100))

# P(X <= 980) = pnorm(a)
# P(X > 980) = 1 - pnorm(a)
1 - pnorm(a)

[1] 0.1870314

Пример 2

Регистрираният багаж на пътниците в даден самолет не трябва да надвишава общо 4000 кг. Количеството регистриран багаж на произволно избран пътник е случайна величина със средно 24 кг и стандартно отклонение 7 кг. Каква е вероятността общото количество регистриран багаж в самолет със 160 пътници да надвиши 4000 кг?

Използваме еквивалентна формула:

a <- (4000 - 160*24)/(7*sqrt(160))

# P(X <= 4000) = pnorm(a)
# P(X > 4000) = 1 - pnorm(a)
1 - pnorm(a)

[1] 0.03537991

Проверка на хипотези при една извадка

z-тест за средно

Използваме z-тест за средно, когато искаме да сравним наблюдаваната средна стойност на случайни величини с предварително зададена стойност. Дисперсията на случайните величини е известна.

Пример

Фирма произвежда електрически крушки. Средното време на живот на една крушка е 2000 часа със стандартно отклонение 300 часа. Предложен е нов тип крушки. Изпробвани са 100 крушки от новия тип. Резултатите показват средно време на живот на новите крушки 2100 часа и същото стандартно отклонение. Може ли да се твърди, че средното време на живот на новия тип крушки е повече от 2000 часа?

# H_0: mu = 2000
# H_1: mu > 2000
x.bar <- 2100
n <- 100
sigma <- 300
mu <- 2000

# пресмятаме наблюдаваната стойност
z.obs <- (x.bar - mu) / (sigma/sqrt(n))

# P(X <= 2000) = pnorm(z.obs)
# P(X > 2000) = 1 - pnorm(z.obs)
p.value <- 1 - pnorm(z.obs)

# p-стойността е по-малка от 0,05
# -> отхвърляме H_0 в полза на H_1]
# -> времето на живот е повече от 2000
p.value

[1] 0.0004290603

t-тест за средно

Използваме t-тест за средно, когато искаме да сравним наблюдаваната средна стойност на случайни величини с предварително зададена стойност. Дисперсията на случайните величини е неизвестна.

Пример

Според исторически данни, средната киселинност на дъждовете в определен индустриален район е 5.2. За да се провери има ли изменение в тази стойност е измерена киселинността на 12 валежа през изминалата година. Получени са следните резултати:

6.1   5.4   4.8   5.8   6.6   5.3   6.1   4.4   3.9   6.8   6.5   6.3

От предишни изследвания е известно, че киселинността на валежите има нормално разпределение. Имаме ли основания да твърдим, че киселинността в района се е променила в сравнение с историческите данни?

# H_0: mu = 5.2
# H_1: mu =\= 5.2
x <- c(6.1, 5.4, 4.8, 5.8, 6.6, 5.3, 6.1, 4.4, 3.9, 6.8, 6.5, 6.3)

t.test(x, mu=5.2)

    One Sample t-test

data:  x
t = 1.7556, df = 11, p-value = 0.1069
alternative hypothesis: true mean is not equal to 5.2
95 percent confidence interval:
 5.081616 6.251717
sample estimates:
mean of x 
 5.666667 

# p-стойността е по-голяма от 0.05
# -> не отхвърляме H_0
# -> средната стойност не се е променила
t.test(x, mu=5.2)$p.value

[1] 0.1069226

z-тест за пропорция

Използваме z-тест за пропроция, когато имаме Бернулиеви опити (с резултат “успех” или “неуспех”), всеки от които с вероятност за успех p. Предполагаме, че имаме “успех” в X% от случаите. Искаме да проверим дали предположението ни отговаря на действителността.

Пример

Известно е, че при 10% от колите от дадена марка се появява сериозен дефект по време на гаранционния срок. От първите 25000 продадени коли от нов модел, 2700 се оказали с дефект. Може ли да се твърди, че вероятността в кола от новия модел да се появи дефект не е 10%?

# H_0: p = 0.1
# H_1: p =\= 0.1
prop.test(x=2700, n=25000, p=0.1, correct=F)

    1-sample proportions test without continuity correction

data:  2700 out of 25000, null probability 0.1
X-squared = 17.778, df = 1, p-value = 2.483e-05
alternative hypothesis: true p is not equal to 0.1
95 percent confidence interval:
 0.1042126 0.1119078
sample estimates:
    p 
0.108 

# p-стойността е по-малка от 0.05
# -> отхвърляме H_0 в полза на H_1
# -> вероятността не е 10%
prop.test(x=2700, n=25000, p=0.1, correct=F)$p.value

[1] 2.482661e-05

Проверка на хипотези при две извадки

t-тест за разлика на средни

Използваме t-тест за разлика на средни, когато имаме 2 независими извадки, чиито средни искаме да сравним.

Пример

Мениджър обмисля въвеждане на допълнителна 15-минутна почивка за работниците си. За да разбере дали такава почивка ще намали броя на грешките, които правят работниците, избрал случайно 2 групи по 10 души. Първата група имала допълнителна почивка, а втората работила по обичайното работно време. Данните за броя на допуснатите грешки от двете групи са следните:

Група 1: 8, 7, 5, 8, 10, 9, 7, 8, 4, 5
Група 2: 7, 6, 14, 12, 13, 8, 9, 6, 10, 9

Приемаме, че данните са приблизително нормално разпределени. Може ли да се заключи, че работниците с допълнителна почивка правят средно по-малко грешки?

# H_0: mu_X = mu_Y
# H_1: mu_X < mu_Y
x <- c(8, 7, 5, 8, 10, 9, 7, 8, 4, 5)
y <- c(7, 6, 14, 12, 13, 8, 9, 6, 10, 9)


t.test(x, y, alt="less")

    Welch Two Sample t-test

data:  x and y
t = -2.1264, df = 15.78, p-value = 0.02481
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
       -Inf -0.4099189
sample estimates:
mean of x mean of y 
      7.1       9.4 

# p-стойността е по-малка от 0.05
# -> отхвърляме H_0 в полза на H_1
# -> средното на група 1 е по-малко от това на група 2
t.test(x, y, alt="less")$p.value

[1] 0.02480789

t-тест при зависими извадки

t-тест, който използваме, когато имаме 2 зависими извадки (т.е. на всеки елемент от 1-вата извадка съответства елемент от втората - двойки наблюдения), чиито средни искаме да сравним.

Пример

За да се изследва ефекта на диета върху нивото на холестерол в кръвта са избрани 15 мъже на възраст между 35 и 50 години. Нивото на холестерола на всеки участник е измерено преди започване на диетата и три месеца след прилагане на диетата. Данните са следните:

Участник Преди След
       1   265  229
       2   240  231
       3   258  227
       4   295  240
       5   251  238
       6   245  241
       7   287  234
       8   314  256
       9   260  247
      10   279  239
      11   283  246
      12   240  218
      13   238  219
      14   225  226
      15   247  233

Приемаме, че нивото на холестерол е нормално разпределено. Дали тези данни дават основание да се твърди, че диетата намалява нивото на холестерол в средно?

# H_0: mu_X = mu_Y
# H_1: mu_X > mu_Y
x <- c(265,240,258,295,251,245,287,314,260,279,283,240,238,225,247)
y <- c(229,231,227,240,238,241,234,256,247,239,246,218,219,226,233)

t.test(x, y, alt="greater", paired=T)

    Paired t-test

data:  x and y
t = 5.4659, df = 14, p-value = 4.158e-05
alternative hypothesis: true mean difference is greater than 0
95 percent confidence interval:
 18.20922      Inf
sample estimates:
mean difference 
       26.86667 

# p-стойността е по-малка от 0.05
# -> отхвърляме H_0 в полза на H_1
# -> средното "Преди" е по-малко от това "След"
t.test(x, y, alt="greater", paired=T)$p.value

[1] 4.157964e-05

z-тест за разлика на пропорции

Използваме z-тест за разлика на пропорции, когато имаме 2 независими извадки, показващи резултати от Бернулиеви опити (т.е. всеки един елемент от извадката да е булева стойност), и искаме да сравним пропорциите на “успех” в двете извадки.

Пример

В проучване участвали 220 жени и 210 мъже. Според резултатите, 71 жени и 58 мъже отговорили, че предпочитат безкофеиново кафе. Може ли да твърдим, че процентът на жените, предпочитащи безкофеиново кафе, е различен от процентът на мъжете, предпочитащи безкофеиново кафе?

# H_0: p_X = p_Y
# H_1: p_X =\= p_Y

# брой "успехи" във всяка от извадките
x <- c(71,58)
# съответстващ общ брой наблюдения във всяка от извадките
n <- c(220,210)

prop.test(x, n, correct=F)

    2-sample test for equality of proportions without continuity correction

data:  x out of n
X-squared = 1.108, df = 1, p-value = 0.2925
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.03991226  0.13298585
sample estimates:
   prop 1    prop 2 
0.3227273 0.2761905 

# p-стойността е по-голяма от 0.05
# -> не отхвърляме H_0
# -> процентът на жените, пиещи безкофеиново кафе, не е по-голям от този на мъжете
prop.test(x, n, correct=F)$p.value

[1] 0.292513

Доверителни интервали

Доверителен интервал за средно при известна дисперсия

Намираме доверителен интервал за средното независими и еднакво разпределени случайни величини. Дисперсията трябва да е известна.

Пример

Фирма произвежда електрически крушки. Средното време на живот на една крушка е 2000 часа със стандартно отклонение 300 часа. Предложен е нов тип крушки. Изпробвани са 100 крушки от новия тип. Резултатите показват средно време на живот на новите крушки 2100 часа и същото стандартно отклонение. Намерете 95-процентен доверителен интервал за средното време на живот на новия тип крушки.

z1.ci <- function(x.bar, sigma, n, alpha) {
b1 <- x.bar - qnorm(1-alpha/2)*(sigma/sqrt(n))
b2 <- x.bar + qnorm(1-alpha/2)*(sigma/sqrt(n))
c(b1, b2)
}

z1.ci(x.bar=2100, sigma=300, n=100, alpha=0.05)

[1] 2041.201 2158.799

Доверителен интервал за средно при неизвестна дисперсия

Намираме доверителен интервал за средното независими и еднакво разпределени случайни величини. Дисперсията трябва да е неизвестна.

Пример

Според исторически данни, средната киселинност на дъждовете в определен индустриален район е 5.2. За да се провери има ли изменение в тази стойност е измерена киселинността на 12 валежа през изминалата година. Получени са следните резултати:

6.1   5.4   4.8   5.8   6.6   5.3   6.1   4.4   3.9   6.8   6.5   6.3

Намерете 99-процентен доверителен интервал за средната киселинност.

x <- c(6.1, 5.4, 4.8, 5.8, 6.6, 5.3, 6.1, 4.4, 3.9, 6.8, 6.5, 6.3)

t.test(x, conf.level=0.99)$conf.int[1:2]

[1] 4.841103 6.492230

Доверителен интервал за пропорция

Намираме доверителен интервал за процента успехи при Бернулиеви опити.

Пример

От първите 25000 продадени коли от нов модел, 2700 се оказали с дефект. Намерете 95-процентен доверителен интервал за вероятността кола от новия модел да е дефектна.

prop.test(x=2700, n=25000, conf.level=0.95, correct=F)$conf.int[1:2]

[1] 0.1042126 0.1119078

Доверителен интервал за медиана

Пример

Направено е проучване с цел да се изследва ефектът на звука от сърдечния ритъм на майката върху новороденото. Бебетата в родилно отделение са разделени на две групи. Първата група е непрекъснато изложена на звука от сърдечния ритъм на възрастен, а втората група не е изложена на такъв звук. Измерена е промяната в теглото на бебетата от раждането до четвъртия ден. Данните са във файла salk.txt . Намерете 95-процентен доверителен интервал за медианата на промяната в теглото на бебетата (за първата и за втората група).

med1.ci <- function(x, alpha=0.05, nboot=1000) {
x <- x[is.finite(x)]
nx <- length(x)
est1 <- median(x)
med1.bt <- replicate( nboot, median( sample( x, size=nx, replace=TRUE ) ) )
ci <- quantile( med1.bt, probs=c(alpha/2, 1-alpha/2), names=FALSE )
list( est.med1=est1, ci=ci )
}

Доверителен интервал за разлика на медиани

Пример

Направено е проучване с цел да се изследва ефектът на звука от сърдечния ритъм на майката върху новороденото. Бебетата в родилно отделение са разделени на две групи. Първата група е непрекъснато изложена на звука от сърдечния ритъм на възрастен, а втората група не е изложена на такъв звук. Измерена е промяната в теглото на бебетата от раждането до четвъртия ден. Данните са във файла salk.txt . Намерете 95-процентен доверителен интервал за медианата на промяната в теглото на бебетата (за първата и за втората група).

med2.ci <- function(x, y, alpha=0.05, nboot=1000) {
    x <- x[is.finite(x)]
    y <- y[is.finite(y)]
    nx <- length(x)
    ny <- length(y)
    est1 <- median(x)
    est2 <- median(y)
    est.dif <- est1 - est2
    med1.bt <- replicate( nboot, median( sample( x, size=nx, replace=TRUE ) ) )
    med2.bt <- replicate( nboot, median( sample( y, size=ny, replace=TRUE ) ) )
    dif.bt <- med1.bt - med2.bt
    ci <- quantile( dif.bt, probs=c(alpha/2, 1-alpha/2), names=FALSE )
    list( est.med1=est1, est.med2=est2, est.dif=est.dif, ci=ci )
}