############################################ # Test pro parametr p v binomickem rozdeleni ############################################ # # Uvazujme n bernoulliovskych alternativnich nezavislych pokusu # (napr. pri hodu kostkou sledujeme, zda padla sestka). # Uspech nastava s pravdepodobnosti p (v nasem pripade 1/6). # Oznacme X_n relativni pocet uspechu. CLV rika, ze X_n ma priblizne # (pro velka n) normalni rozdeleni se stredni hodnotou p a rozptylem p(1-p)/n. # # Interval spolehlivosti pro neznamy parametr p zalozeny na CLV je # (X_n - qnorm(1-alpha/2)*sqrt(X_n*(1-X_n)/n), # X_n + qnorm(1-alpha/2)*sqrt(X_n*(1-X_n)/n)). # Nulova hypoteza: p=p_0 # Test zalozeny na CLV probiha tak, ze nulovou hypotezu zamitneme, # kdyz abs(X_n-p_0) >= qnorm(1-alpha/2)*sqrt(X_n*(1-X_n)/n) # # Lze ovsem zkonstruovat presny interval spolehlivosti a presny test. # Vyuzijeme toho, ze n*X_n ma binomicke rozdeleni s parametry (n,p). # # Ukol: pomoci "binom.test" testujte nulovou hypotezu, ze pravdepodobnost # padnuti sestky na kostce je 1/6, pokud v 600 hodech padlo 80 sestek. # Sestrojte interval spolehlivosti pro p pomoci funkce "qbinom". # ######################## # chi^2-test dobre shody ######################## # # Ukol: predpokladejme, ze pri 600 hodech sestistennou kostkou byly # dosazene tyto empiricke cetnosti: 103,99,91,108,119,80 # Lze kostku povazovat za symetrickou? Pouzijte "chisq.test". # ############################## # # Krome testu, ktere jsme uz potkali, nabizi Rko celou radu dalsich # statistickych testu. # Pomoci help(package=stats) ziskame informace o vsech statistickych funkci # v knihovne 'stats', je mezi nimi rada testu. # ############################## # # Ukol: Mezi testy chybi Fisheruv podmineny test shody s Poissonovym rozdelenim. # Napiste jednoduchou funkci, ktera vypocte testovou statistiku # a spolu s ni vypise kriticke hodnoty prislusneho chi^2 rozdeleni. # # Hint1: Test je zalozen na faktu, ze Q = (n-1)*S^2/{\bar X} ma asymptoticky # chi^2 rozdeleni s n-1 stupni volnosti. V praxi se doporucuje asymptoticke # rozdeleni pouzivat, pokud {\bar X}>5. # # Hint2: Pocet pozorovani ve vektoru X se zjisti funkci length(X), # vyberovy rozptyl S^2 je var(X) a vyberovy prumer {\bar X} je mean(X). # # Hint3: funkce je opet objekt, takze to bude vypadat nejak takhle: # test.poisson <- function(data, alpha=0.05){ # statistika <- ... hodnota test. statistiky Q # krit.hodn <- ... prislusne kriticke hodnoty # vystup <- c(statistika, krit.hodn) # return(vystup) # } #end function # Takto napsana funkce bude na vstupu ocekavat 2 parametry, z toho jeden (data) povinne. # Pokud ten bude chybet, nahradi ho implicitni hodnotou 0.05. # Vystupem je vektor, jehoz prvni slozka je test. statistika a druha je vektor kritickych hodnot. # Jako vystup muzete rovnez pridat rozhodnuti, zda hypotezu zamitame, nebo nezamitame. # # Hint4: funkci si pripravte v nejakem textovem editoru (uplne staci Notepad, # Emacs apod., ale lze pouzit i editor Rka) a do Rka zkopirujte/nactete jako # celek, tj. neodesilejte zvlast kazdy radek kodu. Jinak se to ponekud spatne # ladi :-)) # # Svou funkci pak vyzkousejte na datech navstevy.txt udavajicich pocet navstev # jedne www stranky behem nekolika dnu. # K nacteni dat do vektoru staci pouzit funkci "scan", napr. data <- scan("navstevy.txt") # ############################## # # Jina moznost, jak testovat Poissonovo rozdeleni je nasledujici (Pearson): # rozdelime vyber do trid o dostatecne velkem poctu pozorovani, odhadneme # intenzitu a provedeme chi^2-test dobre shody. lambda <- mean(data) # odhad intenzity (maximalne verohodny) n <- length(data) rozdel <- function(x) { # abychom si rozdelili data do trid v <- x if (x<=8) v <- 8 # pocty <=8 sloucime do jedne tridy if (x>=22) v <- 22 # pocty >=22 do jedne tridy return(v) } tridy <- sapply(data,rozdel) # kdyz chceme aplikovat nejakou funkci # najednou na vektor, hodi se "sapply" tab <- table(tridy) emp.cetnosti <- tab[] # empiricke cetnosti jednotlivych poctu teor.cetnosti <- n*c(ppois(8,lambda),dpois(9:21,lambda),1-ppois(21,lambda)) # teoreticke cetnosti chi2 <- sum((emp.cetnosti-teor.cetnosti)^2/teor.cetnosti) df <- length(teor.cetnosti)-2 p.hodnota <- 1-pchisq(chi2,df) krit.hodnota <- qchisq(0.95,df) # # Pozn.: volba trid neni jednoznacna a test je citlivy na tuto volbu. # Existuji ruzna kriteria a dopuruceni, kdy je aproximace limitnim chi^2 # rozdelenim pouzitelna. Napriklad pokud vsechny teoreticke cetnosti jsou # vetsi nez 5q, kde q je podil trid, pro nez teoreticke cetnosti jsou mensi # nez 5 (to je v nasem pripade splneno). # # Lze take vyuzit funkci goodfit() v knihovne 'vcd': library(vcd) # Pokud nemame knihovnu nainstalovanu, lze to udelat pres menu volbou # Packages - Install Package(s) gf <- goodfit(data,type="poisson",method="MinChisq") # volba method="MinChisq" znamena, ze se pouzije odhad intenzity, ktery minimalizuje # prislusnou chi^2-statistiku ("nejlepsi shoda mezi teoretickymi a empirickymi cetnostmi") summary(gf) plot(gf) # odmocnina z Poissona priblizne jako normalni # # Tento postup vsak nepouziva deleni na tridy (pracuje se vsemi pozorovanymi pocty), # proto je pouziti aproximace chi^2 rozdelenim problematicke. # # ########################################### ######### Korelacni koeficient ########## ########################################### # # Nejdrive natahneme potrebnou knihovnu: # library(MASS) # pro praci s vicerozmernymi nah. velicinami # sigma1 <- diag(2) # jednotkova matice 2x2 pro nekorelovane veliciny sigma2 <- cbind(c(1,.7),c(.7,1)) # variancni matice s korelaci mu <- c(0,0) # vektor strednich hodnot # # Nagenerujeme dvoje data z dvourozmerneho normalniho rozdeleni: # (pro funkci mvrnorm() bylo nutne nacist knihovnu 'MASS') # data1 <- mvrnorm(n=250, mu=mu, Sigma=sigma1) data2 <- mvrnorm(n=250, mu=mu, Sigma=sigma2) # # A podivame se, jak vypadaji korelace mezi slozkami: # cor(data1) # vrati celou korelacni matici cor(data1[,1],data1[,2]) # vrati jen korelacni koeficient X a Y, tj. zde mimodiagonalni prvek # cor(data2) # # Jakpak ono to vypada, kdyz se to nakresli: # plot(data1) points(data2, pch='+', col="red") # # Jak by to vypadalo, kdyby korelace byla zaporna, rekneme -0.8? # # Test(y) o korelacnim koeficientu: # cor.test(data1[,1],data1[,2]) # musime zadat 2 vektory pozorovani, tedy oba sloupce nasich dat cor.test(data2[,1],data2[,2]) # # opet se daji volit parametry jako alternativni hypoteza, hladina vyznamnosti...) cor.test(data2[,1],data2[,2], alternative="less", conf.level=.9) # # Krome klasickeho vyberoveho korelacniho koeficientu (Pearsonova) je # moznost uzit i Spearmanuv a Kendalluv, to se provede volbou parametru # method v 'cor' nebo 'cor.test'. cor(data1,method="spearman") cor(data2,method="kendall") # q() # Koncime ...