#Zadatak 1:
#  Prikupljeni su podaci o visini ljudi:
#  159 188 175 176 177 168 162 188
#183 187 187 162 184 161 180 169
#195 171 170 199 181 169 189 191
#172 182 183 178 180 165 185 202
#183 187 188 182 163 179 178 188
#Ucitajte podatke koristeći naredbu scan, 
#izracunajte karakteristicnu petorku te procijenite 95% p.i. 
#za ocekivanu visinu covjeka.
visine=scan("visine.txt")
#karakteristicna petorka
quantile(visine)
n=length(visine)#mozemo primijeniti CGT
mu_hat=mean(visine)
sigma_hat=sd(visine)
alpha=0.05
z=qnorm(alpha/2, mean=0, sd=1, lower.tail=F)
#granice procijenjenog 95%-tnog pouzdanog intervala
pi_lwr=mu_hat-z*sigma_hat/sqrt(n)
pi_upr=mu_hat+z*sigma_hat/sqrt(n)
pi_lwr# 175.8669
pi_upr# 182.4331

#da nismo koristili CGT vec Studentovu t.distribuciju
pi_lwr=mu_hat-qt(alpha/2, df=n-1)*sigma_hat/sqrt(n)
pi_upr=mu_hat+qt(alpha/2, df=n-1)*sigma_hat/sqrt(n)
pi_lwr# 175.7619
pi_upr# 182.5381

#2.nacin:
t=t.test(visine, conf.level=0.95)
t$conf.int#175.7619 182.5381


#Zadatak 2:
#Ucitajte podatke iz .csv datoteke koristeci read.csv o broju nukleotida u ljudskom genomu
#"04e1HumanGeneLengths.csv". Koliki je broj podataka u toj bazi? 
#Nacrtajte histogram tih podataka te usporedite vrijednosti medijana i artimeticke sredine uzorka. 
#Odredite 99% p.i. za prosječan broj nukleotida u ljudskom genomu. 
#Odredite i 99% p.i. za varijancu koristeći naredbu VarCI iz paketa "DescTools". 
#(Koristite "install.packages("DescTools")" i "library(DescTools)").


#ucitavanje: ImportDataset>FromText(readr)>Browse, preimenovati u gene, npr.
data<-gene$gene.length#gene je data.frame, a s gene$gene.length dobivamo podatke iz stupca s imenom gene.length
#broj podataka
n=length(data)
n#20290
barplot(data)#nepregledno - podatke treba podijeliti u razrede, no najprije pogledajmo
#sto nam vraca iduca naredba
h<-hist(data, prob=T)
h$breaks#[0,5000], <5000, 10000], <1000, 15000], <15000, 20000], <20000, 25000], ...
h$counts#vidimo da ima nekoliko izuzetno velikih vrijednosti (outliera) - pogledajmo boxplot
boxplot(data)
#prosjek i medijan
mean(data)#2622.027 - prosjek je veci od medijana zbog tih iznimno velikih vrijednosti
median(data)#2226.5


#iz h$breaks i h$counts vidimo da imamo svega 6 vrijednosti vecih od 20000
#pokusajmo te vrijednosti izbaciti iz podataka
data1<-data[data<20000]
hist(data1, prob=T)#ljepše se vidi početak histograma
mean(data1)#2611.646
median(data1)#2226

#99% p.i. za prosjecan broj nukleotida
#zbog velikog broja podataka mozemo koristiti CGT
mu_hat=mean(data)
sigma_hat=sd(data)
z=qnorm(0.01/2, 0, 1, lower.tail=F)
pi_lwr=mu_hat-z*sigma_hat/sqrt(n)
pi_upr=mu_hat+z*sigma_hat/sqrt(n)
pi_lwr#2585.192
pi_upr#2658.862

#alternativno:
t=t.test(data, conf.level=0.99)
t$conf.int#2585.189 2658.866


#99% p.i. za varijancu
install.packages("DescTools") 
library(DescTools)
VarCI(data, conf.level=0.99)#4045038 4257324 



#Zadatak 3: duljina vremenskog perioda prije nego se dogodi dioba bakterije 
#eksponencijalno je distribuirana s parametrom \lambda za koju vrijedi. 
#a) Procijenite 90% p.i. za nepoznati parametar \lambda na temelju podataka:
#  0.3600 0.3427 0.4308 0.4177 0.1849 0.1004 0.4310 0.3315 0.1714 0.6693 0.1407 0.7263 1.0726 0.1859
#0.1982 0.3258 0.0700 0.1031 0.3686 0.2581 0.0299 0.3694 0.0824 0.5240 1.6109 0.1437 0.9101 0.3789
#0.2711 0.2790 0.5949 0.7708 0.9700 0.0952 0.1296 0.0174 0.1173 0.5217 0.2715 0.9197

#b) Procijenite 90% p.i. za \lambda na podatcima koje sami generirate. 
#Duljina uzorka neka bude 40, i neka podaci dolaze iz Exp(1) razdiobe. 
#Koristite funkciju rexp. Je li taj pravi paramtetar \lambda=1 unutar vašeg p.i.?

#a)
#ucitavanje podataka
dioba=scan("dioba.txt")
#za parametar ocekivanja mu procjene za pouzdane intervale dobivamo s
#za male uzorke
#[pi_lwr, pi_upr]=[mean(podaci)-t*sd(podaci)/sqrt(n), mean(podaci)+t*sd(podaci)/sqrt(n)], t=qt(alpha/2, df=n-1, lower.tail=F)
#za velike uzorke: uz CGT
#[pi_lwr, pi_upr]=[mean(podaci)-z*sd(podaci)/sqrt(n), mean(podaci)+z*sd(podaci)/sqrt(n)], z=qnprm(alpha/2, 0, 1,lower.tail=F)
#ovdje nam je  EX=1/\lambda=mu,x->1/x je padajuca funkcija pa vrijedi
#pa izvedemo pi_lwr_za_lambda=1/pi_upr_za_mu=1/(mean(podaci)+t*sd(podaci)/sqrt(n)) 
#pi_upr_za_lambda=1/pi_lwr_za_mu=1/(mean(podaci)+t*sd(podaci)/sqrt(n)) (primjer za male uzorke, slicno za velike)

#imamo 40 podataka: ako  koristimo CGT
n=length(dioba)
alpha=0.1
z=qnorm(alpha/2, 0, 1, lower.tail=F)
pi_lwr_za_mu=mean(dioba)-z*sd(dioba)/sqrt(n)
pi_upr_za_mu=mean(dioba)+z*sd(dioba)/sqrt(n)
pi_lwr_za_lambda=1/pi_upr_za_mu
pi_lwr_za_lambda#2.059449
pi_upr_za_lambda=1/pi_lwr_za_mu
pi_upr_za_lambda#3.233544

#alternativno, ako ne koristimo CGT, nego Studentovu t-distribuciju:
t=t.test(dioba, conf.level=0.90)
1/t$conf.int# 2.050392 3.256127

#b)
#parametar lambda se kod eksponencijalne distribucije u R-u naziva rate
#generiramo parametre iz eksponencijalne razdiobe (R-uvod) funkcijom rexp(n, rate), 
#gdje je n broj podataka koje zelimo generirati, a rate parametar lambda
dioba1<-rexp(40, rate=1)
#pravi parametar lambda=1 (tako smo zadali kod generiranja)
#ako rješavamo uz CGT:
alpha=0.1
n=40
z=qnorm(alpha/2, 0, 1, lower.tail=F)
pi_lwr_za_mu=mean(dioba1)-z*sd(dioba1)/sqrt(n)
pi_upr_za_mu=mean(dioba1)+z*sd(dioba1)/sqrt(n)
pi_lwr_za_lambda=1/pi_upr_za_mu
pi_lwr_za_lambda# 0.7621016
pi_upr_za_lambda=1/pi_lwr_za_mu
pi_upr_za_lambda#1.191159
#da smo rjesavali s t-distribucijom
#ne bi bilo dobro da smo samo napisali iduce dvije linije, jer nismo invertirali granice
#t1=t.test(dioba1, conf.level=0.9)
#t1$conf.int


#ako znamo pravu vrijednost parametra, vazna provjera kod procjena pouzdanih intervala je da 
#oni (skoro uvijek - ovisno o vrijednosti alpha) sadrže tu pravu vrijednost,
#a sigurno sadrže procjenu za vrijednost tog parametra na temelju uzorka


#Zadatak 4:
#Ucitajte podatke iris te izracunajte 90% p.i. za sepal length. 
#Izracunajte 90% p.i. za standardnu devijaciju za sepal length, koristite naredbu VarCI.

#ucitavanje podataka koji su vec poznati u R-u - to su obicno data.frameovi
data(iris)#iris je data.frame
podaci=iris$Sepal.Length

n=length(podaci)
n#150
#90%-tni p.i. za ocekivanje
#mozemo koristiti CGT
alpha=0.1
z=qnorm(alpha/2, 0, 1, lower.tail = F)
mu_hat=mean(podaci)
sigma_hat=sd(podaci)
pi_lwr=mu_hat-z*sigma_hat/sqrt(n)
pi_lwr#5.732123
pi_upr=mu_hat+z*sigma_hat/sqrt(n)
pi_upr#5.954544

#alternativno (ako smo koristili t-distribuciju)
t=t.test(podaci, conf.level=.9)
t$conf.int#5.731427 5.955240


#90%-tni p.i. za stnd. dev.
#znamo rijesiti b) zadatak za parmetar varijance
#stnd. dev.=sqrt(varijanca) - tu cemo transformaciju primijeniti na granice p. intervala za varijancu
#x->sqrt(x) je rastuca funkcija, pa vrijedi
#pi_lwr_za_sd=sqrt(pi_lwr_za_var)
#pi_upr_za_sd=sqrt(pi_upr_za_var)
var_ci<-VarCI(podaci, conf.level=0.9)#0.5724186 0.8389097 
pi_lwr=sqrt(0.5724186)
pi_lwr#0.7565835
pi_upr=sqrt(0.8389097)
pi_upr#0.9159201
sd(podaci)# 0.8280661 - procijenjeni pouzdani interval mora sadrzavati procjenu odgovarajuceg parametra


#DZ: Trajanje neke kemijske reakcije mjereno je u sekundama. 
#Na uzorku od 150 mjerenja imamo sljedece
#informacije: x¯ = 0.78, s = 0.05. Procijenite 99% pouzdani interval za prosjecnu duljinu trajanje te
#kemijske reakcije.
n=150#mozemo koristiti CGT
mu_hat=0.78
sigma_hat=0.05
alpha=0.01
z=qnorm(alpha/2, 0, 1, lower.tail = F)
pi_lwr=mu_hat-z*sigma_hat/sqrt(n)
pi_lwr#0.7694842
pi_upr=mu_hat+z*sigma_hat/sqrt(n)
pi_upr#0.7905158

#DZ: Mjeren je gestacijski period kod macaka. Dobiveni su sljedeci podatci:
#  63, 65, 69, 66, 63, 67, 67, 66, 64, 64, 64, 65, 63, 64, 65, 64, 64, 63, 63, 65, 67, 62, 63, 66, 64,
#64, 63, 63, 64, 65, 66, 64, 63, 64, 64, 64, 68, 65, 63, 62, 64, 64, 64, 65, 66, 65, 64, 63, 66, 64.
#a) Procijenite 90%, 95% i 99% pouzdani interval za prosjecno trajanje gestacijskog perioda.
#b) Procijenite 90%, 95% i 99% pouzdani interval za varijancu trajanja gestacijskog perioda. (Uputa:
#koristite funkciju VarCI() iz paketa DecsTools)
#c) Napisite interpretaciju 90% pouzdanih intervala za prosjecno trajanje i varijancu gestacijskog
#perioda.

#ucitamo podatke u vektor
macke=c(63, 65, 69, 66, 63, 67, 67, 66, 64, 64, 64, 65, 63, 64, 65, 64, 64, 63, 63, 65, 67, 62, 63, 66, 64, 64, 63, 63, 64, 65, 66, 64, 63, 64, 64, 64, 68, 65, 63, 62, 64, 64, 64, 65, 66, 65, 64, 63, 66, 64)
n=length(macke)
n#50

#p.i. za parametar ocekivanja
#mozemo s CGT
alpha=0.1#0.05, 0.01
z=qnorm(alpha/2, 0,1, lower.tail=F)
mu_hat=mean(macke)  
sigma_hat=sd(macke)
pi_lwr=mu_hat-z*sigma_hat/sqrt(n)
pi_lwr#64.11398
pi_upr=mu_hat+z*sigma_hat/sqrt(n)
pi_upr#64.80602
#alternativno: uz t-distribuciju
t=t.test(macke, conf.level = 0.9)#0.95, 0.99
t$conf.int

#p.i. za parametar varijance
VarCI(macke, conf.level = 0.9)#0.95, 0.99
#1.634341 3.195374 

#interpretacija npr. 90%-tnog p.i. za parametar ocekivanja:
#vjerojatnost da 90%-tni p.i. za param. ocekivanja sadrzi pravu vrijednost tog parametra je 90%
#ili: ako 100 puta određujemo 90%-tni p.i. za param. ocekivanja, bar 90 njih ce sadrzavati
#pravu (nepoznatu) vrijednost tog parametra

#slicno za varijancu


#DZ: Simulirajte podatke iz eksponencijalne razdiobe s parametrom 2. Na temelju uzorka duljine 200
#procijenite 95% pouzdani interval za ocekivanje te distribucije.
data=rexp(200, rate=2)
n=200#moze CGT
alpha=0.05
z=qnorm(alpha/2, 0, 1, lower.tail=F)
mu_hat=mean(data)
sigma_hat=sd(data)
pi_lwr_mu=mu_hat-z*sigma_hat/sqrt(n)
pi_upr_mu=mu_hat+z*sigma_hat/sqrt(n)
pi_upr_lambda=1/pi_lwr_mu
pi_lwr_lambda=1/pi_upr_mu
pi_lwr_lambda#1.932312
pi_upr_lambda#2.503364

#dobro je da smo dobili da ovaj p.i. sadrzi pravu vrijednost za lambda (2)

#DZ: Simulirajte podatke iz eksponencijalne razdiobe s parametrom 5. Na temelju uzorka duljine 500
#procijenite 99% pouzdani interval za ocekivanje te distribucije. - slicno, drugi n i alpha


#DZ: Simulirajte podatke iz geometrijske razdiobe s parametrom p = 0.4. Na temelju uzorka duljine 80
#procijenite 90% pouzdani interval za ocekivanje te distribucije.
data=rgeom(80, prob=0.4)
data
n=80#moze CGT
#za X iz G(p) vrijedi EX=1/p=mu - okrenuti granice p.i. za mu
alpha=0.01
z=qnorm(alpha/2, 0,1, lower.tail=F)
mu_hat=mean(data)  
sigma_hat=sd(data)
pi_lwr_mu=mu_hat-z*sigma_hat/sqrt(n)#1.095699
pi_upr_mu=mu_hat+z*sigma_hat/sqrt(n)#2.329301
pi_upr_p=1/pi_lwr_mu
pi_lwr_p=1/pi_upr_mu
pi_lwr_p#0.4645664
pi_upr_p#0.7583641
1/mean(data)#0.5839416 - procjena za parametar p na temelju uzorka
#moze se dogoditi da p.i. ne sadrzi pravu vrijednost parametra, ali mora sadrzavati procjenu za taj parametar