# Przygotowanie danych
rm(list=ls())
colors = c("red", "green", "blue")
iris.train = read.csv2("irysy_trening.csv", row.names = 1)
iris.test = read.csv2("irysy_prognoza.csv", row.names = 1)
names(iris.train) = c("SL", "SW", "PL", "PW", "Species")
names(iris.test) = c("SL", "SW", "PL", "PW", "Species")
attach(iris.train)

### Wielowymiarowa analiza wariancji
# (Dla jednowymiarowej anovej jest przygotowany poprzedni skrypt)

# Sprawdzenie jak bardzo, jako całość, różnią się PL, PW, SL, SW w grupach
model = manova(cbind(PL, PW, SL, SW) ~ Species)
summary(model)  # statystyka Lambda jest przerabiana na F (approx F)
summary(model, "Wilks") # można uzyskać statystykę Wilka
# Ogólnie powyższe statystyki pozwalają nam ocenić, czy statystyki różnią się między grupami
# Tzn każda funckja dyskryminacyjna odpowiada wektorowi własnemu pewnej macierzy
# A hipotezą zerową jest, czy wszystkie wartości własne to 0
# Bo jeśli tak jest to wartości PL, PW, SL, SW nie różnią się między grupami

# Sprawdzenie, dla jakich współczynnikach wartość statystyki F będzie największa
library(MASS)
model.lda = lda(cbind(PL, PW), grouping = Species)
lda(Species ~ PL + PW) # znowu, to jest ten sam sposób zapisania tego co wyżej
# Zauważymy, że współczynniki przy LD1 to 1.33 oraz 2.63
summary(aov(PL + PW ~ Species)) # małe F
summary(aov(1.331925*PL + 2.630017*PW ~ Species)) # największe możliwe F
# więc zmienna postaci 1.331925*PL + 2.630017*PW najlepiej rozróżnia grupy

# Test: Sprawdzamy, czy wartości własne (posortowane malejąco) od k+1-ego miejsca są równe 0
# Czyli jak wybieramy k = 1, to sprawdzamy czy od 2 miejsca wzwyż wszystki wartości własne są równe 0
# H0: wartość własna k+1 = wartość własna k+2 = ... = wartość własna n
# Ha: w przeciwnym przypadku
# Statystyka lambda = floor(n - 1 + (p+g)/2)*{suma od k+1 do r}{ln(1+wartość własna i-ta)}
# Stopnie swobody: (p-k)(n-k-1) (na wykładzie jest (p-k)(g-k-1) - prawdopodobnie źle)
# Oznaczenia: g - ilośc grup; p - ilość zmiennych losowych oruginalnych; k - sami wybieramy; n - liczba rekordów (wierszy, obserwacji)

a = floor(model.lda$N - 1 + (4 + 3) / 2) # floor(n - 1 + (p+g)/2)
chi = a * log(1+model.lda$svd) # *{suma od k+1 do r}{ln(1+wartość własna i-ta)}
chi[2] # dla k = 1 sumujemy tylko jeden wyraz, czyli wybieramy ten
st.swobody = (4-1) * (60 - 2) # (p-k)(n-k-1)
# 1 sposób:
1 - pchisq(chi[2], st.swobody) # wychodzi nam, p.value= 1, więc nie można odrzucić hipotezy, że wartość własna dla LD2 wynosi 0, czyli LD2 nie ma znaczenia
# 2 sposób:
chi[2] < qchisq(0.95, 174) # więc nie można odurzucić hipotezy zerowej

### LDA (analiza dyskryminacyjna)
model.lda = lda(Species ~ ., data = iris.train)
prediction = predict(model.lda)
prediction$x[,1] -> LD1
prediction$x[,2] -> LD2
plot(LD1, LD2, col=colors[Species]) #wykres pokazujący, jak wyglądają funkcje dyskryminacyjne
plot(model.lda) # tu dostajemy to samo, tylko brzydko
ldahist(LD1, Species) #histogramy dla każdej z grup względem funkcji LD1, tej najlepszej
ldahist(LD2, Species) #a tu względem tej najgorszej

# Parametry LDA:
model.lda$svd # wartości własne odpowiadającem funkcjom LD1, LD2
model.lda$means # średnie w grupach
model.lda$scaling # parametry funkcji
predict(model.lda, as.data.frame(model.lda$means))$x # średnie w grupach dla funkcji LD1, LD2
predict(model.lda, iris.test)$x # wartościu funkcji dyskryminacyjnych dla danych testowych
predictions.test = predict(model.lda, iris.test)$class # predykcje za pomocą LDA
table(predictions.test, iris.test$Species) # sprawdzenie skuteczności
predict(model.lda, iris.test, prior = c(0.5, 0.49, 0.01))$class # można wspomagać się znając prawdopodobieństwa należenia do grup
# Uwaga: nie stosować powyższego sposobu, gdy jedna z wartości będzie skrajnie mała, np. 0.00001

# Test: Sprawdzamy, czy zmienna jest istotnie statystyczna w LDA
# H0: Statystyka jest nieistotna
# Ha: Statystyka jest istotna
# Statystyka: F = (n-g-j) / (g-1) * ( lambda_bez_zmiennej / lambda_ze_zmienna - 1)
# Gdzie: n - ilość wierszy, obserwacji; g - ilość grup;
# j - ilośc zmiennych objaśniających (ale bez tej którą sprawdzamy, czyli wszystkie zminne objaśniające - 1)
# lambdy - statystyki wilka
# Stopnie swobody: n-g-j, g-1
library(rrcov)
vine.train = read.csv2("wino.trening.csv", row.names = 1)
Wilks.test(typ ~ ., data = vine.train)$wilks-> L13
Wilks.test(typ ~ . - proline, data = vine.train)$wilks -> L12
F <- 30 / 2 * (L12 / L13 - 1) # 30  = 45 - 3 - (13 - 1)
F > qf(p=0.95, df1 = 2, df2 = 30)
# więc odrzucamy hipotezę zerową, nie można powiedzieć, że zmienna proline nie ma wpływu na gatunek wina