# Intervali poverenja

# Zelimo da nadjemo takve U_n, i V_n da nepoznati parametar Teta 
# upadne u interval [U_n, V_n] sa verovatnocom
# beta (unapred zadatom) tj:
# P{U_n <= Teta <= V_n}=beta
# beta nazivamo nivo poverenja.

# (1)
# Naci 95%-tni interval poverenja za parametar "m", ako je 
# nasa slucajna velicina X~N(m, sigma^2). 
# Poznato nam je sigma^2 i n - broj elemenata u
# prostom slucajnom uzorku.

#
n=64
x_sr=118
beta=0.95
sigma_kv= 100 #- poznata vrednost


# Resenje:
c=qnorm((1+beta)/2, 0, 1) # isto sto i qnorm((1+beta)/2) isto sto i qnorm(0.975)
c


interval_1=x_sr-c*sqrt(sigma_kv/n)
interval_2=x_sr+c*sqrt(sigma_kv/n)
interval_1
interval_2
(interval_poverenja=c(interval_1, interval_2))


# (2)
# Naci 95%-tni interval poverenja za m na osnovu uzorka obima n
# ako je X~N(m, sigma^2), sigma^2 - nepoznato.
# ako je
n=25
x_sr=2.6
ocena_sigma_kv=170.36
beta=0.95


# Resenje:
c=qt((1+beta)/2, n-1) #qt - kvantili studentove raspodele, drugi argument predstavlja broj stepeni slobode.
c
(interval_poverenja=c(x_sr-c*sqrt(ocena_sigma_kv/n), x_sr+c*sqrt(ocena_sigma_kv/n) ))


# Da primenimo ovakav interval poverenja i na prvi zadatak. Moramo da dobijemo nesto siri 
n=64
x_sr=118
beta=0.95
sigma_kv= 100 # ovde sada smatramo da je to ocenjena vrednost

c1=qt((1+beta)/2, n-1) 
c1


interval_1=x_sr-c1*sqrt(sigma_kv/n)
interval_2=x_sr+c1*sqrt(sigma_kv/n)
(interval_poverenja=c(interval_1, interval_2))
# I zaista vidimo da je ova malcice siri interval poverenja.


# INTERVALI POVERENJA

# U prva dva zadatka smo pominjali intervale poverenja u slucaju normalne raspodele 
# Norm(m, sigma^2) za parametar "m". 
# I videli smo da se taj slucaj deli
# na dva podslucaju u zavisnosti od toga da li znamo sigma^2 
# ili je ocenjujemo za dati uzorak.
# U slucaju kada znamo sigma^2 trazimo interval poverenja pomocu statistike:
# Z=(X_sr-m)*sqrt(n)/(sigma),
# gde je
# X_sr - srednja vrednost uzorka za dato obelezje.
# m - parametar koji ocenjujemo
# n - obim uzorka
# sigma - koren iz sigma^2 (koja nam je poznata u ovom slucaju)
# I u tom slucaju mozemo reci da je prema Centralnoj 
# Granicnoj Teoremi (CGT): Z~N(0,1) , tj Z je normalna sa parametrima 0 i 1.

# Dok u slucaju kada ne znamo sigma^2, moramo je oceniti i tada imamo 
# sledecu statistiku pomocu koje nalazimo Interval poverenja:
# T=(X_sr-m)*sqrt(n)/s_n
# gde je 
# s_n =sqrt(1/(n-1) * sum[(x_i-X_sr)^2]) - tj ocena korena iz disperzije.
# Tada T~t_{n-1} - tj Studentovu raspodelu sa n-1 stepeni slobode.


# To je sve bilo za slucaj kada se ocenjuje parametar m. Dok ako zelimo da intervalno ocenimo
# sigma^2 tu je prica malo drugacija.

# Statistika koja se koristi za to je:
# W=(n-1)*s_n_kv/sigma^2 ~ Hi-kvadrat sa n-1 stepeni slobode.
# gde je:
# n - obim uzorka
# s_n_kv = var(uzorak) = 1/(n-1)*sum[(x_i-x_sr)^2]
# sigma^2 - parametar koji ocenjujemo.
# I sasvim je prirodno da ta statistika, kojom se ocenjuje parametar sigma^2 
# ima neku nenegativnu raspodelu, kao sto je Hi-kvadrat.

# Za razliku od ocene za "m", gde su u oba slucaja statistike imale 
# simetricne raspodele oko nule, hi-kvadrat nije simetricna, i stavise
# nenagativna je.
# Iz tog razloga na vise nacina se moze prirodno definisati interval poverenja za sigma^2.
# kao na primer:

# Napomena: sve granice intervala koje cu da navodim dole, tj (A_n, B_n, U_n, V_n) su zapravo granice za statistiku(!) ne i za parametar.
# A tek kad se dobiju granice za stistiku - moramo da nadjemo bas Interval poverenja za trazeni parametar sigma^2.

# 1)
# Donji interval, tj (0,A_n], tj da pocinje uvek od nule, pa ide do vrednosti gde se dostigne trazena Beta, gde se A_n nalazi bas iz:
# P{W<=A_n}=beta => A_n=F^(-1)(beta), gde je F^(-1) inverzna funkcija za funkciju raspodele koju ima statistika W.
# sto se u R-u poziva na sledeci nacin:
# A_n=qchisq(beta, n-1)
# gde je qchisq - funkcija pomocu koje dobijamo kvantile hi-kvadrat raspodele. 

# 2) Dvostrani interval:
# [U_n, V_n], gde se granice biraju tako da je levo i desno od tog intervala ista povrsina koja je jednaka (1-beta)/2, za dato "beta".
# Nalazimo ih pomocu:
# P{W<=U_n}=(1-beta)/2 => U_n=qchisq((1-beta)/2, n-1).
# P{W<=V_n}=(1+beta)/2 => V_n=qchisq((1+beta)/2, n-1). Ovde je (1+beta)/2, jer je prosto povrsina koja je pre vrednosti V_n mora da bude beta+(1-beta)/2
# sto je bas jednako (1+beta)/2.

# 3) Gornji interval:
# slicno kao donji, smo sto nije ogranicen sa gornje strane pa je interval [B_n, +besk).
# B_n iz formule:
# P{W<=B_n}=1-beta => B_n=qchisq(1-beta, n-1).




# (3)
# Neka je dato:
# X~ Norm(m, sigma^2).
# n=20
# s_n_kv=21.12.
# Naci dvostrani interval poverenja za parametar sigma^2 za beta=0.99

# resenje:
# Potrebno je naci interval oblika:
# P{U_n<=W<=V_n}=0.99
# Prema gorenavedenoj teoriji:
# vrednost za V_n mora za bude vrednost F^(-1)(1.99/2)
v_n=qchisq(1.99/2, 19)
# vrednost za U_n mora za bude vrednost F^(-1)(0.01/2)
u_n=qchisq(0.01/2, 19)
u_n
n=20
s_na_kvadrat=21.12
# Odakle dobijamo sledeci izrar iz kojeg moramo da nadjemo interval poverenja.
# P{6.85<=(n-1)*S^2_n/sigma^2<=38.6}
# Sto je dalje jednako:
# P{(n-1)*21.12/38.6<=sigma^2<=(n-1)*21.12/6.85}

# donja vrednost intervala poverenja:
a=(n-1)*21.12/38.6

# gornja vrednost intervala poverenja:
b= (n-1)*21.12/6.85
a
b
(interval=c(a,b))
# Sto i jeste trazeni interval poverenja za parametar sigma^2

# Probajte da uradite gornji i donji interval poverenja za ovaj primer.


# Jos jedna vrsta intervala poverenja, koju cemo da razmotrimo su intervali poveranja za parametar p u Bernulijevoj raspodeli. (Podsetite se te slucajne velicine)
# Neka je X~Ber(p)=Bin(1, p) i zelimo da ocenimo parametar p.
# Tada to mozemo oceniti intervalno pomocu statistike:
# Z=(X_sr-p)*sqrt(n)/sqrt(p*(1-p))~ Norm(0,1).
# I dalje je ocenjivanje jako slicno kao za parametar "m" sto smo imali, s tim sto
# ovde imamo "malo vise" javljanja parametra "p" u datoj statistici, nego sto je to bilo u slucaju parametra "m".
# Stoga u opstem slucaju moracemo da resavamo kvadratnu nejednacinu i da pazimo da granice ne izadju izvan intervala (0,1) - jer ocenjujemo parametar
# p - koji predstavlja verovatnocu nekog dogadjaja!

# Dok u nekim specijalnim slucajevima mozemo umesto p*(1-p) koristiti (X_sr*(1-X_sr)). Za to mora da bude ispunjeno:
#  - parametar p - ne sme da bude blizu 0 ili 1.
#  - n>=100, tj da imamo dovoljno mnogo elemenata uzorka.
# U tom slucaju vise nemamo kvadratne nejednacine vec direktne vrednosti koje cemo da ubacimo i brzo izrazimo i dobijemo trazeni interval poverenja.

# (4)
# Od 500 ljudi, 285 se izjasnilo protiv. Naci 90% interval poverenja za verovatnocu da se slucajno odabrana osoba izjasni protiv.

# resenje:
# Kako se ta statistika aproksimira normalnom raspodelom, odmah nadjemo trazenu vrednost za "c".
# primetimo da je (1+beta)/2=(1+0.9)/2=0.95
c=qnorm(0.95)
c

# Dalje imamo:
# P{-c<=(X_{sr}-p)/sqrt(p*(1-p))*sqrt(n)<=c}
# mozemo da aproksimiramo sa 
# X_sr*(1-X_sr), jer n=500>100 i p~0.57
n=500
x_sr=285/500
x_sr

# dalje nadjemo donju, odnosno gornju vrednost intervala:
a=x_sr-c*sqrt(x_sr*(1-x_sr))/sqrt(n)
b=x_sr+c*sqrt(x_sr*(1-x_sr))/sqrt(n)

# trazeni interval poverenja za "p":
(interval=c(a,b))

# (5)
# Istrazivac procenjuje procenat gojaznih osoba sa visokim krvnim pritiskom.
# Od 25 ispitanika - 20 ima visok pritisak.
# Naci 95% interval poverenja za parametar p, koji predstavlja verovatnocu 
# da gojazna osoba ima visok krvni pritisak.

# resenje:
# ovde ne mozemo da primenimo aproksimaciju, jer n=25<100.
# Zato moramo da resimo kvadratnu nejednacinu. Ali pre toga vrednost za c:
beta=0.95
c=qnorm((1+beta)/2)
c

#
x_sr=20/25
# P{-c<=(X_{sr}-p)/sqrt(p*(1-p))*sqrt(n)<=c}
# = P{-1.96<=(0.8-p)/sqrt(p*(1-p))*5<=1.96}
# = P{|0.8-p|<=1.96*sqrt(p*(1-p))/5} sto sada kvadriramo:
# = P{0.8^2-1.6*p+p^2<=1.96^2*(p-p^2)/25}
# odakle klasicnim metodama dobijamo: 
#(vrednosti su zaokruzene na dve decimale)
a=0.61
b=0.91
(interval=c(a,b))

# (6)
# Naci 90% i 80% interval poverenja za istrazivanje u kojem je 
# ucestvovalo 100 ljudi
# i srednja vrednost je 0.45.

# Mozemo iskoristiti aproksimacije p*(1-p) = X_sr(1-X_sr)

# Z = (X_sr - p)/ sqrt[X_sr(1-X_sr)] * sqrt(n)
c_90 = qnorm((1+0.9)/2) # 1.64
c_80 = qnorm((1+.8)/2) # 1.28

(interval_90 = c(0.45 - c_90*sqrt(0.45*.55) / 10, 0.45 + c_90*sqrt(0.45*.55) / 10 ) ) 
(interval_80 = c(0.45 - c_80*sqrt(0.45*.55) / 10, 0.45 + c_80*sqrt(0.45*.55) / 10 ) ) 


# (7)
# Student zeli da istrazi procenat levorukih osoba. 
# Ispitano je 100 ljudi od kojih je samo 5 levorukih.
# Ako se nadje 90% interval poverenja da li ce on sadrzati vrednost 0.1? 
# Ako ne - koliko-procentni interval poverenja treba da trazimo?

# Da li 0.1 pripada (a, b) , ako P(a<=p<=b)=0.9?
# p = 5%, sto jeste blizu nule. Ovde ne smemo da koristimo aproksimaciju 
# za disperziju uzorka. 

# Test statistika nam je Z = sqrt(n)*[X_sr - p] / sqrt[p*(1-p)]

# P(a<=p<=b) = 0.9
# P(-c <= Z <= c) = 0.9
c = qnorm((1+0.90)/2) # 1.64
# P(|X_sr - p| <= c*sqrt(p*(1-p)/n)) =
# P(X_sr^2 - 2_X_sr*p + p^2 <= c^2*p(1-p)/n) = 
# 
a = 0.025
b = 0.099
(interval_1 = c(a, b) )
# ne upada 0.1

# ali za beta = 0.91 vec upada jer je c tada
c = qnorm((1+0.91)/2) # 1.69
# i interval je (0.024, 0.10001) otprilike.

# (8)
# Od poslednjih 10 puta kada je tostirani sendvic pao - devet puta 
# je pao na "ukusnu stranu".
# Ako je uzorak iz Bernulijeve raspodele sa parametrom p, 
# naci 80% interval poverenja da sendvic pada na "ukusnu stranu".

# X_sr = 0.9 nije dovoljno daleko od 1, a i n=10 je bas malo.
# Test statistika je Z = [X_sr - p] / sqrt[p(1-p)/n]
# P{a<=p<=b} = P{-c <= Z <= c} = 0.8
c = qnorm(1.8/2) # 1.28
a = 0.72
b = 0.97
interval = c(a, b)


# Ukoliko ipak primenimo aproksimaciju:
#(interval = c(0.9 - c*sqrt(0.1*0.9/10),0.9 + c*sqrt(0.1*0.9/10) ) )
# Moramo da izbacimo ove vrednost preko 1 iz intervala poverenja.
#(interval = c(max(0, 0.9 - c*sqrt(0.1*0.9/10)), min(1, 0.9 + c*sqrt(0.1*0.9/10))))
# Dobijamo skroz drugacije intervale.


# (9)
# Izvrseno je ispitivanje 250 ljudi, ocenjeno p odgovora je 0.45.
# Isto istrazivanje je ponovljeno na 1000 ljudi i dobijen isti odnos.
# Naci 95% interval poverenja u oba slucaja. Koliko puta je drugi manji od prvog?
# Koliko veliko treba da bude n u istrazivanju da bi 95% interval poverenja 
# bio duzine ne vece od 0.01.

# Z = [X_sr - p] / sqrt[X_sr(1-X_sr)/n]

c_95 = qnorm(0.975) # 1.96

interval_250 = c(0.45 - c_95*sqrt(0.45*0.55/250), 0.45 + c_95*sqrt(0.45*0.55/250))
# (0.39, 0.51)
diff(interval_250) # 0.12334

interval_1000 = c(0.45 - c_95*sqrt(0.45*0.55/1000), 0.45 + c_95*sqrt(0.45*0.55/1000))
interval_1000
# (0.42, 0.48)
diff(interval_1000) # 0.06167

# Bolje u obliku funkcije da se zapise.
duzina_intervala <- function(n)
{
  interval_1000 = c(0.45 - c_95*sqrt(0.45*0.55/n), 0.45 + c_95*sqrt(0.45*0.55/n))
  return(diff(interval_1000))
}
duzina_intervala(38031)

# Trazeno n nalazimo na sledeci nacin:
(duzina_intervala(250)/0.01)^2*250
# ~38031

# (10)
# Bajesovski intervali pokrivanja/prekrivanja (credible interval).

# X iz Bernulijeve sa parametrom p. Dobijen PSU duzine n=250
# Apriorna za p je Beta(2, 3).
# Tada aposteriorna je Beta(sum(x)+2, n - sum(x)+3)
# Naci interval prekrivanja za beta= 0.9

x = sample(x = c(0,1), size = 250, replace = TRUE, prob = c(0.45, 0.55))
suma = sum(x)
suma
# Interval prekrivanja je P(a<= p <= b | X) = 0.9
b = qbeta(0.95, suma + 2, -suma + 3 + 250)
a = qbeta(0.05, suma + 2, -suma + 3 + 250)
interval = c(a, b)
interval
# prava vrednost je bila 0.45. a srednja 0.456.