Laboratorium 2

Podstawowe rozkłady zmiennych losowych

Rozpatrzmy jednowymiarowe zmienne losowe. Nazewnictwo funkcji związanych ze zmiennymi losowymi jest zestandaryzowane. Nazwy funkcji opisane są zgodnie ze schematem.

[prefix] [nazwa rodziny rozkl] ()

prefixy

• Litera r (od random) - generuje próbę prostą o liczebności n. Gdzie n to pierwszy argument funkcji z określonego rozkładu.
• Litera p (od probability) - wartość dystrybuanty danego rozkładu w punktach określonych przez wektor x (pierwszy arg funkcji)
• Litera d (od density) - wartość prawdopodobieństwa (bądź gęstości prawdopodobieństwa) danego rozkładu w punktach określoncyh przez wektor x. -Litera q (od quantile) - wyznacza wartości kwantyli danego rozkładu.

Rozkłady

# rozkład jednostajny
runif(5)

## [1] 0.9913186 0.1564684 0.2823996 0.1561200 0.5552902

# rozkład normalny
rnorm(10)

##  [1] -0.2707108 -0.6530515 -0.2486439 -0.9750074 -0.2680999  1.2086921
##  [7] -0.8203182  0.4834951 -0.8819339  1.3057377

rnorm(10, mean = 4, sd = 2)

##  [1] 5.963445 4.843093 3.096795 2.073831 3.030533 6.704087 1.620970
##  [8] 5.953806 1.181132 2.639664

Zgodnie z regułą 3 sigma prawdopodobieństwo, że zmienna losowa o standardowym rozkładzie normalnym przyjmie wartość mniejszą niż -3 lub większą od 3 jest bardzo małe. Prawdopodobieństwo to można wyznaczyć z dystrybuanty

pnorm(-3) + (1-pnorm(3))

## [1] 0.002699796

Gęstość rozkładu w punktach -1,0,1 otrzymamy z:

dnorm(-1:1, mean=0, sd=1)

## [1] 0.2419707 0.3989423 0.2419707

Kwantyle dla rozkładu normalnego:

qnorm(c(0.001,0.025,0.05,0.5,0.95,0.975,0.999))

## [1] -3.090232 -1.959964 -1.644854  0.000000  1.644854  1.959964  3.090232

Inne rozkłady:

• rbeta Rozkład Beta
• rchisq Rozkład Chi-kwadrat
• rgamma Rozkład Gamma
• rlnorm Rozkład Log normalny
• rt Rozkład t-studenta
• rexp Rozkład Wykładniczy
• rbinom Rozkład Dwumianowy
• rnbinom Rozkład Ujemny dwumianowy
• rpois Rozkład Poissona

Estymacja parametrów rozkładu

Esytmację parametrów wybranego rozkładu otrzymujemy korzystając z funkcji fitdistr()

test <- rlnorm(100)
test[1:10]

##  [1] 2.5490572 3.6282405 1.5514653 1.9711813 0.7476180 1.5169339 5.0771396
##  [8] 1.4182751 0.7468536 2.2687698

MASS::fitdistr(test, "normal")

##      mean         sd    
##   1.7812724   1.7022056 
##  (0.1702206) (0.1203641)

MASS::fitdistr(test, "gamma")

##      shape       rate   
##   1.3069499   0.7337154 
##  (0.1664717) (0.1133715)

Statystyki opisowe

Zmienne w analizie danych można scharakteryzować najczęściej jako:

1. Zmienne jakościowe - przyjmujące określoną liczbę wartości (najczęściej nieliczbowych)

• binarne
• nieuporządkowane (nominalne)
• uporządkowane (porządkowe)

2. Zmienne ilościowe

W R zmienne ilościowe reprezentowane są przez typ numeric. Zmienne jakościowe, natomiast przez typ factor. Czasem zmienne jakościowe binarne reprezentowane są przez typ logica.

Zmienne ilościowe

Dane do analizy

library("Przewodnik")

# define random vector
x <- rnorm(1000, mean = 10, sd = 4)
# data frame
head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

# others df
head(daneSoc)

##   wiek wyksztalcenie st.cywilny      plec             praca
## 1   70      zawodowe  w zwiazku mezczyzna uczen lub pracuje
## 2   66      zawodowe  w zwiazku   kobieta uczen lub pracuje
## 3   71      zawodowe    singiel   kobieta uczen lub pracuje
## 4   57       srednie  w zwiazku mezczyzna uczen lub pracuje
## 5   45       srednie  w zwiazku   kobieta uczen lub pracuje
## 6   48       srednie  w zwiazku mezczyzna       nie pracuje
##   cisnienie.skurczowe cisnienie.rozkurczowe
## 1                 143                    83
## 2                 123                    80
## 3                 167                    80
## 4                 150                    87
## 5                 130                    83
## 6                 138                    75

1. Jaka jest najmniejsza i największa wartość przyjmowana przez zmienną ? Wartości skrajne mogą pomagać wykrywać błędy danych.

min(x)

## [1] -3.930409

max(x)

## [1] 20.04522

range(x)

## [1] -3.930409 20.045217

c(min(x),max(x))

## [1] -3.930409 20.045217

range(iris$Sepal.Length)

## [1] 4.3 7.9

range(daneSoc$wiek)

## [1] 22 75

2. Średnia zmiennej

Zdefiniujmy własną funkcję

my.mean <- function(x){
    return(sum(x)/length(x))
}

i obliczmy średnią za pomocą funkcji z R i własną.

# wynik z R
mean(x)

## [1] 9.864133

# wynik z naszej funkcji
my.mean(x)

## [1] 9.864133

mean(iris$Sepal.Length)

## [1] 5.843333

mean(daneSoc$wiek)

## [1] 43.16176

3. Mediana - jeżli mediana jest bliska średniej to ?

median(x)

## [1] 9.803442

median(iris$Sepal.Length)

## [1] 5.8

median(daneSoc$wiek)

## [1] 45

Dla odważnych

my.median <- function(x) { sorted.x <- sort(x)
if (length(x) %% 2 == 0) {
indices <- c(length(x) / 2, length(x) / 2 + 1)
return(mean(sorted.x[indices])) }
else {
index <- ceiling(length(x) / 2)
return(sorted.x[index]) }
}

4. Odchylenie standardowe

Zdefiniujmy wariancje

my.var <- function(x) { m <- mean(x)
return(sum((x - m) ^ 2) / length(x)) }

Ale czy to dobry wzór ?

my.var <- function(x) { m <- mean(x)
return(sum((x - m) ^ 2) / (length(x) - 1)) }

Jak zdefiniować odchylenie ?

my.sd <- function(x) { 
    return(sqrt(my.var(x)))
}

sd(x)

## [1] 4.079838

sd(iris$Sepal.Length)

## [1] 0.8280661

sd(daneSoc$wiek)

## [1] 13.8471

5. Kurtoza (miara spłaszczenia)

6. Podsumowanie

summ.x <- summary(x)
summ.iris <- summary(iris$Sepal.Length)
summ.wiek <- summary(daneSoc$wiek)

Zadanie !!! Z otrzymanego podsumowania wybierz piątkę Tukeya.

Inne przydatne

# ilosc elementow
length(x)

## [1] 1000

# rozstęp miedzykwartylowy q0.75-q0.25
IQR(x)

## [1] 5.344031

# variancja
var(x)

## [1] 16.64508

# kwantyle
quantile(x, c(0.1,0.25,0.5,0.75,0.9))

##       10%       25%       50%       75%       90% 
##  4.576085  7.325584  9.803442 12.669615 14.963532

Zmienne jakościowe

Table - w przypadku wartości NA pomija w zliczaniu.

Summary - w przypadku wartości NA pokazuje ich ilość.

table(daneSoc$wyksztalcenie)

## 
## podstawowe    srednie     wyzsze   zawodowe 
##         93         55         34         22

table(iris$Species)

## 
##     setosa versicolor  virginica 
##         50         50         50

# dla wielu zmiennych jakościowych
table(daneSoc$wyksztalcenie, daneSoc$praca)

##             
##              nie pracuje uczen lub pracuje
##   podstawowe          22                71
##   srednie             16                39
##   wyzsze               6                28
##   zawodowe             8                14

summary(daneSoc$wyksztalcenie)

## podstawowe    srednie     wyzsze   zawodowe 
##         93         55         34         22

summary(iris$Species)

##     setosa versicolor  virginica 
##         50         50         50

x <- rnorm(1000, mean = 10, sd = 4)

Analiza graficzna

barplot

Funkcja barplot() pozwala przedstawić dane w postaci słupków. Za argument możemy podać wektor lub macierz.

Przydatne argumenty:

• horiz - pion lub poziom dla słupków
• las - kierunek opisów
• funkcja legend() dodaje legendę
• col - kolory słupków

tab <- table( daneSoc$wyksztalcenie)
barplot(tab, horiz = TRUE, las=1)

tab2 <- table( daneSoc$plec, daneSoc$wyksztalcenie)
barplot(tab2, las=1, beside = TRUE, col=c('blue','red'))
legend("topright",c('kobieta','mezczyzna'),fill=c('blue','red'))

Histogram

Histogram generowany jest za pomocą funkcji hist() - przedstawia rozkład wartości zmiennych ilościowych. Argumentem podstawowym funkcji jest wektor liczb.

hist(x)

• breaks - określa sposób podziału zakresu zmienności
• probability - pozwala wyświetlić frakcje zamiast liczebności w przedziałach

hist(daneSoc$wiek, breaks = 15, main='Histogram zm wiek', las=1,ylab='Liczba', xlab="wiek")

hist(daneSoc$wiek, breaks = 25, main='Histogram zm wiek', las=1, probability=TRUE, ylab='Czestosc', xlab="wiek", col="grey", border = "white" )

Pakiet ggplot

ggplot(daneSoc, aes(x = wiek))+geom_histogram(binwidth = 5) + labs(x = "Wiek", y = "Liczba")

ggplot(daneSoc, aes(x = wiek))+geom_histogram(binwidth = 1) + labs(x = "Wiek", y = "Liczba")

Boxplot

Do porównania rozkładów wartości w grupach bardzo często stosuje się wykres pudełkowy boxplot. Do jego wyznaczenia (dla zmiennej ilościowej) służy funkcja boxplot().

Głównym argumentem może być wektor ale również i zbiór wektorów.

boxplot(x)

y=rnorm(1000, mean = 0, sd =2)
boxplot(x,y)

boxplot(daneSoc$cisnienie.rozk, daneSoc$cisnienie.skur, horizontal = TRUE, names = c('Skurczowe','Rozkurczowe'))

Argumentem głównym może być również formuła określająca sposób grupowania.

boxplot(wiek~wyksztalcenie, data=daneSoc, ylab='Wiek', las=1,col='lightgrey', varwidth=TRUE)

Rozkład gęstości

Gładka wersja histogramu generowana jest przy pomocy funkcji density(). Argumentem głównym jest wektor liczb dla których chcemy otrzymać jądrowy estymator gęstości.

plot(density(daneSoc$wiek), main="Rozklad wieku")

Za pomocą parametru bw możemy regulować tzw. “szerokość okna”

plot(density(daneSoc$wiek, bw=1.5), main="Rozklad wieku")

ggplot(daneSoc, aes(x = wiek)) + geom_density() + labs(x = "Wzrost", y = "gęstość")

ggplot(daneSoc, aes(x = wiek, fill = plec)) + geom_density() + labs(x = "Wzrost", y = "gęstość", fill = "Płeć")

ggplot(daneSoc, aes(x = wiek, fill = plec)) + geom_density() + labs(x = "Wzrost", y = "gęstość", fill = "Płeć")+facet_grid(plec ~ .)

Zadanie Narysuj rozkład Couchego

Dystrybuanta

Dystrybuante generujemy za pomocą funkcji ecdf().

plot(ecdf(x))

plot(ecdf(daneSoc$wiek), main = "Dystrybuanta wieku")

ZADANIE

Narysuj dystrybuanty wieku dla kobiet i mężczyzn (zdefiniowane w kolumnie plec)

Data presentation ggplot2

Load library

# load library ggplot 
library(ggplot2)

Nasze dane dostępne bezpośrednio w R

# our data
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

# info o typach danych
str(mtcars)

## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

Including Plots

mtcars %>% 
    ggplot(aes(x=cyl, y=mpg))

mtcars %>% 
    ggplot(aes(x=cyl, y=mpg))+geom_point()

Ponieważ liczba cylindrów (cyl) jest zmienną kategoryjną oznaczmy ja tak w ggplot i zobaczmy jak zmieni to wykres:

mtcars %>% 
    ggplot(aes(x=factor(cyl), y=mpg))+geom_point()

Zobaczmy jak wygląda zależność ilości mil na galon (mpg) w zależności od masy pojazdu (weight - in thousands of pounds)

ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point()

Ciekawym sposobem wizualizacji kolejnych wymiarów (3-d) jest dodanie koloru. Np zobaczmy jak powyższy wykres przedstawia się wraz ze zmienną opisującą pojemność skokową silnika.

ggplot(mtcars, aes(x = wt, y = mpg, color=disp)) +
  geom_point()

oprócz koloru można zmienić również rozmiar punktów w zależności od kolejnej zmiennej co daje nam możliwość przedstawienia kolejnego wymiaru (4-d)

ggplot(mtcars, aes(x = wt, y = mpg, size=disp)) +
  geom_point()

W przypadku zmiennej kategoryjnej można również użyć shape aby odróżnić kategorie kształtem.

GGPlot Warstwy

1. First ggplot2 Layers - DATA

Iris (R.A. Fischer) 3 gatunku Irysa - setosa, versicolor 5 zmiennych

2. Second Layers - Aesthatics co z X i Y

3. Third Layers - Geometries

wystarczy aby zrobić plot

4. Facets

5. Statistics Modelowanie

6. Coordinates

7. Themes Wygląd

Exploring ggplot

Dane - 50 tyś diamentów

zawiera informacje o karacie - carat oraz cenie prize

head(diamonds)

## # A tibble: 6 x 10
##   carat cut       color clarity depth table price     x     y     z
##   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23  Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43
## 2 0.21  Premium   E     SI1      59.8    61   326  3.89  3.84  2.31
## 3 0.23  Good      E     VS1      56.9    65   327  4.05  4.07  2.31
## 4 0.290 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63
## 5 0.31  Good      J     SI2      63.3    58   335  4.34  4.35  2.75
## 6 0.24  Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48

str(diamonds)

## Classes 'tbl_df', 'tbl' and 'data.frame':    53940 obs. of  10 variables:
##  $ carat  : num  0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
##  $ cut    : Ord.factor w/ 5 levels "Fair"<"Good"<..: 5 4 2 4 2 3 3 3 1 3 ...
##  $ color  : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 2 2 6 7 7 6 5 2 5 ...
##  $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 2 3 5 4 2 6 7 3 4 5 ...
##  $ depth  : num  61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
##  $ table  : num  55 61 65 58 58 57 57 55 61 61 ...
##  $ price  : int  326 326 327 334 335 336 336 337 337 338 ...
##  $ x      : num  3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
##  $ y      : num  3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
##  $ z      : num  2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...

ggplot(diamonds, aes(x = carat, y = price))

ggplot(diamonds, aes(x = carat, y = price))+ geom_point()

ggplot(diamonds, aes(x = carat, y = price)) + geom_point() + geom_smooth()

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(diamonds, aes(x = carat, y = price)) + geom_smooth()

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(diamonds, aes(x = carat, y = price, color=clarity)) + geom_smooth()

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(diamonds, aes(x = carat, y = price, color=clarity)) + geom_point(alpha=0.4)

dia_plot <- ggplot(diamonds, aes(x = carat, y = price))
dia_plot <- dia_plot + geom_point(alpha=0.2)
dia_plot + geom_smooth(se = FALSE) # bez przedziałów ufności

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'