#Rozkład normalny
# 4 podstawowe funkcje przydane w opisie rozkładu
# rnorm - generuje losowe liczby z rozkładu normalnego

# 1) losuję 2 próby: 1000 wartości ze standardowego rozkładu normalnego (średnia = 0, odch.stand =1)
# oraz 500 wartości z r.n. o średniej 3 i odch.stand.0.7

losnorm <-rnorm(1000,0,1) #1 liczba - ile wartości? 2liczba - jaka średnia? 3liczba - jakie odch.stand?
losnorm1 <- rnorm(500, 3, 0.7)

#Sprawdźmy jakie w rzeczywistości wylosowaliśmy próby:
#Przypomnienie: policz średnią (mean())i odchylenie standardowe (sd()) dla obu prób
mean(losnorm)
mean(losnorm1)

#Porównajmy te próby na wykresie pudełkowym
boxplot(losnorm,losnorm1, main="Losowe punkty", 
        names =c("Standard","Mean=3"), ylab="Wartość x", varwidth=TRUE)

#oraz na wykresie skrzypcowym
#install.packages("vioplot")
library(vioplot)
vioplot(losnorm, losnorm1, names=c("Standard", "Mean=3"), 
        col=c("gold","green"))
title("Wykres skrzypcowy")

###################################################################################
#Ćwiczenie 1. Wylosuj 2 próby: 100 wartości z rozkładu normalnego
#o średniej 4 i odchyleniu standardowym 1.5 oraz 100 wartości
# z r.n. o średniej 4 i odch. stand 2.8. Porównaj te próby na 
#wykresie skrzypcowym



###################################################################################
#Ćwiczenie 2. Wylosuj 2 próby: 100 wartości z rozkładu normalnego
#o średniej 1 i odchyleniu standardowym 1.5 oraz 100 wartości
# z r.n. o średniej 4 i odch. stand 1.5. Porównaj te próby na 
#wykresie skrzypcowym


#2) porównanie r.norm. z innymi rozkładami - boxplot vs. violin plot
par(mfrow=c(1,2)) # dzielenie obszaru wykresu na dwie kolumny
mu<-2 #średnia
si<-0.6 #odchylenie standardowe
bimodal<-c(rnorm(1000,-mu,si),rnorm(1000,mu,si)) #losowanie wartości z rozkładu bimodalny
uniform<-runif(2000,-4,4) #losowanie wartości z rozkładu unimodalnego
normal<-rnorm(2000,0,3) #losowanie wartości z rokładu...........<- tu wpisz właściwą nazwę
boxplot(bimodal,uniform,normal) #porównanie rozkł. z użyciem boxplot
vioplot(bimodal,uniform,normal) #porównanie rozkł z użyciem wykresu skrzypcowego

#histogram jako narzędzie do oceny rozkładów
par(mfrow=c(1,1)) #wracam do pojedynczej kolumny
hist(normal, main="Losowe wartości z rozkładu normalnego standardowego", ylab="Częstość",
     cex.axis=.8, xlim=c(-10,10))

########################################################################################
#Ćwiczenie 3. Utwórz histogramy dla rozkładu dwumodalnego oraz unimodalnego





#przykład - chcę wybrać 100 liczb losowych z rozkładu normalnego 
#o średniej 2 i odchyleniu standardowym 0.5 np. jako kontrola do pewnych pomiarów
par(mfrow=c(1,2))
x <- rnorm(100,2,.5)
mean(x)
par(mfrow = c(1,1))
# tworzymy histogram dla wygenerowanych punktów
hist(x, main="Losowe wartości", ylab="Częstość",
     cex.axis=.8, xlim=c(-2,6))

# tworzymy histogram w którym mamy wiecej słupków
hist(x, breaks=10, main="Losowe wartości", ylab="Częstość",
     cex.axis=.8, xlim=c(-2,6))

# tworzymy histogram z przerwami dokładnie tam gdzie chcemy
brks <- seq(0,4,0.1) # najpierw tworzymy wektor z przerwami co 0.1
brks
# teraz możemy utworzyć histogram z odstępami co 0.1
hist(x, breaks=brks, main="Losowe wartości", ylab="Częstość",
     cex.axis=.8, xlim=c(-2,6))
par(mfrow = c(1,1))
# tworzymy histogram na którym zamiast częstości widać gęstość prawdopodobieństwa
hist(x, breaks=brks, freq=FALSE, main="Losowe wartości", ylab="Gęstość prawdopodobieństwa",
     cex.axis=.8, xlim=c(-2,6))
curve(dnorm(x, 2, .5), col = 2, lty = 2, lwd = 2, add = TRUE)

##################################################################################
#Ćwiczenie 4. Używając parametru breaks w funkcji hist narysuj histogram o słupku, 
#który odpowiada wartości 1 oraz taki na którym słupek ma wartość 0.2.
# Który histogram Twoim zdaniem lepiej spełnia swoją funkcję?
x
x <- rnorm(10,2,.5)
brks1 <- seq(0,4,1)
brks2 <- seq(0,4,0.2)
par(mfrow=c(1,2))
hist(x, breaks=brks1, freq=FALSE, main="Losowe wartości", ylab="Gęstość prawdopodobieństwa",
     cex.axis=.8, xlim=c(-2,6))
hist(x, breaks=brks2, freq=FALSE, main="Losowe wartości", ylab="Gęstość prawdopodobieństwa",
     cex.axis=.8, xlim=c(-2,6))

# Do opisu rozkładu stosujemy także inne parametry takie jak
# funkcja gęstości prawdopodobieństwa oraz dystrybuanta
xseq <- seq(-5,5,.01) #najpierw utwórzmy sekwencję liczb w przedziale od -5 do +5 w odstępach co 0.01
ges <- dnorm(xseq, 0, 1) #funkcja gęstości prawdopodobieństwa
dys <- pnorm(xseq, 0, 1) #dystrybuanta rozkładu normalnego

# narysujemy teraz wykresy powyższych funkcji
plot(xseq, ges, col="darkgreen",xlab="", 
     ylab="Gęstość", type="l",lwd=2, cex=2, 
     main="Funkcja gęstości prawdopodobieństwa",ylim=c(0,1), cex.axis=0.8) #Gęstość

plot(xseq, dys, col="darkorange", xlab="", 
     ylab="Skumulowane prawdopodobieństwo",type="l",lwd=2, cex=2, 
     main="Dystrybuanta", cex.axis=.8) #Dystrybuanta

################################################################################################
#Ćwiczenie 5. Oblicz wartości funkcji gęstości prawdopodobieństwa dla rozkładu normalnego, 
# który ma średnią średnią 0 i odchylenie standardowe 3 i zwizualizuj ją za pomocą funkcji plot.
# Porównaj z funkcją gęstości dla roz. norm. standardowego.
par(mfrow = c(1,2))
xseq<- seq(-6,6,.01)
ges<- dnorm(xseq, 0, 3)
ges1<- dnorm(xseq, 0, 1)

plot(xseq, ges, col="darkorange", xlab="", ylab= "Gęstość", 
     lwd=2, cex=2, main="Funkcja gęstości prawdopodobieństwa", cex.axis=.8)
plot(xseq, ges1, col="red", xlab="", ylab= "Gęstość", 
     lwd=2, cex=2, main="Funkcja gęstości prawdopodobieństwa", cex.axis=.8)

dys<- pnorm(xseq, 0, 3)
dys1<- pnorm(xseq, 0, 1)

plot(xseq, dys, col="red", xlab="",
     ylab="Skumulowane prawdopodobieństwo",lwd=2, cex=2,
     main="Dystrybuanta", cex.axis=.8)

plot(xseq, dys1, col="darkorange", xlab="",
     ylab="Skumulowane prawdopodobieństwo",lwd=2, cex=2,
     main="Dystrybuanta", cex.axis=.8)

################################################################################################
#Ćwiczenie 6. Porównaj dystrybuanty rozkładu standardowego i rozkładu normalnego z ćw. 5
xseq <- seq(-5,5,.01) #najpierw utwórzmy sekwencję liczb w przedziale od -5 do +5 w odstępach co 0.01
ges <- pnorm(xseq, 0, 1) #
ges1 <- pnorm(xseq, 0, 3) #

par(mfrow = c(1,2))
# narysujemy teraz wykresy powyższych funkcji
plot(xseq, ges, col="darkgreen",xlab="", 
     ylab="Gęstość", type="l",lwd=2, cex=2, 
     main="Dystrybuanta", cex.axis=0.8) #Gęstość

plot(xseq, ges1, col="darkgreen",xlab="", 
     ylab="Gęstość", type="l",lwd=2, cex=2, 
     main="Dystrybuanta", cex.axis=0.8) #Gęstość

# losowanie wartości z populacji o różnych rozkładach 
require(graphics)
y <- rt(200, df = 5) # losowanie wartości z rozkładu t Studenta
z <- rweibull(200, shape=0.5) # losowanie wartości z rozkładu Weibulla
chi <- rchisq(200, df = 4) #losowanie wartości z rozkładu chi kwadrat
norm <- rnorm(200, 0, 1) # losowanie wartości z rozkładu normalnego
# porównanie różnych rozkładów za pomocą wykresu Quantile-Quantile plots
par(mfrow=c(1,2))
qqnorm(norm); qqline(norm, col = 2)
qqnorm(y); qqline(y, col = 2)


#################################################################################################
#Ćwiczenie 7. Porównaj ze sobą rozkład wartości losowanych z populacji
#o rozkładzie normalnym z wartościami z populacji o rozkładzie Weibulla.
#Użyj funkcji par(mfrow) do podzielenia okna na dwie kolumny.
par(mfrow=c(1,1))


#Przeanalizujmy teraz realne dane biologiczno-medyczne, przykładowe dane 
#zapisane w r

dane <- iris #pomiary kwiatów 3 gatunków kosaćca, więcej informacji po wpisaniu ?iris
head(dane) #wyświetlamy pierwsze wiersze zestawu danych
tail(dane) #wyświetlamy ostatnie wiersze zestawu danych

levels(dane$Species) # sprawdzamy jakie gatunki opisano w zestawie danych
names(dane) # sprawdzamy jakie mamy zmienne w zestawie danych

dim(dane) # sprawdzamy ile mamy wierszy i kolumn w danych
class(dane) #sprawdzamy jakiej klasy jest nasz obiekt dane

#Spróbujmy wyznaczyć histogram dla długości płatków korony Iris setosa
setosa <- subset(dane, dane$Species=="setosa") #wydzielanie danych tylko dla I.setosa
hist(setosa$Sepal.Length, main="I. setosa", ylab="Częstość")

#############################################################################################
#Ćwiczenie 8.Utwórz histogramy dla I. virginica oraz I. versicolor
#Porównaj je z użyciem par(mfrow)
versicolor <- subset(dane, dane$Species=="versicolor")
virginica <- subset(dane, dane$Species=="virginica")
# a teraz ocena rozkładów tych długości z użyciem quantile-quantile plot
qqnorm(versicolor$Petal.Length); qqline(versicolor$Petal.Length, col = 2)
qqnorm(virginica$Petal.Length); qqline(virginica$Petal.Length, col = 2)

#Czy jednak na pewno wiemy czy te dane mają rozkład normalny czy też nie?
# Ostateczną odpowiedź z punktu widzenia statystyki nie można udzielić 
# w 100% na podstawie niektórych wykresów. Stosuje się specjalne testy
# które pozwalają ocenić czy rozkład danej zmiennej jest normalny czy też nie.
# Jednym z nich jest test Kołmogorowa-Smirnowa - sprawdźmy czy powie on nam 
# czy długości płatków są zgodne z rozkładem normalnym u tych 3 gatunków kosaćca.
mean(setosa$Petal.Length)
sd(setosa$Petal.Length)
ks_setosa <- ks.test(setosa$Petal.Length, "pnorm", mean=1.462, sd=0.173664)
ks_setosa
#data:  setosa$Petal.Length
#D = 0.1534, p-value = 0.19
#alternative hypothesis: two-sided
# powyższe wartości wskazują, że nie ma podstaw do odrzucenia hipotezy 
#o zgodności z rozkładem normalnym
# p= 0.05      jeżeli   p = 0.19 > p=0.05 - wniosek: nasze dane zgodne z rozkładem normalnym
###########################################################################################
#Ćwiczenie 9. Sprawdź, czy długości płatków I. virginica oraz I. versicolor 
#odbiegają od rozkładu normalnego? Użyj testu Kołmogorova-Smirnova.

#Drugim znanym i powszechnie używanym testem zgodności z rozkładem normalnym 
#jest test Shapiro-Wilka. Nie trzeba obliczać dystrubuanty rozkładu normalnego.
sh_setosa <- shapiro.test(setosa$Petal.Length)
sh_setosa
#Shapiro-Wilk normality test
#data:  setosa$Petal.Length
#W = 0.95498, p-value = 0.05481 
#Wyniki testu wskazują, że również w tym przypadku nie ma podstaw do odrzucenia 
#hipotezy zerowej - jednakże było blisko - wartość p tylko nieznacznie przekracza 
#umowny poziom istotności alfa - jest tutaj równe 0.05481

############################################################################################
#Ćwiczenie 10. Sprawdź, czy długość płatków pozostałych gatunków kosaćca ma 
#rozkład normalny.Użyj testu Shapiro-Wilka