getwd()
## [1] "/Users/vladimirbenic/Dropbox/math/probability/R/rossiter/scripts"
setwd("rossiter/scripts/")
## Error: cannot change working directory
mhw = read.csv("../mhw.csv")
Pogledajmo strukturu data framea mhw.
str(mhw)
## 'data.frame': 500 obs. of 4 variables:
## $ r : int 1 2 3 4 5 6 7 8 9 10 ...
## $ c : int 1 1 1 1 1 1 1 1 1 1 ...
## $ grain: num 3.63 4.07 4.51 3.9 3.63 3.16 3.18 3.42 3.97 3.4 ...
## $ straw: num 6.37 6.24 7.05 6.91 5.93 5.59 5.32 5.52 6.03 5.66 ...
Grafički prikažimo točke sa kordinatama grain i straw.
plot(mhw$straw ~ mhw$grain,
xlab = "Prinos zrna",
ylab = "Prinos slame",
pch = 20, col="darkblue", cex=1.2,
sub="Median je prikazan pomocu crvenih isprekidanih linija")
title (main="Slama vs. zrnevlje, Mercet-Hall eksperiment")
grid()
abline(v=median(mhw$grain), lty=2, col="red")
abline(h=median(mhw$straw), lty=2, col="red")
Nacrtajmo “istu” sliku pomoću paketa ggplot2.
ggplot(mhw, aes(grain, straw)) +
geom_point(col="darkblue", size=3) +
ggtitle("Slama vs. zrnevlje, Mercet-Hallov eksperiment") +
geom_hline(aes(yintercept=median(straw)), col="red", lty=2) +
geom_vline(aes(xintercept=median(grain)), col="red", lty=2) +
xlab("Prinos zrna\nCrveno je naznacen polozaj mediana") + ylab("Prinos slame") +
theme_bw()+
theme(plot.title = element_text(face="bold", colour="#990000", size=20))
Dva su najznačajnija mjerila korelacije
Ovdje ćemo upotrebiti Pearsonov koeficijent. Ako je dvodimenzionala razdioba normalna tada taj koeficijent ima neka dobra svojstva. To je procjenitelj korelacije populacije po metodi maksimalne vjerodostojnosti, zatim je asimptotski nepristrasan i efikasan procjenitelj.
Stoga ćemo naša mjerenja usporediti sa uzorkom dvodimenzionalne normalne razdiobe sa istim parametrima kao što su parametri empirijske razdiobe.
Zadatak: Vizualiziraj dvodimenzionalnu normalnu razdiobu sa parametrima iz prinosa zrna i slame; vizualno usporedi uzorak sa podacima iz mhw.
Za dobivanje podataka dvodimenzionalne korelirane normalne razdiobe
koristiti ćemo funkciju mvrnorm iz biblioteke MASS. Za to nam
trebaju podaci o srednjim vrijednostima i variance-covariance matrica.
Vrijednosti srednjih vrijednosti ćemo dobiti pomoću funkcije colMeans().
mu = colMeans(mhw[, c("grain", "straw")])
mu
## grain straw
## 3.949 6.515
Varijace-covariance matricu ćemo dobiti pomoću vektorizirane var()
funkcije
sigma = var(mhw[, c("grain", "straw")])
sigma
## grain straw
## grain 0.2100 0.3004
## straw 0.3004 0.8070
Sada mižemo generirati uzorak iz dvodimenzionalne korelirane normalne razdiobe.
require(MASS)
## Loading required package: MASS
sim.sample = mvrnorm(length(mhw$grain), mu = mu, Sigma = sigma)
head(sim.sample)
## grain straw
## [1,] 3.792 7.125
## [2,] 4.497 7.498
## [3,] 4.536 8.350
## [4,] 3.719 6.222
## [5,] 3.503 5.377
## [6,] 4.276 6.721
Usporedimo osnovne podatke genriranih i mjerenih veličina.
summary(sim.sample)
## grain straw
## Min. :2.38 Min. :3.00
## 1st Qu.:3.63 1st Qu.:5.95
## Median :3.98 Median :6.47
## Mean :3.97 Mean :6.52
## 3rd Qu.:4.29 3rd Qu.:7.14
## Max. :5.33 Max. :9.05
summary(mhw[, c("grain", "straw")])
## grain straw
## Min. :2.73 Min. :4.10
## 1st Qu.:3.64 1st Qu.:5.88
## Median :3.94 Median :6.36
## Mean :3.95 Mean :6.51
## 3rd Qu.:4.27 3rd Qu.:7.17
## Max. :5.16 Max. :8.85
par(mfrow = c(2, 2))
grain.lim = c(min(sim.sample[, "grain"], mhw$grain ),
max(sim.sample[, "grain"], mhw$grain))
straw.lim = c(min(sim.sample[, "straw"], mhw$straw),
max(sim.sample[, "straw"], mhw$straw))
hist(mhw$grain, xlim = grain.lim, main = "Grain (actual)",
col = "lightyellow", breaks = seq(grain.lim[1], grain.lim[2],
length = 17))
hist(sim.sample[, "grain"], xlim = grain.lim, main = "Grain (simulated)",
col = "antiquewhite", breaks = seq(grain.lim[1], grain.lim[2],
length = 17))
hist(mhw$straw, xlim = straw.lim, main = "Straw (actual)",
col = "lightblue", breaks = seq(straw.lim[1], straw.lim[2],
length = 17))
hist(sim.sample[, "straw"], xlim = straw.lim, main = "Straw (simulated)",
col = "lavender", breaks = seq(straw.lim[1], straw.lim[2],
length = 17))
par(mfrow = c(1, 1))
Nacrtajmo “istu” sliku pomoću bibilioteke ggplot2.
Za crtanje više figura na jednoj slici treba nam funkcija multiplot()
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols: Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
require(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
ssdf = data.frame(sim.sample)
grain.lim = c(min(sim.sample[, "grain"], mhw$grain ),
max(sim.sample[, "grain"], mhw$grain))
straw.lim = c(min(sim.sample[, "straw"], mhw$straw),
max(sim.sample[, "straw"], mhw$straw))
binwidth.grain = (grain.lim[2] - grain.lim[1])/16
binwidth.straw = (straw.lim[2] - straw.lim[1])/16
p1 = ggplot(mhw, aes(grain)) +
geom_histogram(binwidth=binwidth.grain, col="black", fill="lightyellow") +
coord_cartesian(xlim = grain.lim) +
ggtitle("Grain (actual)") + theme_bw() +
theme(panel.grid.major=element_line(size = 0.3, color="black",
linetype = "dotted"))
p2 = ggplot(ssdf, aes(grain)) +
geom_histogram(binwidth=binwidth.grain, col="black", fill="antiquewhite") +
coord_cartesian(xlim = grain.lim) +
ggtitle("Grain (simulated)") + theme_bw()+
theme(panel.grid.major=element_line(size = 0.3, color="black",
linetype = "dotted"))
p3 = ggplot(mhw, aes(straw)) +
geom_histogram(binwidth=binwidth.straw, col="black", fill="lightblue") +
coord_cartesian(xlim = straw.lim) +
ggtitle("Straw (actual)") + theme_bw()+
theme(panel.grid.major=element_line(size = 0.3, color="black",
linetype = "dotted"))
p4 = ggplot(ssdf, aes(straw)) +
geom_histogram(binwidth=binwidth.straw, col="black", fill="lavender") +
coord_cartesian(xlim = straw.lim) +
ggtitle("Straw (simulated)") + theme_bw() +
theme(panel.grid.major=element_line(size = 0.3, color="black",
linetype = "dotted"))
multiplot(p1, p2, p3, p4, cols=2)
## Loading required package: grid
\[ \phi(x) = \int_0^x \frac{t\,dt}{1+\sin^3(t)} \]
library(xtable)
print(xtable(head(mhw, 6)), type = "html", include.rownames = F)
| r | c | grain | straw |
|---|---|---|---|
| 1 | 1 | 3.63 | 6.37 |
| 2 | 1 | 4.07 | 6.24 |
| 3 | 1 | 4.51 | 7.05 |
| 4 | 1 | 3.90 | 6.91 |
| 5 | 1 | 3.63 | 5.93 |
| 6 | 1 | 3.16 | 5.59 |
par(mfrow = c(1, 2))
plot(sim.sample, main = "Slama vs. zrnevlje (simulirano)",
xlab = "Zrnevlje", ylab="Slama",
xlim = grain.lim, ylim = straw.lim,
pch = 20, col= "blue" )
abline(v = median(sim.sample[, 1]), lty = 2, col = "red")
abline(h = median(sim.sample[, 2]), lty = 2, col = "red")
plot(mhw$grain, mhw$straw, main = " Slama vs. zrno (izmjereno)",
xlab = "Zrnevlje", ylab="Slama",
xlim = grain.lim, ylim = straw.lim,
pch = 20, col= "black" )
abline(v = median(mhw$grain), lty = 2, col = "red")
abline(h = median(mhw$straw), lty = 2, col = "red")
par(mfrow = c(1,1))
Analognu sliku nacrtajmo pomoću biblioteke ggplot2.
p1 = ggplot(ssdf, aes(grain, straw)) +
geom_point(color = "blue") +
labs(title="Slama vs zrno (simulirano)", x = "Zrno", y = "Slama") +
geom_hline(aes(yintercept=median(straw)), col="red", lty=2) +
geom_vline(aes(xintercept=median(grain)), col="red", lty=2) +
theme_bw() +
theme(panel.grid.major=element_line(size = 0.3, color="black",
linetype = "dotted"),
plot.title = element_text(face="bold", colour="#990000", size=15),
axis.title.x = element_text(face="bold", size=12),
axis.title.y = element_text(face="bold", size=12))
p2 = ggplot(mhw, aes(grain, straw)) +
geom_point(color = "black") +
labs(title="Slama vs zrno (izmjereno)", x = "Zrno", y = "Slama") +
geom_hline(aes(yintercept=median(straw)), col="red", lty=2) +
geom_vline(aes(xintercept=median(grain)), col="red", lty=2) +
theme_bw() +
theme(panel.grid.major=element_line(size = 0.3, color="black",
linetype = "dotted"),
plot.title = element_text(face="bold", colour="#990000", size=15),
axis.title.x = element_text(face="bold", size=12),
axis.title.y = element_text(face="bold", size=12))
multiplot(p1, p2, cols=2)
Budući da smo gotovi sa vizualnom usporedbom izbrisati cemo nepotrebne objekte.
rm(sim.sample, grain.lim, straw.lim, binwidth.grain, binwidth.straw)
Kako se iz predhodnog vidi empirijska dvodimenzionalna razdioba zrna i slame je približno normalna pa će Pearsonov koeficijent korelacije imati mnoga dobra svojstva.
Zadatak Izračunajmo koeficijent korelacije između zrnevlja i slame i njegov 95% interval povjerenja.
cor(mhw$grain, mhw$straw)
## [1] 0.7298
Usporedimo dobiveni koeficijent korelacija sa rezultatom simuliranih podataka.
cor(ssdf$grain, ssdf$straw)
## [1] 0.7469
Vidimo da su dobiveni rezultati slični.
Odredimo interval povjerenja korelacije.
rez = cor.test(mhw$grain, mhw$straw)
rez
##
## Pearson's product-moment correlation
##
## data: mhw$grain and mhw$straw
## t = 23.82, df = 498, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6860 0.7683
## sample estimates:
## cor
## 0.7298
Ako nas zanima samo interval povjerenja a ne cjelokupni output funkcije cor.test()
izvedemo naredbu
rez$conf.int
## [1] 0.6860 0.7683
## attr(,"conf.level")
## [1] 0.95