# 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")