######################################################### ######################################################### ### Cviceni NSTP097 9.12.2008 ### Uvod do R, pravdepodobnostni rozdeleni atd. ######################################################## ######################################################### ######################## ### Ukol 1 ####################### ################ # UKOL 1 a) # chceme si nakreslit hustotu exponencialniho rozdeleni Exp(4) # ulozim si do parametru lam hodnotu 4 a dal budeme pracovat s timto parametrem # kdybyste chteli vsecko prepocitat pro jinou hodnotu lambda -> staci zmenit zde a pak znovu projet kod lam=4 xx=seq(0,2,length=500) # v techto bodech budeme hustotu pocitat hust=dexp(xx,rate=lam) # hodnoty hustoty, parametr rate je lambda v hustote (jak to znate) plot(xx,hust,type="l") # vykresli # zkuste zmenit nastaveni type, co se stane.. type="p", type="b", type="s" atd # kresleni obrazku lze plot(xx,hust,type="l") nebo (je to totez) plot(hust~xx,type="l") # NEBO to jde cele jen jednim prikazem: curve(dexp(x,rate=lam),0,2) # curve vykresli zadanou krivku (nejaka funkce x) od 0 do 2 # necham na Vas, co se Vam libi vic.... # v obrazku je mozne volit ruzne parametry: barva, tloustka cary, typ cary, popisky, velikost popisku atd... plot(xx,hust,type="l",col="red",lwd=3,ylab="f(x)",xlab="x",cex.lab=1.5,main=paste("Hustota Exp(",lam,")")) ############################################# ## VSUVKA: ## ?? Jak ulozit obrazek, kdyz uz ho mame ?? ####### ## a) Pokud chceme obrazek potom pouzit do neceho jako word atd, tak se nam asi hodi format bmp nebo jpeg ## b) Pro LaTeX chceme postscript nebo pdf ## v podstate existuji dva zpusoby: klikaci a prikazem #### Jednoduchy (ale nekdy hloupy) pristup: kliknout na Export a vybrat zpusob, jakym chceme obrazek ulozit #### Sofisikovanejsi (a flexibilnejsi) postup #### napr. pro format pdf pdf("hustota.pdf", width=7, height=10) # vytvori soubor "hustota.pdf" v pracovnim adresari; kresli do nej plot(xx,hust,type="l",col="red",lwd=3,ylab="f(x)",xlab="x",cex.lab=1.5,main=paste("Hustota Exp(",lam,")")) dev.off() # uzavre soubor; odted kresli zase na obrazovku # analogicky postup pro jine formaty; pozor na jina nastaveni; napr. # jpeg("bin10.jpg", width=1280, height=1024) # postscript("bin10.jpg", width=7, height=10, horizontal=FALSE) atd #################################################### # UKOL 1 b) # nyni chceme distribucni funkci Exp(lam) dfce=pexp(xx,rate=lam) # spocitame ji v bodech xx plot(dfce~xx,type="l") # nakresli # pripomente si vlastnosti distribucni funkce.... # nebo to jde zase prikazem curve(pexp(x,rate=lam),0,2) curve(pexp(x,rate=lam),-3,5) # co kdyz zmenime min a max v obrazku... # je vsem jasne, proc je distr. fce nulova pro x<0? ################################## # UKOL 1 c) # generovani nahodneho vyberu z exp. rozdeleni # chceme napr. 25 velicin -> ulozime si zase 25 do nejake promenne n (kdybysme chteli casem vic, bude jednoduche to zmenit...) lam=4 n=25 nvyb=rexp(n,rate=lam) # vygeneruje n hodnot z exp. rozdeleni s parametrem lambda=lam (v nasem pripade 25 cisel z Exp(4) summary(nvyb) # zakladni popisne statistiky pro dany vyber mean(nvyb) # chceme-li pouze pr?m?r # vsimnete si hlavne prumeru. vime, ze prumer je jakousi "analogii" (resp. odhadem) stredni hodnoty... # jaka je stredni hodnota v Exp(lambda)? hist(nvyb) # histogram - pripomina nam hustotu, co jsme videli? # spolecne s hustotou: hist(nvyb,prob=T) curve(dexp(x,rate=lam),0,2,add=T) ##################3 # poznamka na okraj pro zajemce: kdyz chcete pripadne Vasi simulaci opakovat, nastavte pred spustenim parametr seed # napr. set.seed(2008) -> potom vzdycky dostanete totez (hodi se u bakalarek, diplomek atd) ##################### ## UKOL 1 d) # Dalsi rozdeleni: ## ZKUSTE SAMI!!! -> hustoty, distribucni funkce a nahodny vyber ##################################################################### ## Nejake moje obrazky: # Ruzne gamma rozdeleni curve(dgamma(x,rate=2,shape=0.5),0,3,lwd=2,xlab="x",ylab="f(x)",cex.lab=1.5) curve(dgamma(x,rate=0.5,shape=0.1),0,3,lwd=2,add=T,col="red") curve(dgamma(x,rate=2,shape=2),0,3,lwd=2,add=T,col="blue") curve(dgamma(x,rate=1,shape=2),0,3,lwd=2,add=T,col="green") legend("topright",lty=c(1,1,1,1),col=c("black","red","blue","green"),legend=expression(Gamma(2,0.5),Gamma(0.5,0.1),Gamma(2,2),Gamma(1,2)),inset=0.05) ######## ## Ruzne normalni rozdeleni curve(dnorm(x,0,1),-5,5,lwd=2) curve(dnorm(x,0,2),-5,5,lwd=2,col="red",add=T) curve(dnorm(x,0,3),-5,5,lwd=2,col="blue",add=T) curve(dnorm(x,-2,2),-5,5,lwd=2,add=T,col="green") legend("topright",lty=c(1,1,1,1),col=c("black","red","blue","green"),legend=expression(N(0,1),N(0,4),N(0,9),N(-2,4)),inset=0.05) ################################################################## ################### ## Priklad 1 e) ################## # cauchyho rozdeleni n=25 vyber.cauchy=rcauchy(n) summary(vyber.cauchy) hist(vyber.cauchy) n=100 vyber.cauchy=rcauchy(n) summary(vyber.cauchy) hist(vyber.cauchy) # vidite tezke chvosty? # vsimnete si rozdilu mezi medianem a prumerem! ################# ## Priklad 1 f) ################ # N(0,1) a C(0,1) do jednoho obrazku: curve(dnorm(x,0,1),-5,5) curve(dcauchy(x,0,1),-5,5,add=T,col="red") legend("topleft",lty=c(1,1),col=c("red","black"),legend=c("C(0,1)","N(0,1)"),inset=0.05) # uplne totez lze provest i takto (vice prikazu...): x=seq(-5,5,length=500) hustc=dcauchy(x,0,1) hustn=dnorm(x,0,1) plot(x,hustc,type="l",ylim=1.2*range(c(hustc,hustn)),ylab="Hustota") lines(x,hustn,lty=2) legend("topleft",lty=c(1,2),legend=c("C(0,1)","N(0,1)"),inset=0.05) ################## ## Priklad 1 g) ################# # ZKUSTE SAMI! # jaka rozdeleni vlastne chceme porovnat? # (vzpomente si, jak vypada stredni hodnota u exp a gamma rozdeleni...) ######################################################### # UKOL 2: ######################################################### # a) Chceme Gamma(0.5,2) curve(dgamma(x,rate=0.5,shape=2),0,10,lwd=2) # b) bud muzeme pouzit for cyklus - coz ale v Rku neni uplne vhodne - # nebo to udelat pres matice (coz udelame) # usporadame vybery do matice o 25 sloupcich a 50 radcich n=25 nvyb.50=rgamma(50*n,rate=0.5,shape=2) nvyb.50=matrix(nvyb.50,ncol=n) prumery=apply(nvyb.50,1,mean) # nyni spocteme prumer pres radky prumery=rowMeans(nvyb.50) # nebo takto (dela to totez) var(prumery) ## odhad rozptylu: hist(prumery) ############# ## ZOPAKUJTE -SAMI # histogramy kreslete takto: hist(prumery, main=paste("Pocet pozorovani n =",n),xlim=c(2.5,5.5)) ###################################### ## nejake moje obrazky moje.funkce= function(n,opak=50){ vyber=rgamma(opak*n,rate=0.5,shape=2) vyber=matrix(vyber,ncol=n) prumery=apply(vyber,1,mean) print(var(prumery)) hist(prumery, main=paste("Pocet pozorovani n =",n),xlim=c(2.5,5.5)) } par(mfrow=c(2,2)) moje.funkce(25) moje.funkce(250) moje.funkce(2500) moje.funkce(25000) ##################### ## c) ## cauchy.funkce= function(n,opak=50){ vyber=rcauchy(opak*n) vyber=matrix(vyber,ncol=n) prumery=apply(vyber,1,mean) print(var(prumery)) hist(prumery,main=paste("Pocet pozorovani n =",n)) } par(mfrow=c(2,2)) cauchy.funkce(25) cauchy.funkce(250) cauchy.funkce(2500) cauchy.funkce(25000) ########################################## ### Uloha 3 ######################################### # zkuste sami nebo si spustte nasledujici funkci moje.funkce2= function(n,opak=100){ vyber=rgamma(opak*n,rate=0.5,shape=2) vyber=matrix(vyber,ncol=n) prumery=apply(vyber,1,mean) z=sqrt(n)*(prumery-4)/sqrt(8) hist(z,prob=T,main=paste("Pocet pozorovani n =",n)) curve(dnorm(x),add=T,lwd=2,col="red") } par(mfrow=c(2,2)) moje.funkce2(25) moje.funkce2(250) moje.funkce2(2500) moje.funkce2(25000) ### c) cauchy.funkce2= function(n,opak=100){ vyber=rcauchy(opak*n) vyber=matrix(vyber,ncol=n) prumery=apply(vyber,1,mean) z=(prumery-mean(prumery))/sd(prumery) hist(z,prob=T,main=paste("Pocet pozorovani n =",n)) curve(dnorm(x),add=T,lwd=2,col="red") } par(mfrow=c(2,2)) cauchy.funkce2(25) cauchy.funkce2(250) cauchy.funkce2(2500) cauchy.funkce2(25000)