3  Rozkłady teoretyczne i ich nieparametryczna estymacja - wybrane biblioteki programu R

“All models are wrong, but some are useful.”

“Wszystkie modele są błędne, ale niektóre są użyteczne.”

George E. P. Box (Box, 1979, s. 202)

Współczesna analiza danych statystycznych wykracza poza same obliczenia numeryczne, kładąc coraz większy nacisk na zaawansowaną wizualizację, która jest kluczowa dla zrozumienia natury badanych zjawisk i weryfikacji założeń modelowych. Graficzna prezentacja wyników wnioskowania pozwala nie tylko na skuteczniejszą interpretację skomplikowanych zależności, ale także ułatwia komunikację wniosków odbiorcom o różnym poziomie zaawansowania. Choć w środowisku R standardem stał się pakiet ggplot2, specyficzne zagadnienia statystyczne – takie jak wizualizacja obszarów krytycznych, operacje na zmiennych losowych czy estymacja jądrowa – często wymagają zastosowania wyspecjalizowanych rozszerzeń i dodatkowych bibliotek.

Niniejszy rozdział koncentruje się na przeglądzie i praktycznym zastosowaniu wybranych pakietów środowiska R, które wspomagają prezentację rozkładów teoretycznych oraz ich oszacowań. Omówione zostaną narzędzia do automatyzacji wykreślania funkcji gęstości i dystrybuant (m.in. ggfortify, ggnormalviolin), pakiety umożliwiające pracę na mieszankach i funkcjach zmiennych losowych (convdistr), a także metody nieparametrycznej estymacji gęstości jedno- i wielowymiarowej z wykorzystaniem bibliotek takich jak KernSmooth, MASS czy lattice. Celem rozdziału jest zaprezentowanie zestawu narzędzi pozwalających na kompleksową wizualizację probabilistycznych aspektów analizy danych.

3.1 Rozkłady teoretyczne i ich estymacja - wybrane biblioteki programu R

Kończak (2024) wskazuje na różnorodne biblioteki pozwalające rozszerzyć możliwości pakietu graficznego ggplot2. W tej części skupiono się na pakietach, które mogą być zastosowane w graficznych prezentacjach rezultatów wnioskowania statystycznego. Wykaz pakietów, które zostały wykorzystane w tym rozdziale przedstawiono w tabeli tabela 4.1.

Tabela 3.1. Wybrane biblioteki programu R wspomagające prezentację rozkładów teoretycznych, ich funkcji oraz nieparametrycznej estymacji gęstości.
Biblioteka Opis
ggfortify Automatyzacja wykreślania wybranych wykresów rozkładów teoretycznych.
ggnormalviolin Konstrukcja normalnych wykresów wiolinowych z określonymi średnimi i odchyleniami standardowymi.
convdistr Łączenie rozkładów prawdopodobieństwa przy użyciu funkcji generatora liczb losowych każdego rozkładu.
ggplot2 Ogólne środowisko graficzne oparte na Grammar of Graphics; tworzenie elastycznych wykresów rozkładów (histogramy, gęstości, QQ‑ploty, wykresy rozrzutu z gęstościami brzegowymi itp.).
lattice System grafiki trellis do wielowymiarowych wizualizacji, w tym prezentacji rozkładów warunkowych i porównywania gęstości w wielu grupach.
MASS Zbiór funkcji i danych do klasycznej książki Venablesa i Ripleya; generowanie zmiennych z różnych rozkładów, dopasowanie modeli i rozkładów, m.in. funkcje pomocne przy analizie rozkładów empirycznych.
KernSmooth Implementacja jądrowej estymacji gęstości i regresji; umożliwia wygładzoną, nieparametryczną estymację gęstości oraz jej wizualizację.
see Pakiet do atrakcyjnej wizualizacji wyników modeli statystycznych, rozkładów a posteriori oraz przedziałów niepewności.
data.table Wysokowydajne narzędzie do przetwarzania i agregacji danych, ułatwiające przygotowanie dużych zbiorów danych do wizualizacji rozkładów i estymacji gęstości.
ggmosaic Rozszerzenie ggplot2 do tworzenia wykresów mozaikowych (mosaic plots) prezentujących rozkłady wielowymiarowych zmiennych jakościowych.
tidyverse Zestaw spójnych pakietów (m.in. dplyr, tidyr, readr) do pracy z „uporządkowanymi” danymi, ułatwiający przygotowanie danych do wizualizacji w ggplot2 oraz do estymacji gęstości.

Na wstępie konieczne jest załadowanie bibliotek, które będą wykorzystywane w tym rozdziale. Z zakresu wnioskowania statystycznego są to biblioteki:

library(ggfortify)
library(ggnormalviolin)
library(convdistr)

library(ggplot2)
library(lattice)
library(ggpubr)

library(MASS)
library(KernSmooth)
library(see)
library(data.table)
library(ggmosaic)
library(tidyverse)
library(gridExtra)

3.2 Rozkłady teoretyczne zmiennych losowych - gęstość i dystrybuanta

Pakiet ggfortify to narzędzie pozwalające na konstrukcję wykresów z wynikami dla różnych analiz statystycznych. Pakiet umożliwia między innymi przedstawienie wykresów dla modeli liniowych, metod klasyfikacji oraz szeregów czasowych z uwzględnieniem prognoz. W tej części uwaga jednak zostanie skoncentrowana wyłącznie na wykreślaniu funkcji gęstości i dystrybuanty rozkładów teoretycznych wybranych zmiennych losowych. Rozkłady teoretyczne stanowią podstawę konstrukcji przedziałów ufności jak również parametrycznych testów statystycznych. Do najczęściej wykorzystywanych rozkładów teoretycznych we wnioskowaniu statystycznym należy zaliczyć m.in. rozkład normalny, rozkład \(t\) Studenta oraz rozkład chi-kwadrat.

Wykres gęstości rozkładu normalnego standardowego uzyskuje się za pomocą następującego kodu, którego wynik realizacji przedstawiono na rys. 3.1.

ggdistribution(dnorm, seq(-4, 4, 0.1), mean = 0, sd = 1,linetype=1,fill = 'blue',alpha=0.1) + theme_minimal()
Rysunek 3.1. Gęstość rozkładu normalnego standardowego

Często prezentując rezultat testowania hipotez statystycznych przydatna jest prezentacja wartości statystyki testowej na tle odpowiedniego obszaru krytycznego. Pakiet ggfortify umożliwia dla określonego rozkładu i przyjętego poziomu istotności \(\alpha\) zaznaczenie na wykresie obszaru krytycznego. Taką możliwość w przypadku rozkładu normalnego standardowego przedstawia poniższy kod.

alpha=0.05
q=qnorm(1-alpha/2)
r1=ggdistribution(dnorm, seq(-q, q, 0.1), mean = 0, sd = 1,fill = 'yellow',alpha=0.3)
r2=ggdistribution(dnorm, seq(-4, 4, 0.1), mean = 0, sd = 1,linetype=1,fill = 'blue',alpha=0.1,p=r1)
r2+xlim(-4,4)+ylim(0,0.4)+labs(title="Obszar krytyczny",  x="x", y="f(x)")+theme_minimal()
Rysunek 3.2. Gęstość rozkładu normalnego standardowego z wyróżnionym obszarem krytycznym

W wyniku realizacji powyższego kodu otrzymuje się rys. 3.2. Wykres przedstawia dwustronny obszar krytyczny dla poziomu istotności \(\alpha=0,05\) dla zmiennej losowej o rozkładzie normalnym standardowym.

Pakiet ggfortify umożliwia przedstawienie na jednym wykresie kilku funkcji gęstości. Poniższy kod prezentuje umieszczenie na wykresie gęstości rozkładu normalnego standardowego oraz dwóch funkcji gęstości rozkładu t-Studenta dla stopni swobody 3 (linia czerwona) i 30 (linia niebieska).

r1=ggdistribution(dnorm, seq(-q, q, 0.1), mean = 0, sd = 1,fill = 'blue',alpha=0.1)
r2=ggdistribution(dnorm, seq(-4, 4, 0.1), mean = 0, sd = 1,linetype=1,fill = 'yellow',alpha=0.1,p=r1)
r3=ggdistribution(dt, seq(-4, 4, 0.1), df=3,colour='red',p=r2)
ggdistribution(dt, seq(-5, 5, 0.1), df=30,colour='blue',p=r3)+theme_minimal()
Rysunek 3.3. Gęstości rozkładu normalnego standardowego i rozkładów t-Studenta

Na rysunku rys. 3.3 przedstawiono gęstość rozkładu normalnego oraz gęstości 2 zmiennych losowych o rozkładzie \(t\) Studenta z liczbą stopni swobody 3 oraz 30. Zauważalne jest, że gęstość rozkładu \(t\) Studenta dla \(df\) = 30 stopni swobody prawie pokrywa się z gęstością rozkładu normalnego standardowego.

Dla wykreślenia dystrybuanty zmiennej losowej o rozkładzie normalnym standardowym należy wykonać następujący kod:

ggdistribution(pnorm, seq(-3, 3, 0.1), mean = 0, sd = 1, colour = 'red',xlab='x',ylab='F(x)')+theme_minimal()
Rysunek 3.4. Dystrybuanta rozkładu normalnego standardowego

Na rysunku rys. 3.4 przedstawiono dystrybuantę rozkładu normalnego standardowego.

Poniższy kod pozwala na jednym wykresie zamieścić funkcje gęstości trzech zmiennych losowych o rozkładzie chi-kwadrat odpowiednio z liczbą stopni swobody \(df = 3, \ 5 \text{ oraz } 10\). Wynik realizacji kodu został przedstawiony na rys. 3.5.

r1=ggdistribution(dchisq, seq(0, 15, 0.1), df=3,colour='red',fill='red',alpha=0.1)
r2=ggdistribution(dchisq, seq(0, 15, 0.1), df=5,colour='blue',p=r1,fill='blue',alpha=0.1)
ggdistribution(dchisq, seq(0, 15, 0.1), df=10,colour='green',p=r2,fill='green',alpha=0.1)+theme_minimal()
Rysunek 3.5. Gęstości rozkładów chi-kwadrat

3.3 Rozkłady teoretyczne zmiennych losowych - wykresy wiolinowe

Nieco inne podejście do prezentacji graficznej teoretycznych rozkładów prawdopodobieństwa umożliwia pakiet ggnormalviolin. Do prezentacji graficznej wykorzystywane są wykresy wiolinowe pozwalające w sposób czytelny porównać charakterystyki kilku rozkładów. Taki przykład z wykresami gęstości trzech zmiennych losowych o rozkładzie normalnym z różnymi parametrami prezentuje poniższy kod.

d <- data.frame(
  zmienna = c("X", "Y", "Z"),
  zmienna_mean = c(80, 110, 130),
  zmienna_sd = c(15,  20, 5))
p <- ggplot(data = d, 
            aes(x = zmienna,
                mu = zmienna_mean,
                sigma = zmienna_sd,
                fill = zmienna)) +
  theme(legend.position = "none")+
  coord_flip()
p + geom_normalviolin()+theme_minimal()
Rysunek 3.6. Gęstość trzech zmiennych o rozkładach normalnych z wartościami oczekiwanymi: 80, 110 i 130 oraz odchyleniami standardowymi: 15, 20 i 5

W efekcie realizacji tego kodu uzyskuje się wykres jak na rys. 3.6. Na tym wykresie za pomocą graficznej reprezentacji w postaci wiolin przedstawiono gęstości trzech zmiennych losowych o rozkładach normalnych z wartościami oczekiwanymi kolejno \(\mu\) = 80, 110 i 130 oraz odpowiednio odchyleniami standardowymi \(\sigma\) = 15, 20 i 5.

Podobnie jak w przypadku wcześniej omawianego pakietu ggfortify tak i w tym przypadku jest możliwość zaznaczenia obszarów krytycznych. Przykład taki dla wcześniej ujętych zmiennych losowych przy przyjęciu poziomu istotności \(\alpha\) = 0,05 realizuje poniższy kod, a wynik jego realizacji przedstawiono na rys. 3.7.

p + geom_normalviolin(p_tail = 0.05)+theme_minimal()
Rysunek 3.7. Gęstość czterech zmiennych o rozkładach normalnych z wartościami oczekiwanymi: 80, 110 i 130 oraz odchyleniami standardowymi: 15, 20 i 5 z zaznaczonymi obszarami dwustronnymi krytycznymi dla \(\alpha\) = 0,05

3.4 Wizualizacja weryfikacji normalności rozkładu zmiennej

Poniższe rozważania zostaną przeprowadzone z wykorzystaniem danych ze zbioru iris (Fisher, 1936; Kończak, 2024). Zbiór danych iris to jeden z najbardziej znanych i najczęściej używanych zbiorów danych w dziedzinie statystyki. Został on wprowadzony przez brytyjskiego statystyka i biologa Ronalda Aylmera Fishera w jego pracy z 1936 roku (Fisher, 1936), jako przykład w zastosowaniu analizy dyskryminacyjnej.

Zbiór składa się ze 150 obserwacji i 5 następujących zmiennych:

  • Sepal.Length: Długość kielicha (działki kielicha) w centymetrach.
  • Sepal.Width: Szerokość kielicha (działki kielicha) w centymetrach.
  • Petal.Length: Długość płatka w centymetrach.
  • Petal.Width: Szerokość płatka w centymetrach.
  • Species: Gatunek irysa (zmienna jakościowa o trzech wariantach: setosa, versicolor i virginica).

W pierwszym kroku zostanie sprawdzona normalność rozkładu długości płatków kielicha wraz z graficzną wizualizacją w formie wykresu \(QQ\)-plot.

ggplot(iris, aes(sample = Petal.Length)) +
  stat_qq() + 
  stat_qq_line(color = "red", linetype = "dashed") + 
  labs(x = "Kwantyle teoretycznego rozkładu normalnego",
    y = "Kwantyle empiryczne") +
  theme_minimal() + 
  coord_fixed(ratio = .4) 
Rysunek 3.8. Rozkład długości płatka (Petal.Length) a rozkład normalny

Na rys. 3.8 widoczne jest, że analizowana zmienna ma inny rozkład niż normalny. W analizowanym przypadku wyróżnione są trzy gatunki kwiatu iris. W takiej sytuacji niezasadne jest sprawdzanie normalności rozkładu charakterystyk dla gatunków łącznie, ale dla każdego z gatunków osobno. Zrealizowano to w formie graficznej następująco:

ggplot(iris, aes(sample = Petal.Length)) +
  stat_qq() + 
  stat_qq_line(color = "red", linetype = "dashed") + 
  facet_wrap(~ Species, ncol = 3) + 
  labs( x = "Kwantyle teoretycznego rozkładu normalnego", y = "Kwantyle empiryczne") +
  theme_minimal() + 
  coord_fixed(ratio = 1) 
Rysunek 3.9. Rozkład długości płatka (Petal.Length) a rozkład normalny według gatunku

Na rys. 3.9 przedstawiono trzy wykresy \(QQ\)-plot, dla zmiennej Petal.Length dla każdego z gatunków z osobna. Widoczne jest, że we wszystkich tych przypadkach punkty na wykresie są położone w pobliżu linii teoretycznej. Dla uzyskania formalnego potwierdzenia wskazane jest przeprowadzenie testu normalności

shapiro.test(iris$Petal.Length[iris$Species=='setosa'])$p.value
# [1] 0.05481147
shapiro.test(iris$Petal.Length[iris$Species=='versicolor'])$p.value
# [1] 0.1584778
shapiro.test(iris$Petal.Length[iris$Species=='virginica'])$p.value
# [1] 0.1097754

Przeprowadzenie testu normalności we wszystkich przypadkach, przy przyjęciu poziomu istotności \(\alpha=0,05\), prowadzi do stwierdzenia braku podstaw do odrzucenia hipotezy głoszącej normalność rozkładu zmiennej Petal.Length.

Dla testowania normalności i odpowiedniej wizualizacji niekiedy wskazane jest przedstawienie wyników na kilku wykresach. Poniższy kod określa funkcję diag_nor() i prezentuje wyniki na wykresie.

diag_nor <- function(dane) {
  shapiro_test <- shapiro.test(dane)
  p1 <- ggplot(data.frame(x = dane), aes(x = x)) +
    geom_histogram(aes(y = after_stat(density)), 
                   bins = 10,           
                   color = "black",     
                   fill = "skyblue", 
                   alpha = 0.7) +
    stat_function(fun = dnorm, 
                  args = list(mean = mean(dane), sd = sd(dane)),
                  color = "red", size = 1) +
    theme_bw()
  p2 <- ggqqplot(dane) +
    theme_bw()
  
  wyniki <- data.frame(
    Test = "Shapiro-Wilk",
    Statystyka = round(shapiro_test$statistic, 4),
    p_wartość = round(shapiro_test$p.value, 4),
    Wniosek = ifelse(shapiro_test$p.value > 0.05, 
                     "Rozkład normalny", 
                     "Brak normalności")
  )
  p3 <- tableGrob(wyniki)
  grid.arrange(p1, p2, p3, ncol = 2,
               layout_matrix = rbind(c(1, 2), c(3, 3)))
}

dane <- iris %>%
  filter(Species == "setosa") %>%
  select(Petal.Length) %>%
  unlist(use.names = FALSE) 
diag_nor(dane)
Rysunek 3.10. Diagnoza normalności rozkładu długości płatka (Petal.Length)

Funkcja diag_nor() przeprowadza kompleksową diagnostykę normalności rozkładu zmiennej numerycznej, tworząc panel trzech elementów wizualizacyjnych. Generuje histogram z nałożoną teoretyczną krzywą rozkładu normalnego (górny lewy), wykres kwantyl-kwantyl (\(QQ\) plot) porównujący rozkład empiryczny z teoretycznym (górny prawy) oraz tabelę z wynikami testu Shapiro-Wilka zawierającą statystykę testową, \(p\)-wartość i wniosek o normalności rozkładu (dolny pas). W przykładzie użycia funkcja analizuje długość płatków (Petal.Length) dla gatunku setosa z zestawu danych iris (por. rys. 3.10).

3.5 Nieparametryczna estymacja gęstości

Jednym z ważniejszych zagadnień nieparametrycznego wnioskowania statystycznego jest estymacja funkcji gęstości. Zagadnienie jest szeroko opisywane w literaturze (Baszczyńska, 2016; Domański & Pruska, 2000).

W pierwszej kolejności estymacja nieparametryczna długości płatka kwiatu iris zostanie wykonana z pakietem KernSmooth.

x <- iris$Petal.Length
plot(bkde(x), type = "l", xlab='x',ylab='kde(x)')
Rysunek 3.11. Estymacja gęstości z pakietem KernSmooth. Gęstość długości płatka (Petal.Length)

Na rys. 3.11 przedstawiono rezultat powyższego kodu. Wykres przedstawia estymację gęstości zmiennej Petal.Length.

Podobne możliwości jak pakiet KernSmooth w zakresie graficznej prezentacji wyników estymacji nieparametrycznej zapewnia pakiet ggplot2. Dla tych samych danych graficzna prezentacja zrealizowana w oparciu o ten pakiet jest następująca:

ggplot(iris, aes(x = Petal.Length)) +
  geom_density(fill = "lightblue", alpha = 0.7) + # Dodajmy trochę estetyki
  labs( x = "Długość płatka (cm)", y = "Gęstość"  ) +
  theme_minimal()
Rysunek 3.12. Estymacja długości płatka (Petal.Length)

Na rys. 3.12 przedstawiono rezultat powyższego kodu. Wykres przedstawia estymację gęstości zmiennej Petal.Length. W analizowanym przypadku wskazane jest przeprowadzenie estymacji gęstości zmiennej dla poszczególnych gatunków z osobna. Takie zadanie realizuje poniższy kod:

ggplot(iris, aes(x = Petal.Length)) +
  geom_density(fill = "lightblue", alpha = 0.7) + 
  facet_wrap(~ Species, ncol=1,scales = "free_y") +
  labs(x = "Długość płatka (cm)", y = "Gęstość" ) +
  theme_minimal() 
Rysunek 3.13. Estymacja gęstości długości płatka (Petal.Length) według gatunku

Na rys. 3.13 przedstawiono oszacowania gęstości zmiennej Petal.Length według gatunku.

Uzyskanie w bardzo wygodny sposób jedno- i dwuwymiarowych ocen gęstości i ich wizualizację umożliwia także pakiet lattice. Dla wszystkich czterech zmiennych liczbowych ze zbioru iris według gatunku realizuje to poniższy kod:

densityplot(~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width | Species,
            data = iris,
            auto.key = list(columns = 3), 
            plot.points = FALSE, 
            col = c("blue", "red", "darkgreen"), 
            xlab = "Zmienna") 
Rysunek 3.14. Estymacja gęstości zmiennych numerycznych zbioru iris według gatunku

Na rys. 3.14 przedstawiono estymację gęstości zmiennych Petal.Length, Petal.Width, Sepal.Length i Sepal.Width ze zbioru iris według gatunku.

Nieco inną formę prezentacji tych samych oszacowań przedstawia poniższy kod, którego rezultat prezentuje rys. 3.15.

iris_long <- data.frame(
  value = unlist(iris[, 1:4]),
  variable = rep(names(iris[, 1:4]), each = nrow(iris)),
  Species = rep(iris$Species, times = ncol(iris[, 1:4])))
densityplot(~ value | variable * Species,
            data = iris_long,
            auto.key = list(columns = 3),
            plot.points = FALSE,
            col = c("blue", "red", "darkgreen"),
            xlab = "Wartość zmiennej",
            scales = list(x = list(relation = "free")))
Rysunek 3.15. Estymacja gęstości zmiennych numerycznych zbioru iris według zmiennej i gatunku

Pakiet lattice umożliwia także uzyskanie ocen dla dwuwymiarowych zmiennych. Przykład takiego rozwiązania przedstawia poniższy kod:

get_density_data <- function(data, x_var, y_var, species_name) {
  x <- data[[x_var]]
  y <- data[[y_var]]
  if (length(x) < 2 || length(y) < 2) {
    return(NULL) 
  }
  dens <- kde2d(x, y, n = 50) 
  data.frame(
    x = rep(dens$x, each = length(dens$y)),
    y = rep(dens$y, times = length(dens$x)),
    z = as.vector(dens$z),
    Species = species_name)
}
density_list <- lapply(levels(iris$Species), function(s) {
  subset_iris <- subset(iris, Species == s)
  get_density_data(subset_iris, "Sepal.Length", "Sepal.Width", s)
})
density_df <- do.call(rbind, density_list)
contourplot(z ~ x * y | Species,
            data = density_df,
            cuts = 10,
            labels = FALSE,
            col.regions = topo.colors,
            xlab = "Sepal.Length",
            ylab = "Sepal.Width",
            panel = function(x, y, z, subscripts, ..., B = NULL, A = NULL) {
              current_species_name <- levels(density_df$Species)[panel.number()]
              panel.contourplot(x, y, z, subscripts, ...)
              points_data_filtered <- subset(iris, Species == current_species_name)
              panel.points(points_data_filtered$Sepal.Length, points_data_filtered$Sepal.Width,
                           col = "grey", pch = 16, cex = 0.6)
            })
Rysunek 3.16. Estymacja dwuwymiarowej gęstości (Sepal.Length i Sepal.Width)

Na rys. 3.16 przedstawiono dwuwymiarowe oceny gęstości zmiennych Sepal.Length i Sepal. Width według gatunków. Nieco inną formę graficzną tych samych ocen prezentuje kolejny kod, którego wynik przedstawia rys. 3.17.

levelplot(z ~ x * y | Species,
          data = density_df,
          col.regions = heat.colors(100), 
          xlab = "Sepal.Length",
          ylab = "Sepal.Width")
Rysunek 3.17. Mapa ciepła - estymacja dwuwymiarowej gęstości (Sepal.Length i Sepal.Width)

Wykorzystanie pakietu KernSmooth dla estymacji gęstości dwuwymiarowych przedstawia poniższy kod:

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
  geom_density_2d() + 
  geom_point(alpha = 0.3, size = 1) + 
  facet_wrap(~ Species) + 
  labs(x = "Długość kielicha (cm)", y = "Szerokość kielicha (cm)" ) +
  theme_minimal()
Rysunek 3.18. Estymacja gęstości dwuwymiarowej (długość i szerokość kielicha) według gatunku

Na rys. 3.18 przedstawiono wizualizację estymacji gęstości rozkładu dwuwymiarowego, dla zmiennych długość (Sepal.Length) i szerokość kielicha (Sepal.Width) według gatunku.

3.6 Funkcje i mieszanki zmiennych losowych

W badaniach naukowych bardzo często należy się odwołać nie do zmiennej losowej, której rozkład jest znany, ale do kombinacji kilku zmiennych losowych lub do ich mieszanek. Często prowadzi to do bardzo złożonych rozważań teoretycznych, a w wielu przypadkach uzyskanie postaci rozkładu zmiennej losowej na drodze rozważań teoretycznych jest wręcz niemożliwe. W takim przypadku pomocny może być pakiet convdistr, który wykorzystując metody symulacyjne pozwala na przedstawienie graficznego obrazu kombinacji lub mieszanek zmiennych losowych. Pakiet ten umożliwia również generowanie obserwacji losowych z takiego rozkładu, co ma duże znaczenie w analizach symulacyjnych. Dla uzyskania przybliżonego graficznego obrazu zmiennej losowej o rozkładzie normalnym standardowym można wykorzystać następujące komendy.

x <- new_NORMAL(0,1)
ggDISTRIBUTION(x)+theme_minimal()
Rysunek 3.19. Estymacja gęstości rozkładu normalnego standardowego

Na rysunku rys. 3.19 przedstawiono wizualizację estymacji gęstości zmiennej losowej o rozkładzie normalnym standardowym.

Pakiet convdistr umożliwia bardzo wygodną ocenę rozkładów funkcji zmiennych losowych. Rozkład sumy dwóch zmiennych losowych o rozkładach jednostajnych na przedziale (0, 1) jest rozkładem trójkątnym określonym na przedziale (0, 2). Symulacyjną estymację takiego rozkładu przedstawia rys. 3.20, który został uzyskany po wykonaniu następującego kodu:

d1 <- new_UNIFORM(0,1)
d2 <- new_UNIFORM(0,1)
dsum <- new_SUM(list(d1,d2))
ggDISTRIBUTION(dsum)+theme_minimal()
Rysunek 3.20. Estymacja gęstości sumy dwóch zmiennych losowych o rozkładzie jednostajnym

Dla uzyskania rozkładu sumy trzech zmiennych losowych o rozkładach normalnym, jednostajnym i Poissona należy wykonać poniższy kod, którego rezultat przedstawia rys. 3.21.

d1 <- new_NORMAL(1,1)
d2 <- new_UNIFORM(2,8)
d3 <- new_POISSON(5)
dsum <- new_SUM(list(d1,d2,d3))
ggDISTRIBUTION(dsum)+theme_minimal()
Rysunek 3.21. Estymacja gęstości sumy zmiennych losowych o rozkładach normalnym, jednostajnym i Poissona

Dla uzyskania graficznego obrazu zmiennej losowej \(Y\) będącej sumą zmiennych losowych o rozkładzie normalnym (\(X_1\)) oraz zmiennej losowej będącej iloczynem zmiennych losowych o rozkładach Poissona (\(X_2\)) oraz beta (\(X_3\)) odpowiednio z parametrami jak poniżej:

\(Y = X1 + X2*X3\)

gdzie \(X_1 \sim N(1, \ 0,5)\), \(X_2 \sim Poisson(5)\), \(X_3 \sim Beta(10,20)\)

należy wykonać następujący kod:

X1 <- new_NORMAL(1,0.5)
X2 <- new_POISSON(5)
X3 <- new_BETA(10,20)
res <- X1 + X2 * X3
ggDISTRIBUTION(res) + ggtitle("X1 + X2 * X3")+theme_minimal()
Rysunek 3.22. Estymacja gęstości sumy zmiennych losowych o rozkładach normalnym oraz iloczynu zmiennych losowych o rozkładach Poissona i beta

W wyniku realizacji powyższego kodu otrzymuje się graficzną prezentację estymacji gęstości zmiennej losowej \(Y\) (zob. rys. 3.22).

Pakiet convdistr umożliwia także uzyskanie wizualizacji mieszanek rozkładów zmiennych losowych. Dla uzyskania mieszanki trzech rozkładów zmiennych losowych o rozkładach normalnych o różnych wartościach oczekiwanych i wariancjach należy wykonać następujący kod

d1 <- new_NORMAL(2,2)
d2 <- new_NORMAL(8,1)
d3 <- new_NORMAL(12,3)
dmix <- new_MIXTURE(list(d1,d2,d3))
ggDISTRIBUTION(dmix)+ggtitle("")+theme_minimal()
Rysunek 3.23. Estymacja gęstości mieszanki trzech zmiennych losowych o rozkładach normalnych

rys. 3.23 przedstawia wynik realizacji powyższego kodu. Jest na nim zaprezentowana gęstość mieszanki trzech zmiennych losowych o rozkładach normalnych z różnymi parametrami.

Ogromną zaletą rozważanego pakietu jest możliwość uzyskania losowych wartości z otrzymanych rozkładów zmiennych losowych - funkcji lub mieszanek. Uzyskanie 1000 losowych wartości dla ostatniej mieszanki i przedstawienie wykresu jest realizowane następująco:

dane <- rfunc(dmix, 1000)
ggplot(dane,aes(rvar))+geom_histogram()+theme_minimal()
Rysunek 3.24. Histogram 1000 losowych wartości z rozkładu mieszanki zmiennych d1, d2 i d3
Baszczyńska, A. (2016). Parametr wygładzania w estymacji jądrowej funkcji gęstości dla zmiennych losowych w badaniach ekonomicznych (I). Wydawnictwo Uniwersytetu Łódzkiego.
Box, G. E. P. (1979). Robustness in the strategy of scientific model building. W R. L. Launer & G. N. Wilkinson (Red.), Robustness in Statistics (s. 201–236). Academic Press.
Domański, C., & Pruska, K. (2000). Nieklasyczne metody statystyczne (s. 368). Polskie Wydawnictwo Ekonomiczne.
Fisher, R. A. (1936). The Use of Multiple Measurements in Taxonomic Problems. Annals of Eugenics, 7(2), 179–188.
Kończak, G. (2024). Wizualizacja wyników badań naukowych. Zasady, metody i narzędzia. Uniwersytet Ekonomiczny w Katowicach.