Rko.txt

#######################
########################
## Pravdepodobnost a statistika NMAI059
## Pouziti programu R

## http://artax.karlin.mff.cuni.cz/~novap4dm/Uvod.r

# spustime pres Rgui
# programy obvykle piseme a ukladame jako scripty .r
# radky programu ze scriptu spoustime pomoci ctrl+r
# volane prikazy a vysledky se objevuji v konzoli

#######################
## R funguje jako kalkulacka

3*8^2+10/log(5)+exp(-20)+factorial(6)

#######################
## ALOKOVANI A VOLANI PROMENNYCH

# !pozor, vsechny nazvy promennych i funkci jsou CASE SENSITIVE! 

# vektor o delce 1 obsahujici hodnotu 3
x=3				 # lze i x<-3
# vektor obsahujici 3 cisla 
y=c(1,5,9)
# spojeni x a y do vektoru delky 4
z=c(x,y)
#vypsani obsahu jednotlivych promennych
x
y
z
# vektor obsahujici desetkrat nulu
x=rep(0,10)
x
#vektor obsahujici cisla 1 az 10
x=1:10
x
#zavolani ctvrteho prvku vektoru
x[4]
#zavolani ctvrteho a osmeho prvku vektoru
x[c(4,8)]

#matice 2x3 obsahujici cisla 1,2,3,4,5,6
X=matrix(1:6,2,3)
X
X[1,2]

# logicke promenne
a=T			# jde i a=TRUE
b=F			# jde i a=False
			# s logickymi promennymi jde pracovat i jako s nulami a jednickami, tj. jimi nasobit, scitat je apod.
			# and a or pomoci & a | 
a*a

# stringove a faktorove promenne
z=c(rep("M",5),rep("F",5))
z

class(z)		## trida objektu
z=as.factor(z)	## prevedeme na faktor (promenna s nekolika moznymi urovnemi)
z
class(z)
class(x)
class(a)

# data frame - tabulka dat
x=1:10;y=rep(0,10);

dat=data.frame(x,y,z)
dat
names(dat)=c("poradi","plat","pohlavi")
dat
dat[1,]	
dat$pohlavi	
dat$pohlavi[1:6]

#######################
## POUZIVANI VESTAVENYCH VYPOCTOVYCH FUNKCI

# v argumentu je vektor cisel
sum(x)		#soucet
mean(x)		#prumer
var(x)		#S^2
sd(x)			#sqrt(S^2)
length(x)		#vrati delku vektoru x


# kombinovani vektoru a cisla
x^2		#kazdy prvek umocni
log(x)	#kazdy prvek zlogaritmuje
x-5		#od kazdeho prvku odecte 5
x-mean(x)   #od kazdeho prvku odecte prumer x
x<=5		#vrati vektor stejne dlouhy jako x, obsahujici TRUE nebo FALSE

sum(x<=5)	#vrati pocet pripadu kdy to nastalo

#priklad vypocet S^2 (stejne cislo dostanem volanim var(x) ):
sum( (x-mean(x))^2 )/9

## CYKLY

for(i in 1:10){
print(x[i])
}

for(i in 10:1){
print(x[i])
}

z=rep(0,10)
for(i in 1:10){z[i]=mean(x[-i])}
z

i=1;z=0
while(z<=20){
z=z+log(i)
i=i+1
}
z

##PODMINKY

if(x[1]==6){z=4}else{z=6}	#  < > <= >= !=

for(i in 1:10){
 if((x[i]<=5)&sum(x[-i]<=20)){
  print("Dobry den")
 }
}

####################
## STATISTICKE FUNKCE

t.test(x,mu=5,alternative="less")		##t-test o stredni hodnote rozdeleni
t.test(x,mu=1)					##default je proti oboustranne alternative mu<>mu0

# hypotezu ze opravdova stredni hodnota je mu zamitame na hladine 0.05, pokud p-value je mensi nez 0.05
# zaroven dostavame i intervalovy odhad stredni hodnoty

# intervalovy odhad dostaneme i jako
n=length(x)

c(
mean(x)-qt(0.975,df=n-1)*sd(x)/sqrt(n),
mean(x)+qt(0.975,df=n-1)*sd(x)/sqrt(n)
)

#qt(0.975,df=n-1) je kvantil (neboli inverze distribucni funkce) t-rozdeleni o n-1 stupnich volnosti v bode 1-alpha/2 pro alpha=0.05

#priblizny interval spolehlivosti dostaneme, pokud misto t rozdeleni pouzijeme CLV a kvantily normalniho rozdeleni, tj phi^-1:
c(mean(x)-qnorm(0.975)*sd(x)/sqrt(n),mean(x)+qnorm(0.975)*sd(x)/sqrt(n))

## CHARAKTERISTIKY ROZDELENI

?qnorm
## lze volat hodnotu hustoty/pravdepodobnosti, distribucni funkce, kvantily a nahodna cisla
?qexp()

#lze pocitat i kvantily pro diskretni rozdeleni
# q(u)=inf(x:F(x)>=u)

qpois(0.95,lambda=5)

ppois(8,lambda=5)
ppois(9,lambda=5)
# jinak receno, pravdepodobnost ze velicina z poissonova rozdeleni s lambda=5 bude <=9 je vetsi nez 95%
#	        a pravdepodobnost ze velicina z poissonova rozdeleni s lambda=5 bude <=8 je mensi nez 95%
		 
####################
## SIMULACE

## nastaveni pocatku simulace (vhodne aby slo simulaci zopakovat)
set.seed(12345678)

#vygenerujeme 50 pozorovani z normalniho rozdeleni N(100,15^2)
n=50
x=rnorm(n,mean=100,sd=15)		# lze i rnorm(n,100,15^2)
x

#totez, ale z exponencialniho s lambda=5
x=rexp(n,rate=5)

## Priklad 1: chceme simulacemi najit takove cislo, ze na 95% jej soucet padesati nezavislych velicin z Exp(5) neprekroci.
## Tj. horni intervalovy odhad
## Vyrobime N krat padesatici cisel, zapamatujeme si vzdy jejich soucet.
## Vezmeme 95%-vyberovy kvantil, tj. cisla seradime a podivame se na to, ktere cislo je v 95%-poradi

N=1000
z=rep(0,N)
for(i in 1:N){
 z[i]=sum(rexp(50,5))
}
quantile(z,0.95)
sort(z)[0.95*N]

## vyslo 12.5, pro jiny seed by to mohlo byt jine

## Priklad 2: ve stejne situaci chceme najit pravdopodobnost, ze soucet takovych padesati cisel nebude vetsi nez 13
## Jednoduse spocitame, v kolika pripadech z N toto nastalo.
## Vyrobime vektor TRUE/FALSE hodnot z<=13, s tim se da pracovat jako s jednickami a nulami, ktere zprumerujeme. 
## Podivame se, v kolika procentech pripadu se tak stalo.

mean(z<=13)

mean(z<=8.5)

## Priklad 3: chceme zjistit, zda dvouvyberovy t-test zamitne nebo nezamitne hypotezu o shode strednich hodnot
## pokud generujeme jeden vyber z N(0,5^2) a druhy z N(1,5^2).
## zkusime to pro N=20 a N=1000

N=20
x=rnorm(N,0,5)
y=rnorm(N,1,5)
t.test(x,y)

N=1000
x=rnorm(N,0,5)
y=rnorm(N,1,5)
t.test(x,y)

## vidime, ze zde pro N=20 test hypotezu nezamitl, pro N=1000 uz jo.
## pro lepsi nahled fungovani t-testu bychom mohli zopakovat simulaci mnohokrat

####################
## FUNKCE

rozdil=function(x,y){
 rozdil=x-y;return(rozdil)
}
rozdil(5,2)
z=rozdil(5,2)

##
x=rnorm(50)
prum1doN=function(N){
 prum1doN=ifelse(N<=length(x),mean(x[1:N]),mean(x))
 return(prum1doN)
}
prum1doN(6)


####################
## GRAFY

#frekvencni diagram (histogram)
n=100
x=rnorm(n,100,15)
hist(x)

## history->recording zapne zaznamenavani grafu

#graf zavislosti dvou promennych
x=1:20
y=rnorm(20,x,1)
plot(y~x)
plot(y~x,type="l")
points(y~x,pch="+")


# qqnorm - slouzi k porovnani s normalnim rozdelenim
# melo by byt rovne okolo osy kvadrantu
x=rexp(20)
y=rnorm(20)
qqnorm(x)
qqline(x)
qqnorm(y)
qqline(y)