# sprawdzamy dostępność potrzebnych pakietów: mlbench, boot, Hmisc
      # jeśli któregoś brakuje, instalujemy
check.inst <- c("mlbench", "boot", "Hmisc") %in% row.names(installed.packages())
if (! check.inst[1])
  install.packages("mlbench")
if (! check.inst[2])
  install.packages("boot")
if (! check.inst[3])
  install.packages("Hmisc")

      # ładujemy pakiet mlbench ze zbiorami danych
library(mlbench)
      # obejrzymy zbiór Glass
data(Glass)

      # charakterystyka atrybutu RI
summary(Glass$RI)
      # średnia, odchylenie standardowe, wariancja
mean(Glass$RI)
sd(Glass$RI)
var(Glass$RI)
      # policzmy wariancję sami
myvar <- function(x) { sum((x-mean(x))^2)/(length(x)-1) }
      # i zobaczmy jak bardzo się różni
myvar(Glass$RI)-var(Glass$RI)

      # mediana i kwartyle
median(Glass$RI)
quantile(Glass$RI)
      # pierwszy kwartyl
quantile(Glass$RI)[2]
quantile(Glass$RI, probs=0.25)
      # trzeci kwartyl
quantile(Glass$RI)[4]
quantile(Glass$RI, probs=0.75)

      # parę przykładów "posługiwania się" rozkładem normalnym
      # dla jakiego u prawdopodobieństwo wartości poniżej u wynosi 0.025
qnorm(0.025)
      # dla jakiego u prawdopodobieństwo wartości powyżej u wynosi 0.025
qnorm(1-0.025)
      # jakie jest prawdopodobieństwo wartości poniżej 1.96
pnorm(1.96)
      # jakie jest prawdopodobieństwo wartości między -1.96 a 1.96
pnorm(1.96)-pnorm(-1.96)
pnorm(1.96)-(1-pnorm(1.96))

      # zobaczmy, czy rozkład atrybutu RI wygląda na normalny
qqnorm(Glass$RI)
      # to samo dla wszystkich atrybutów z wyjątkiem ostatniego
par(mfrow=c(3,3))
for (a in names(Glass)[-10]) qqnorm(Glass[[a]], main=a)

      # detekcja wartości odstających
is.outlier <- function(x, b=1.5) { q <- quantile(x); iqr <- q[4]-q[2]; x < q[2]-b*iqr | x > q[4]+b*iqr }
      # ile jest wartości odstających atrybutu RI?
sum(is.outlier(Glass$RI))

      # boxplot
boxplot(Glass$RI)
      # porównajmy naszą detekcję wartości odstających z tą boxplotową
boxplot(Glass$RI)$out
Glass$RI[is.outlier(Glass$RI)]
boxplot(Glass$RI)$out == Glass$RI[is.outlier(Glass$RI)]

      # histogram
hist(Glass$RI)
      # można zmienić przedziały
hist(Glass$RI, breaks=seq(1.510, 1.535, 0.0025))

      # charakterystyka atrybutu dyskretnego Type
summary(Glass$Type)
table(Glass$Type)

      # wyznaczanie mody (dominanty) -- niezbyt eleganckie, ale krótkie
modus <- function(x) { names(table(x))[table(x)==max(table(x))] }
modus(Glass$Type)

      # m-estymacja prawdopodobieństwa wartości v w wektorze x
m.est <- function(x, v, m=length(levels(x)), p=1/m) { (sum(x==v, na.rm=T)+m*p)/(length(x)+m) }
m.est(Glass$Type, "1")
m.est(factor(1:100), 100)
many1.single2 <- factor(c(rep(1, 99), 2))
m.est(many1.single2, 1)
m.est(many1.single2, 2)

      # załadujmy dane HouseVotes84
data(HouseVotes84)
      # pakiet boot zawiera funkcję wyznaczania przedziału ufności
library(boot)
      # przedział ufności dla średniej atrybutu RI  przy założeniu rozkładu normalnego
      #   t0 -- parametr, dla którego wyznaczamy przedział (tu średnia)
      #   t0.var -- jego wariancja (tu wariancja średniej)
norm.ci(t0=mean(Glass$RI), var.t0=var(Glass$RI)/length(Glass$RI))
      # dla porównania przedział ufności wyznaczany metodą bootstrap
boot.ci(boot(Glass$RI, weighted.mean, 100, stype="w"), type="basic")

      # prosty klasyfikator z przykładów do wykładu 1
ruleV4 <- function(data) { factor(ifelse(data$V4 == "y", "republican", "democrat"), levels=c("democrat", "republican")) }
      # przedział ufności dla błędu klasyfikacji (prawdopodoństwa pomyłki)
n1 <- sum(ruleV4(HouseVotes84)!=HouseVotes84$Class, na.rm=T)
n <- sum(!is.na(ruleV4(HouseVotes84)))
p <- n1/n
norm.ci(t0=p, var.t0=p*(1-p)/n)

      # sprawdzimy to jeszcze funkcją z pakietu HMisc
library(Hmisc)
binconf(n1, n, method="all")