# PRVI CAS

# Prvo nesto kratko o MonteKarlo metodi. 
# Ideja je da se generisu neki slucajni brojevi, na odredjeni nacin ispitaju neophodni uslovi (tj fakticki gledamo koliko cesto se nesto "realizovalo")
# I vratimo taj broj. 

# Ta metoda specijalno dobro prolazi u slucaju integracije i pogodna je ako imamo neke nezgodne podintegralne funkcije, koje je tesko ili nemoguce izraziti u elementarnim
# funkcijama. U slucaju da se integracija ne vrsi na skupu [0,1] vec na nekom drugom, moramo da pazimo kako biramo granice za generisanje pseudoslucajnih
# brojeva. Promenljiva "x" se svakako krece po granicama integracije. dok "y" od "min f(x)" do "max f(x)" po skupu svih mogucih vrednosti za "x". 
# Zato je neophodan uslov da podintegralna funkcija f(x) bude ogranicena na datom skupu.



#Pokazacemo na primeru nekih prostih funkcija.


# Zelimo da odredimo povrsinu ispod neke odredjene funkcije na intervalu [0, 1] pomocu Monte Karlo metode.

Montekarlo <- function(N,s) # s - stepen funkcije, s>0.
{
  brojac=c()
  for(i in 1:N)
  {
    x=runif(1, 0, 1)
    y=runif(1, 0, 1)
    # To su nase koodinate tacke
    
    if(y<x^s) #a mogli smo da stavimo i cos(x), sin(x) itd. 
    {
      brojac=c(brojac, 1)
    }
    else
    {
      brojac=c(brojac, 0)
    }
  }
  return(sum(brojac)/N)
}
Montekarlo(10000, 5)

# A sve smo to mogli i u par redova da napisemo, kako je i bolje raditi u R-u generalno
n=100000
s=5

Montekarlo2<-function(n,s)
{
  x=runif(n, 0, 1)
  y=runif(n, 0, 1)
  z=ifelse(y<x^s, 1, 0)
  return(mean(z))
}
Montekarlo2(n*10, s)

# Mozete da primetite da ova druga funkcija i radi cak mnogo brze nego li prva. Zato po mogucstvu uvek cemo da izbegavamo "for" petlje i generalno 
# pokusacemo da radimo vektorski.

# LINEARNA REGRESIJA

# Kako nam je cilj da odradimo neku logisticku regresiju, tacnije vise njih, prvo moramo da se upoznamo sa malo laksim idejama. A to su linearne regresije.
# Radicete to detaljno na kursu LSM, ali opet zbog kompletnosti price, bice i ovde malo.

# 1) Linearne regresije:
# Neka su data dva vektora "x" i "y", koji predstavljaju sledece:
# protok saobracaja u hiljadama automobila u toku dana
x=c(8.3, 9.3, 12.1, 12.3, 17.0, 17.3, 24.3, 24.5, 33.6)

# prisustva olova u kori drveca:
y=c(227, 312, 362, 521, 539, 728, 945, 1000, 1263)
# Oceniti koeficijente modela y=a*x+b, gde su x, y - vektori, dok a i b - neki brojevi. 

# Resenje:
# primenom linearne regresije sa srednje kvadratnom greskom:
model=lm(y~x)
summary(model)
plot(model)

a=41.318 # a mogli smo i da pozovemo a=model$coefficients[2], odnosno b=model$coefficients[1]
b=-73.356

plot(x, y)
# regresiona prava
abline(b,a, col="red")
# Tako dodajemo regresionu pravu na grafik. Da biste se detaljnije upoznali sa funkcijom abline, ukucajne ?abline.

# predvidjanje na osnovu ovog modela za vrednosti x=5, 13.3, 28, 30 i 100.
predict(model, newdata=data.frame(x=c(5, 13.3, 28, 30, 100)))

# 2)
# broj sati ucenja
x=c(4,9, 10, 14, 4, 7, 12, 22, 1, 3, 8, 11, 5, 6, 10, 11, 16, 13, 13, 10)

# rezultati testa.
y=c(390, 580, 650, 730, 410, 530, 600, 790, 350, 400, 590, 640,450, 520, 690, 690, 770,700, 730, 640)
# Potrebno je opet pokrenuti linearnu regresiju i odrediti koeficijente modela, kao i prokomentarisati da li je model kvalitetan.

# Resenje:
model=lm(y~x)
summary(model)
# Kao sto vidimo u ovom slucaju oba koeficijenta su od znacaja (tj nisu jednaki nuli, prema testu). 
# Kao i da je vrednost za R^2 poprilicno visoka, sto pokazuje da je model dovoljno kvalitetan. Ili kao sto se cesto koristi u literaturi: mozemo reci da 
# nas model na 86% opisuje promenljivu Y.

plot(x, y)

# regresiona prava
a=25.326
b=353.165
abline(b,a, col="red")
# Vidimo da postoji skoro da linearna zavisnost medju brojem sati ucenja i rezultatima.



# 3) Jos jedan primer, kada linearna regresije u svom pocetnom obliku ne prolazi dobro.
# Neka su date q i w na sledeci nacin:
q=seq(-1, 1, 0.05)
w=q^2+rnorm(length(q),0, 0.05)
# naci regresionu pravu za w u odnosu na q.

# Resenje:
model=lm(w~q)
summary(model)
# Vidimo da nemamo "zvezdice" za q - nije znacajna, i p-vrednost je blizu 1
# Vrednosti za R-squared treba da budu sto blizi jedinici dok u ovom slucaju su oko 0, sto ukazuje na totalni promasaj linearne regresije.

# Pogledajmo kako to izgleda na grafiku:
plot(q, w)
a=   model$coefficients[2]
b= model$coefficients[1]
abline(b,a, col="red", pch=20)
# Ovakve stvari se desavaju iz prostog razloga sto veza nije linearna.


# Probamo kvadratni model:
model2=lm(w~I(q)+I(q^2))
summary(model2)
plot(w~q)
points(q, fitted(model2), col="red",pch=20)
# Ovde je od znacaja samo parametar uz q^2, sto i ima smisla skroz (prosto smo tako kreirali vektor w). I vidimo da je regresija ispala dovoljno dobra.


# Nadalje bice nam potreban paket MASS, inace je od velikog znacaja, mnogo razlicitih testova, nacina ocenjivanja i slicnog se nalazi bas tu.
library(MASS)

# I pokusacemo da instaliramo paket ISLR (koji detaljno prati knjigu Hastie-ja i Tibshirani-ja).
install.packages("ISLR")
library(ISLR)

# Obicna linearna regresija na primeru baze- Boston.
Boston

# funkcija "names" vraca nazive promenljivih u datoj bazi.
names(Boston)

# ?Boston=help(Boston), tj opisuje sta predstavlja ta baza. Vezana je za ucestalost krsenja zakona u Bostonu. 
?Boston

# plotujemo sledece i dobijamo grafik raspsenosti promenljive medv u odnosu na lstat. 
plot(medv~lstat,Boston)

# Sledeci poziv se vec odnosi na linearni model, funkcija lm = linear model vraca najbolji linearni model koji odgovara datim podacima.
fit1=lm(medv~lstat,data=Boston)
fit1
# Od znacaja su nam ocene koeficijenata.
# Da bismo detaljnije videli kako i sta se desava, pozovemo summary().
summary(fit1)
# I vidimo da su obe promenljive znacajne.

# Sledecim pozivom mozemo da dodamo regresionu pravu na grafik, a ne samo kao sto to radili na prvom casu.
abline(fit1,col="red")
names(fit1)

# Ovde mozemo da vidimo intervali poverenja za ocene nasih koeficijenata modela.
confint(fit1)

# funkcija pairs() vraca sve moguce grafike rasprsenosti nase baze. Tj. pomocu nje mozemo ponekad i da ustanovimo gde odmah imamo zgodnu linearnu zavisnost
# i samim tim da primenimo linearnu regresiju. Jedini je problem, sto je zgodno predstaviti samo po 2 promenljive na grafiku,
# vec zavinost od 3 promenljive
# nije lako i lepo uocljiva na grafiku, a vece dimenzije cak i ne mozemo da nacrtamo. No, svejedno mozemo pokretati linearne regresije.

### Nadalje, linearna regresija se moze koristiti i na vise promenljivih.
fit2=lm(medv~lstat+age,data=Boston)
summary(fit2)

# Pri ovom zapisu znaci da smo uzeli u obzir sve moguce promenljive koje ima baza.
fit3=lm(medv~.,Boston)
summary(fit3)

# funkcija update(model, ...) menja model na zadati nacin. Prosto da ne bismo uvek pravili od nule.
fit4=update(fit3,~.-age-indus)
summary(fit4)

# I ne mora samo plus da bude izmedju prediktora, moze se javiti, recimo i proizvod: 
### NAPOMENA: ako stavimo x*y, tada se u model ukljucuje sledeci prediktori: x, y i x*y(koji se u pozvanom summary oznacava kao x:y).
# Ako bismo zeleli da pozovemo samo za prediktor x*y, tada bismo morali da pokrenemo:
# lm(z~ I(x*y)) ili u nasem slucaju fit=lm(medv~I(lstat*age))
fit5=lm(medv~lstat*age,Boston)
summary(fit5)

# Kad stavimo I(x^2) - time trazimo da x^2 posmatra totalno nezavisno od x. 
fit6=lm(medv~lstat +I(lstat^2),Boston); summary(fit6)
attach(Boston)
plot(medv~lstat)
points(lstat,fitted(fit6),col="red",pch=20)
# vidimo da je kvadratna regresija ovde bila preciznija od obicne. 

# Hocemo da probamo sa vecim stepenom mogli smo u fit6 da stavimo poly(lstat, 2) dobicemo isto.
fit7=lm(medv~poly(lstat,4))
points(lstat,fitted(fit7),col="blue",pch=20)
# Vidimo da je neznacajno bolja, ali zato mozda overfituje model, sto svakako bolje da izbegnemo.
summary(fit7)