sadII_lab4_rciagle

2429 days ago by macieksk

# Tryb %r ^ #Set output widths better .adjustWidth <- function(...){ options(width=10000) ; TRUE} .adjustWidthCallBack <- addTaskCallback(.adjustWidth) 
       
.adjustWidth <- function(...){ options(width=10000) ; TRUE} 
#Pakiet zawierajacy funkcje rdirichlet #?rdirichlet #install.packages("MCMCpack") library(MCMCpack) 
       
library(lattice) D1<-as.data.frame(rdirichlet(1000,c(13,2))) densityplot( ~ V2, data=D1) a=dev.off() head(D1) 
       
There were 24 warnings (use warnings() to see them)

         V1         V2
1 0.7901627 0.20983728
2 0.8373695 0.16263054
3 0.8873720 0.11262795
4 0.8997226 0.10027744
5 0.7304461 0.26955387
6 0.9072827 0.09271729
plot(dbeta(seq(0,1,by=0.1),2,13)) dev.off() 
       
Warning messages:
1: In plot.xy(xy, type, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
2: In plot.xy(xy, type, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
3: In axis(side = side, at = at, labels = labels, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
4: In axis(side = side, at = at, labels = labels, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
5: In axis(side = side, at = at, labels = labels, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
6: In axis(side = side, at = at, labels = labels, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
7: In box(...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
8: In box(...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
9: In title(...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
10: In title(...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
null device 
          1 
# Co to za rozklad? # http://en.wikipedia.org/wiki/Dirichlet_distribution#Marginal_distributions # Przetestuj na niego testem Chi2 
       
# Wykonaj dwa testy Chi2: # 1) Uzywajac qbeta(...) ze znanymi Ci parametrami (odczytaj je odpowiednio z parametrow uzytych do wylosowania probki z r.Dirichleta) # 2) Wpierw dopasuj odpowiedni rozklad Beta do danych uzywajac np. fitdistr. Wykonaj Chi2 do rozkladu z dopasowanymi parametrami. # Jakiego parametru degrees.of.freedom powinnas/powinienes uzyc? 
       
p<-seq(0,1,by=0.1) 
       
       
 [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
p<-seq(0,1,by=0.1) diff(pbeta(p,13,2)) 
       
 [1] 1.270000e-12 9.337610e-09 1.600927e-06 5.744553e-05
8.564715e-04 7.182103e-03 3.937798e-02 1.504365e-01 3.867170e-01
4.153709e-01
summary(cut(D1$V1,p)) 
       
  (0,0.1] (0.1,0.2] (0.2,0.3] (0.3,0.4] (0.4,0.5] (0.5,0.6]
(0.6,0.7] (0.7,0.8] (0.8,0.9]   (0.9,1] 
        0         0         0         0         1         2       
42       151       374       430 
D2<-as.data.frame(rdirichlet(1000,c(5,2,1))) xyplot(V1 ~ V2, data=D2) a=dev.off() head(D2) 
       
There were 22 warnings (use warnings() to see them)

         V1        V2         V3
1 0.3323772 0.2437777 0.42384506
2 0.4561205 0.1389690 0.40491056
3 0.5246771 0.4375764 0.03774650
4 0.6612822 0.3249982 0.01371958
5 0.2570513 0.2645281 0.47842062
6 0.6079470 0.3164471 0.07560591
# Przeczytaj o tescie Kolmogorov-Smirnov # ?ks.test # 1. Przetestuj zmiennÄ… D2$V2 na odpowiednia dystrybucje uzywajac ks.test z argumentem y="pbeta" i odpowiednimi pozostalymi parametrami # 2. Wykonaj KS test na dwoch probkach: jako parametr y= podaj wektor wylosowany z rbeta z odpowiednimi parametrami 
       
### Testy na normalnosc ## Generate two data sets ## First Normal, second from a t-distribution X1 = rnorm(100); X2 = rt(100, df=3) suppressWarnings({ ## Have a look at the densities plot.density(density(X1)); plot.density(density(X2)) ## Plot using a qqplot qqnorm(X1);qqline(X1, col = 2) qqnorm(X2);qqline(X2, col = 2) }) ## Perform the test shapiro.test(X1); shapiro.test(X2) dev.off() 
       
	Shapiro-Wilk normality test

data:  X1 
W = 0.9803, p-value = 0.1409


	Shapiro-Wilk normality test

data:  X2 
W = 0.7779, p-value = 5.481e-11

null device 
          1 



# Uzywajac testow z biblioteki nortest library(nortest) # http://cran.r-project.org/web/packages/nortest/nortest.pdf # Sprawdz zaimplementowane tam testy na czulosc i specyficznosc. # W tym celu wykonaj na losowych probkach wielkosci 100: # 1. rnorm(100) # 2. rt(100,df=3) # ... po 100(lub wiecej) roznych testow i zanotuj p-value. # (Podobnie jak na poprzednich zajeciach mozesz wykonac test dla roznych wartosci df i wielkosci probki, ale nie jest to konieczne :) #norm.test.data<-... 
       
norm.test.data <- data.frame(norm="ad.test",t="ad.test",stringsAsFactors=FALSE) norm.test.data <- rbind(norm.test.data, rep("cvm.test",2)) norm.test.data <- rbind(norm.test.data, rep("lillie.test",2)) norm.test.data <- rbind(norm.test.data, rep("pearson.test",2)) norm.test.data <- rbind(norm.test.data, rep("sf.test",2)) norm.test.data norm.test.data<-stack(norm.test.data) colnames(norm.test.data) <- c("rodzaj.testu", "zrodlo.danych") norm.test.data 
       
          norm            t
1      ad.test      ad.test
2     cvm.test     cvm.test
3  lillie.test  lillie.test
4 pearson.test pearson.test
5      sf.test      sf.test


   rodzaj.testu zrodlo.danych
1       ad.test          norm
2      cvm.test          norm
3   lillie.test          norm
4  pearson.test          norm
5       sf.test          norm
6       ad.test             t
7      cvm.test             t
8   lillie.test             t
9  pearson.test             t
10      sf.test             t
norm.test.data <- do.call(rbind, replicate(1000,norm.test.data,simplify=FALSE) ) table(norm.test.data ) 
       
              zrodlo.danych
rodzaj.testu   norm    t
  ad.test      1000 1000
  cvm.test     1000 1000
  lillie.test  1000 1000
  pearson.test 1000 1000
  sf.test      1000 1000
do.test<-function(size)function(test.fun,zr) { if(zr=="norm") s<-rnorm(size) if(zr=="t") s<-rt(size,df=2) get(as.character(test.fun))(s)$p.value } 
       
Vectorize(do.test(20))("ad.test","t") 
       
  ad.test 
0.1043804 
norm.test.alldata<-transform(norm.test.data, p.value=Vectorize(do.test(20))(rodzaj.testu,zrodlo.danych) ) 
       
library(lattice) png(width=1000,height=500) histogram(~p.value|rodzaj.testu+zrodlo.danych, data=norm.test.alldata, type="percent", breaks=20, scale=list(x=list(log=2))) dev.off() 
       
There were 50 or more warnings (use warnings() to see the first 50)
null device 
          1 
# Narysuj krzywa ROC, czyli sensitivity vs specificity # http://en.wikipedia.org/wiki/Receiver_operating_characteristic # Koniecznie popatrz na ten obrazek # http://en.wikipedia.org/wiki/File:Receiver_Operating_Characteristic.png # Wygeneruj 100 punktow krzywej ROC: # pval.thresh<-seq(0,1,by=0.01) 
       
# W tym celu zaprogramuj nast. funkcje get.point.for.ROC.curve<-function(class.vector, # Wektor 0 i 1 w zaleznosci od tego, czy dana jest Negative , czy Positive parameter.vector, # Wektor parametru (np. pvalue) po ktorym dokonujemy odciecia threshold ) # Punkt decyzji: parameter.vector<threshold ==> decyzja na Positive { tpr<- ... fpr<- ... # return c( true.positive.rate= tpr, # True.Positives / Positives false.positve.rate= fpr, # False.Positives / Negatives specificity = 1-fpr, # True.Negatives / Negatives sensitivity = tpr) # Dla jasnosci } 
       
# Zaaplikuj te funkcje do danych o testach na normalnosc: norm.test.data # i narysuj wykres ROC pval.thresh <-seq(0,1,by=0.01) roc.data <- sapply(pval.thresh, function(thr)get.point.for.ROC.curve(... ... , thr) ) 
       
## .... albo uzyj pakietu ROCR # # http://rocr.bioinf.mpi-sb.mpg.de/ 
       
#install.packages("ROCR") library(ROCR) 
       
pred<- with(subset(norm.test.alldata,rodzaj.testu=="ad.test"),{ pred<-prediction(-p.value,zrodlo.danych) pred }) perf <- performance(pred, measure = "tpr", x.measure = "fpr") plot(perf, colorize=TRUE, col=rainbow(10)) dev.off() 
       
There were 50 or more warnings (use warnings() to see the first 50)
null device 
          1 
png(width=600,height=600) perf <- performance(pred, measure = "sens", x.measure = "spec") plot(perf, col=rainbow(10),colorize=TRUE, print.cutoffs.at=-seq(0,1,by=0.1)^3, text.adj=c(1.2,1.2),lwd=3, # Averaging over several runs avg="threshold", spread.estimate="boxplot", show.spread.at=-seq(0,1,by=0.1)^3, ) dev.off() 
       
There were 50 or more warnings (use warnings() to see the first 50)

null device 
          1 
png(width=600,height=600) lvls <- levels(factor(norm.test.alldata$rodzaj.testu)) colors <- rainbow(length(lvls)) a<-mapply(function(i,lvl){ pred<- with(subset(norm.test.alldata, rodzaj.testu == lvl), prediction(-p.value,zrodlo.danych) ) perf <- performance(pred, measure = "tpr", x.measure = "fpr") plot(perf, add = i!=1, col=colors[i], colorize=FALSE, print.cutoffs.at=-seq(0,1,by=0.1)^3, text.adj=c(1.2,1.2),lwd=3) }, seq(along=lvls),lvls) legend(0.8,0.8, lvls, fill=colors) dev.off() 
       
There were 40 warnings (use warnings() to see them)
Warning messages:
1: In rect(left, top, r, b, angle = angle, density = density, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
2: In rect(left, top, r, b, angle = angle, density = density, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
3: In rect(left, top, r, b, angle = angle, density = density, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
4: In rect(left, top, r, b, angle = angle, density = density, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
5: In text.default(x, y, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
6: In text.default(x, y, ...) :
  X11 protocol error: BadMatch (invalid parameter attributes)
null device 
          1 
perf <- performance(pred, "pcmiss","lift") plot(perf,colorize=TRUE,col=rainbow(10)) dev.off() 
       
There were 50 or more warnings (use warnings() to see the first 50)
null device 
          1 
 
       
# Przeczytaj o rozkladzie Gammma str. 96 # http://biecek.pl/RinHasselt/manuals/RandStatisticsWithBiologicalExamples.pdf # Wykonaj wykres ze strony 98 za pomoca pakietu lattice # Nastepnie przeczytaj o rozkladzie Beta # Wykonaj podobny wykres dla rozkladu Beta dla roznych parametrow rozkladu: a=seq(0,3,by=0.5), b=seq(0,4,by=0.5) 
       
 
       
# Ze skryptu http://biecek.pl/RinHasselt/manuals/RandStatisticsWithBiologicalExamples.pdf # przeczytaj o Wielowymiarowym rozkladzie normalnym, str. 127 
       
 
       
# Przeczytaj i wykonaj cwiczenia ze skryptu dotyczace The Multinomial Distribution # str.121 # Jesli masz czas, do danych ze skryptu dopasuj model Multinomial # library(nnet) # fit.mn <- multinom(, data=twoje.dane)