Metody probabilistyczne i statytyka - laboratorium (realizacja)

¹ Zakład Bioinformatyki, Instytut Informatyki, Uniwersytet w Białymstoku

^✉ Correspondence: Jarosław Kotowicz <j.kotowicz@uwb.edu.pl>

rm(list = ls())

1 Laboratorium z dnia 23 kwietnia.

2 Powtórzenie i utrwalenie wiadomości odnośnie testowania hipotez.

set.seed(20200423)
x <- rnorm(100, mean = 3, sd = 5)
set.seed(20200423)
y <- runif(100, -2, 2)
set.seed(20200423)
z <- rnorm(100, mean = 1, sd = 2)

2.1 Test t.test (średnie)

t.test(x)

    One Sample t-test

data:  x
t = 8.9586, df = 99, p-value = 2.039e-14
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 3.245063 5.091516
sample estimates:
mean of x 
 4.168289 

Interpretacja wyniku!

t.test(x, mu = 3)

    One Sample t-test

data:  x
t = 2.5109, df = 99, p-value = 0.01366
alternative hypothesis: true mean is not equal to 3
95 percent confidence interval:
 3.245063 5.091516
sample estimates:
mean of x 
 4.168289 

Interpretacja wyniku!

2.1.1 Różne hioptezy alternatywe

t.test(x, mu = 3, alternative = "less")

    One Sample t-test

data:  x
t = 2.5109, df = 99, p-value = 0.9932
alternative hypothesis: true mean is less than 3
95 percent confidence interval:
     -Inf 4.940844
sample estimates:
mean of x 
 4.168289 

Interpretacja wyniku!

t.test(x, mu = 3, alternative = "greater")

    One Sample t-test

data:  x
t = 2.5109, df = 99, p-value = 0.006831
alternative hypothesis: true mean is greater than 3
95 percent confidence interval:
 3.395734      Inf
sample estimates:
mean of x 
 4.168289 

Interpretacja wyniku!

2.1.2 Dwie średnie

t.test(x[1:50], x[51:100])

    Welch Two Sample t-test

data:  x[1:50] and x[51:100]
t = 0.16905, df = 97.965, p-value = 0.8661
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.697727  2.013916
sample estimates:
mean of x mean of y 
 4.247336  4.089242 

Interpretacja wyniku!

2.2 Testy nieparametryczne zgodności z rozkładem normalnym w bibliotece nortest

library(nortest)

ad.test(x)

    Anderson-Darling normality test

data:  x
A = 0.33084, p-value = 0.5095

Interpretacja wyniku!

cvm.test(x)

    Cramer-von Mises normality test

data:  x
W = 0.051848, p-value = 0.4809

Interpretacja wyniku!

lillie.test(x)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  x
D = 0.066694, p-value = 0.3357

Interpretacja wyniku!

pearson.test(x)

    Pearson chi-square normality test

data:  x
P = 19.86, p-value = 0.03061

Interpretacja wyniku!

3 Testowanie hipotez nieparametrycznych

3.1 Test Kołogorowa-Smirnowa (testowanie zgodności z rozkładem wzorocowym)

3.1.1 Rozkład normalny

Niepoprawny sposób wywołania! Testuje z domyślnymi wartościami parametrów!

ks.test(x, "pnorm")

    One-sample Kolmogorov-Smirnov test

data:  x
D = 0.64989, p-value < 2.2e-16
alternative hypothesis: two-sided

Poprawne sposób wywołania!

ks.test(x, "pnorm", 3, 5)

    One-sample Kolmogorov-Smirnov test

data:  x
D = 0.11549, p-value = 0.1388
alternative hypothesis: two-sided

Interpretacja wyniku!

mean(x)

[1] 4.168289

sd(x)

[1] 4.65285

ks.test(x, "pnorm", 4.1, 4.6)

    One-sample Kolmogorov-Smirnov test

data:  x
D = 0.062321, p-value = 0.8321
alternative hypothesis: two-sided

Interpretacja wyniku!

3.1.2 Rozkład jednostajny

ks.test(y, "punif", -2, 2)

    One-sample Kolmogorov-Smirnov test

data:  y
D = 0.10912, p-value = 0.1847
alternative hypothesis: two-sided

Interpretacja wyniku!

3.1.3 Rozkład beta

set.seed(20200423)
v <- rbeta(100, 2, 1)

ks.test(v, "pbeta", 2, 1)

    One-sample Kolmogorov-Smirnov test

data:  v
D = 0.11789, p-value = 0.1241
alternative hypothesis: two-sided

Interpretacja wyniku!

3.2 Test Kołogorowa-Smirnowa (testowanie zgodności dwóch rozkładów)

ks.test(x, z)

    Two-sample Kolmogorov-Smirnov test

data:  x and z
D = 0.42, p-value = 4.366e-08
alternative hypothesis: two-sided

Interpretacja wyniku!

ks.test(x, z, alternative = "less")

    Two-sample Kolmogorov-Smirnov test

data:  x and z
D^- = 0.42, p-value = 2.183e-08
alternative hypothesis: the CDF of x lies below that of y

Interpretacja wyniku!

ks.test(x, z, alternative = "greater")

    Two-sample Kolmogorov-Smirnov test

data:  x and z
D^+ = 0.08, p-value = 0.5273
alternative hypothesis: the CDF of x lies above that of y

Interpretacja wyniku!

3.3 Test chi-kwadrat zgodności z rozkładem jednostajnym

Trzeba odpowiednio przygotować dane!

chisq.test(x)

Błąd w poleceniu 'chisq.test(x)':
  wszystkie wpisy 'x' muszą być nieujemne oraz skończone

Wstęp do kubełkowania, czyli jak dobrze doprać punkty podziału zbioru wartości próbki!

seq(min(x)-.1, max(x)+.1, by = 1)

 [1] -7.7008496 -6.7008496 -5.7008496 -4.7008496 -3.7008496 -2.7008496 -1.7008496 -0.7008496  0.2991504  1.2991504
[11]  2.2991504  3.2991504  4.2991504  5.2991504  6.2991504  7.2991504  8.2991504  9.2991504 10.2991504 11.2991504
[21] 12.2991504 13.2991504 14.2991504 15.2991504

Kubełkowanie, czyli jak ze zmiennej ciągłej zrobić zmienną czynnikową

cut(x, seq(min(x)-.1, max(x)+1, by = 1))

  [1] (1.3,2.3]      (-0.701,0.299] (5.3,6.3]      (-0.701,0.299] (14.3,15.3]    (7.3,8.3]      (-2.7,-1.7]   
  [8] (1.3,2.3]      (6.3,7.3]      (1.3,2.3]      (-0.701,0.299] (7.3,8.3]      (-2.7,-1.7]    (10.3,11.3]   
 [15] (-2.7,-1.7]    (13.3,14.3]    (15.3,16.3]    (9.3,10.3]     (3.3,4.3]      (10.3,11.3]    (3.3,4.3]     
 [22] (6.3,7.3]      (3.3,4.3]      (11.3,12.3]    (1.3,2.3]      (2.3,3.3]      (6.3,7.3]      (-1.7,-0.701] 
 [29] (7.3,8.3]      (2.3,3.3]      (4.3,5.3]      (6.3,7.3]      (3.3,4.3]      (8.3,9.3]      (2.3,3.3]     
 [36] (0.299,1.3]    (1.3,2.3]      (6.3,7.3]      (0.299,1.3]    (-2.7,-1.7]    (1.3,2.3]      (8.3,9.3]     
 [43] (-0.701,0.299] (-1.7,-0.701]  (1.3,2.3]      (5.3,6.3]      (1.3,2.3]      (5.3,6.3]      (2.3,3.3]     
 [50] (-1.7,-0.701]  (5.3,6.3]      (4.3,5.3]      (7.3,8.3]      (-0.701,0.299] (-0.701,0.299] (4.3,5.3]     
 [57] (2.3,3.3]      (10.3,11.3]    (13.3,14.3]    (0.299,1.3]    (2.3,3.3]      (3.3,4.3]      (12.3,13.3]   
 [64] (9.3,10.3]     (-5.7,-4.7]    (4.3,5.3]      (-2.7,-1.7]    (6.3,7.3]      (-2.7,-1.7]    (3.3,4.3]     
 [71] (5.3,6.3]      (7.3,8.3]      (0.299,1.3]    (2.3,3.3]      (0.299,1.3]    (-7.7,-6.7]    (6.3,7.3]     
 [78] (7.3,8.3]      (2.3,3.3]      (11.3,12.3]    (1.3,2.3]      (4.3,5.3]      (3.3,4.3]      (1.3,2.3]     
 [85] (5.3,6.3]      (5.3,6.3]      (-1.7,-0.701]  (-0.701,0.299] (9.3,10.3]     (5.3,6.3]      (2.3,3.3]     
 [92] (-2.7,-1.7]    (-5.7,-4.7]    (6.3,7.3]      (10.3,11.3]    (7.3,8.3]      (1.3,2.3]      (7.3,8.3]     
 [99] (2.3,3.3]      (5.3,6.3]     
24 Levels: (-7.7,-6.7] (-6.7,-5.7] (-5.7,-4.7] (-4.7,-3.7] (-3.7,-2.7] (-2.7,-1.7] ... (15.3,16.3]

xx <- table(cut(x, seq(min(x)-.1, max(x)+1, by = 1)))

chisq.test(xx)

Aproksymacja chi-kwadrat mo戼㹦e by攼㸶 niepoprawna

    Chi-squared test for given probabilities

data:  xx
X-squared = 65.12, df = 23, p-value = 6.727e-06

Interpretacja wyniku!

3.4 Testowanie zgodność z rozkładem dyskretnym

3.4.1 Biblioteka vcd

library(vcd)

愼㸳adowanie wymaganego pakietu: grid

set.seed(20200423)
a <- rbinom(100, size = 1, prob = .3)

goodfit(a,  type = "binomial", method = "MinChisq", par = list(prob = .3, size = 1))

Observed and fitted values for binomial distribution
with fixed parameters 

 count observed fitted pearson residual
     0       60     70        -1.195229
     1       40     30         1.825742

Interpretacja wyniku!

3.4.2 Test chisq.test

chisq.test(a)

Aproksymacja chi-kwadrat mo戼㹦e by攼㸶 niepoprawna

    Chi-squared test for given probabilities

data:  a
X-squared = 60, df = 99, p-value = 0.9993

Interpretacja wyniku!

df <- cbind(x, a)
colnames(df) <- c("wartosci", "grupa")

df %>% head(10)

        wartosci grupa
 [1,]  1.9370821     0
 [2,] -0.3536600     0
 [3,]  5.7540513     0
 [4,] -0.5273656     1
 [5,] 14.9846370     1
 [6,]  7.4107111     0
 [7,] -2.1444124     0
 [8,]  1.9466011     1
 [9,]  6.7705272     1
[10,]  1.4558207     1

4 Laboratorium z dnia 30 kwietnia.

4.1 Wyznaczanie parametrów rozkładu metodą najwięszej wiarogodności

library(MASS)

Do戼㸳戼㸹czanie pakietu: 㤼㸱MASS㤼㸲

Nast攼㹡puj戼㸹cy obiekt zosta戼㸳 zakryty z 㤼㸱package:dplyr㤼㸲:

    select

Nast攼㹡puj戼㸹cy obiekt zosta戼㸳 zakryty z 㤼㸱dane㤼㸲:

    Insurance

set.seed(20200430)
x.norm <- rnorm(100, mean = 1, sd =2)
set.seed(20200430)
x.lnorm <- rlnorm(100, meanlog = .1, sdlog = 2)

set.seed(20200430)
x.gamma <- rgamma(100, 2, 3)

4.1.1 Funkcja x.fitt

x.fitt <- fitdistr(x.norm, "normal")
x.fitt$estimate

     mean        sd 
0.8082532 1.8567991 

y.fitt <- fitdistr(x.lnorm, "log-normal")
y.fitt$estimate

    meanlog       sdlog 
-0.09174678  1.85679914 

(y.fitt <- fitdistr(x.lnorm, "log-normal"))

     meanlog        sdlog   
  -0.09174678    1.85679914 
 ( 0.18567991) ( 0.13129553)

x.fitt <- fitdistr(x.lnorm, "gamma", list(shape = 3, rate = 2))

wyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaN

x.fitt

     shape         rate   
  0.36923824   0.06654028 
 (0.04194398) (0.01330221)

4.2 Porównywanie średnich (więcej niż dwóch) czyli ANOVA

4.2.1 Podczytanie danych

library(readr)

daneMieszkania <- read_delim("http://www.biecek.pl/R/dane/daneMieszkania.csv", 
                             ";", escape_double = FALSE, trim_ws = TRUE)

Parsed with column specification:
cols(
  cena = col_double(),
  pokoi = col_double(),
  powierzchnia = col_double(),
  dzielnica = col_character(),
  `typ budynku` = col_character()
)

daneMieszkania

4.2.2 Zmiana zmiennych dzielnica i typ budynku na zmienne czynnikowe

library(tidyverse)

daneMieszkania <- daneMieszkania %>%
  mutate_if(is.character, list(factor))

4.2.3 Podsumowanie obserwacji

daneMieszkania %>% summary

      cena            pokoi       powierzchnia         dzielnica      typ budynku
 Min.   : 83280   Min.   :1.00   Min.   :17.00   Biskupin   :65   kamienica :61  
 1st Qu.:143304   1st Qu.:2.00   1st Qu.:31.15   Krzyki     :79   niski blok:63  
 Median :174935   Median :3.00   Median :43.70   Srodmiescie:56   wiezowiec :76  
 Mean   :175934   Mean   :2.55   Mean   :46.20                                   
 3rd Qu.:208741   3rd Qu.:3.00   3rd Qu.:61.40                                   
 Max.   :295762   Max.   :4.00   Max.   :87.70                                   

colnames(daneMieszkania)[1] <- "cena"

4.2.4 Budowa modelu liniowego (regresja)

model <- lm(cena~dzielnica, data = daneMieszkania)

4.2.5 Sprawdzanie róWności średnich w podgrupach

anova(model)

Analysis of Variance Table

Response: cena
           Df     Sum Sq    Mean Sq F value   Pr(>F)   
dzielnica   2 1.7995e+10 8997691613  5.0456 0.007294 **
Residuals 197 3.5130e+11 1783263361                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretacja wyniku!

attach(daneMieszkania)

4.2.6 Wykres pudełko-wąsy

boxplot(cena~dzielnica)