### Poredjenje kvaliteta dvije ocjene za isti parametar ###

# Pretostavimo da imamo uzorak X=(X1,X2,...,Xn) iz neke raspodjele koja zavisi
# od nepoznatog parametra theta. Neka su date dvije ocjene, Tn i Vn za taj
# parametar theta. Na predavanjima iz statistike ste culi o asimptotskoj
# efikasnosti ocjena. Medjutim, kad imamo male obime uzorka nemamo asimptotska
# svojstva, pa moramo koristiti neku drugu metodu. Poznato nam je da je jedna od
# mjera kvaliteta ocjene srednjekvadratna greska (MSE), pa cemo nju koristiti. 
# Metod se sastoji u sledecem:

# 1. Za neko fiksirano theta generisemo uzorak obima n iz date raspodjele.
# 2. Racunamo vrijednosti statistike Tn i Vn za taj uzorak.
# 3. Postupak ponovimo veliki broj puta (N) i racunamo srednje kvadratnu gresku
# ocjena Tn i Vn
# 4. Uporedimo greske, ocjena koja ima manju srednjekvadratnu gresku je bolja za
# taj parametar.

# Ovim dobijamo koja je greska bolja za jednu konkretnu vrijednost parametra
# theta, pa onda postupak treba ponoviti i za neke druge vrijednosti i uporediti
# rezultate.

# Napomena: Srednje kvadratnu gresku racunamo iz uzorka pa se ovaj metod koristi
# najvise kad ne mozemo lako odrediti raspodjelu statistike.


# Primjer: Neka je X iz U(0,theta) raspodjele, gdje je theta nepoznati
# parametar. Predlozene su dvije statistike za ocjenu theta : Tn=2Xn_sr i
# Vn=(n+1)/n*X_(n) i hocemo da ispitamo koja je bolja.

# Obim uzorka

n<-20

# broj simulacija

N<-10000

theta.0<-1

Tn<-function(x) 2*mean(x)

Vn<-function(x) ((length(x)+1)/length(x))*max(x)

# Pravimo matricu koja u svakom redu ima jedan uzorak obima n i sastoji se od 
# N takvih redova, tj. N puta generisemo uzorak iz U(0,theta) 
m<-matrix(runif(n*N), ncol = n, nrow = N)

# Racunamo srednje kvadratne greske za Tn i Vn

MSE.Tn<-mean((apply(m, 1,function(x) Tn(x))-theta.0)^2)
MSE.Vn<-mean((apply(m, 1,function(x) Vn(x))-theta.0)^2)

MSE.Tn>MSE.Vn

# Ponovimo sad isti postupak za razlicite vrijednosti theta0

poredjenje<-function(theta.0,n=20,N=10000){
  m<-matrix(runif(n*N,0, theta.0), ncol = n, nrow = N)
  MSE.Tn<-mean((apply(m, 1,function(x) Tn(x))-theta.0)^2)
  MSE.Vn<-mean((apply(m, 1,function(x) Vn(x))-theta.0)^2)
  MSE.Tn>MSE.Vn
}

theta.0<-seq(0,10)

sapply(theta.0, function(x) poredjenje(x))

# Napomena: Ovaj primjer je bio samo ilustracija ovog metoda. Ove dvije ocjene
# su nepristrasne i njihove disperzije se mogu lako izracunati. Disperzija druge
# ocjene je znato manja pa smo zbog toga dobili da je u svim slucajevima druga
# ocjena bila bolja. Za neki naredni domaci ili seminarski rad radicete primjere
# kod kojih ima vise smisla koristiti ovu metodu .



### Metoda maksimalne vjerodostojnosti ###

# Za neke klase raspodjela u R-u postoji funkcija koja za proslijedjeni uzorak
# daje ocjenu metodom mv

# Trazena funkcija je fitdistr(x,) iz paketa MASS

x <- rpois(1000, lambda = 3)
install.packages("MASS")
library(MASS)
fitdistr(x, "Poisson")
mean(x)


?fitdistr

# ``Distributions "beta", "cauchy", "chi-squared", "exponential", "f", "gamma",
# "geometric", "log-normal", "lognormal", "logistic", "negative binomial",
# "normal", "Poisson", "t" and "weibull" are recognised, case being ignored.``

n<-1000

x<-rnorm(n)
fitdistr(x,"normal")
sd(x)

# razlikuje se jer ocjena mmv za sigma2 nije popravljena uzoracka disperzija

sd.mmv<-sqrt((n-1)/n*var(x))
sd.mmv


# Nas interesuje slucaj kad ne mozemo ekspicitno izracunati ocjenu mmv 


#Njutn-Rafson (Newton-Raphson)
#U 1D slucaju zove se samo Njutnova metoda ili metoda tangente
#pogledati knjigu "Numericke metode", Desanka Radunovic, 159-177 str

#PRIMER 1
#Kosijeva raspodela

#generisemo uzorak iz Kosijeve raspodele sa parametrom polozaja 10

x<-rcauchy(100,location=10)
hist(x) 
median(x) 
#da vidite bolje osobine uzorka, parametar lokacije (polozaja) je u stvari 
#medijana raspodele, ocekivanje ne postoji 
#ovo je ocena tog parametra metodom zamene mi hocemo metodom mv

#mi sada znamo da smo generisali tacno iz Kosijeve raspodele sa parametrom 10 
#ali cemo se praviti da to ne znamo i pokusati Njutnovom metodom da ocenimo
#parametar raspodele trebalo bi da dobijemo nesto blizu 10 (ne mozemo ocekivati
#tacno 10 zbog varijabilnosti uzorka) tako sto cemo naci nulu prvog izvoda
#log-likelihood funkcije (score function) (ona zavisi od uzorka) koja se inace
#ne moze izraziti eksplicitno

#prvo mozemo nacrtati funkciju ciji maksimum trazimo
cauchy_loglike<-function(theta,x)
{
  return(sum(-log(1+(x-theta)^2))-length(x)*pi)
}
n<-50
loglike_array<-vector()
theta<-seq(from=8,to=12,length.out=50)
for(i in 1:n)
{
  loglike_array<-c(loglike_array,cauchy_loglike(theta[i],x))
}

plot(theta,loglike_array,type="l")
#pitanje moze biti zasto onda samo ne nadjemo maksimum ovog niza i to proglasimo
#ocenom MMV mozemo, ali mi mozda nadjemo precizniju ocenu maksimuma, posto
#loglike_array ima samo nekoliko tacaka

#inace ne moramo ni crtati ovu funkciju, ili mozemo u jos manje tacaka (manje od
#50) da bi stedeli vreme ali ona moze biti korisna za izbor inicijalnog elementa
#niza


#score funkcija i izvod score funkcije koji se koriste u iteracijama se moraju naci rucno
#pravimo prvo funkcije u Ru koje ce racunati te vrednosti
cauchy_score<-function(theta,x)
{
  return(sum(2*(x-theta)/(1+(x-theta)^2)))
}
#cauchy_score(10,x)

d_cauchy_score<-function(theta,x)
{
  return(2*sum(((x-theta)^2-1)/(1+(x-theta)^2)^2))
}
#d_cauchy_score(10,x)

#mozemo nacrtati i score funkciju, umesto loglikelihood, da bismo videli gde je
#otprilike nula izvoda sto moze onda biti tacka lokalnog maksimuma ili minimuma
#loglikelihood funkcije

n<-50
score_array<-vector()
theta<-seq(from=8,to=12,length.out=50)
for(i in 1:n)
{
  score_array<-c(score_array,cauchy_score(theta[i],x))
}

plot(theta,score_array,type="l")
#vidimo da je score funkcija 0 u nekoj tacki blizu 10

#primetimo da smo ovde crtali na malom intervalu (od 8 do 12) jer smo znali iz
#cega smo generisali uzorak trebalo bi na nekom vecem intervalu ako nemamo
#nikakve informacije o parametru

#ovde, posto je theta bas medijana - mozemo uzoracki lako oceniti parametar i to
#koristiti kao pocetnu vrednost theta0, ne treba nam crtanje loglikelihood i
#score funkcije

eps<-0.003 #tacnost
theta<-vector() #vektor iteracija, poslednji ce biti krajnje resenje
theta0<-11 #pocetna vrednost
theta[1]<-theta0

i=1
repeat
{
  theta[i+1]<-theta[i]- cauchy_score(theta[i],x)/d_cauchy_score(theta[i],x)
  i=i+1
  if(abs(theta[i]-theta[i-1])<eps) break
} #ponavljamo dok razlika izmedju uzastopnih ne bude mala
#sto se desava zbog konvergencije metode

#KONACNO:
theta[i] #OCENA MMV ZA PARAMETAR POLOZAJA KOD KOSIJEVE RASPODELE DOBIJENA NUMERICKOM
#METODOM

#Jedna modifikacija Njutnove metode, koja smanjuje potreban racun je da se
#umesto racunanja izvoda fje u svakoj iteraciji (odnosno izvoda score funkcije u
#nasem slucaju) koristi izracunati izvod u inicijalnoj tacki probajte da
#napravite funkciju koja na taj nacin racuna ocenu metodom MMV


#Fiserov algoritam

#specijalno nas problem, koriscenje za nalazenje ocena MMV umesto inverza
#matrice drugih parcijalnih izvoda log-likelihood funkcije koristimo njenu
#ocekivanu vrednost (ako je stvarno vrednost parametra teta - to koje se
#pojavljuje u izrazu za loglikelihood)

#dakle racunanje matrice u svakoj iteraciji i njenog inverza menjamo racunanjem
#jednog broja po iteraciji

#PRIMER
#Kosijeva raspodela sa parametrom lokacije, ocekivana vrednost Fiserove
#informacione matrice je n/2 gde je n obim uzorka

eps<-0.00001 #tacnost
theta<-vector() #vektor iteracija, poslednji ce biti krajnje resenje
theta0<-9 #pocetna vrednost
theta[1]<-theta0

i=1
repeat
{
  theta[i+1]<-theta[i]+ 4*sum((x-theta[i])/(1+(x-theta[i])^2))/length(x)
  i=i+1
  if(abs(theta[i]-theta[i-1])<eps) break
} 
theta[i]
i

i=1
repeat

{
  theta[i+1]<-theta[i]+ cauchy_score(theta[i],x)*2/n
  i=i+1
  if(abs(theta[i]-theta[i-1])<eps) break
} 
theta[i]
i