##############################################
# Metoda momentu - denni srazky
##############################################
set.seed(95478)
alpha=0.35
beta=10
X=rgamma(1000,alpha,scale=beta) #nahodny vyber z gamma rozdeleni
plot(X,type="l")
hist(X,prob=T,n=50)
curve(dgamma(x,alpha,scale=beta),add=T,col=2,lwd=2)

alpha_vlnka=mean(X)^2/var(X) #odhad parametru alpha metodou momentu
beta_vlnka=var(X)/mean(X)    #odhad parametru beta metodou momentu
alpha_vlnka
beta_vlnka

curve(dgamma(x,alpha,scale=beta),add=T,col=3)



###################################################
# MLE pro intervalova data
################################################
Y=c(99,42,29,28,17,9,3)
c=c(0,7500,17500,32500,67500,125000,300000,Inf)

#verohodnostni funkce pro nase data
f=function(x){
f=0
for (i in 1:7){
f=f+Y[i]*log(pexp(c[i+1],1/x)-pexp(c[i],1/x))
}
return(f)
}

opt=optimize(f,c(0,1000000),maximum=T)$max  # hledani maxima v intervalu (0,1 000 000)
opt  # MLE parametru lambda, odhad prumerne vyse plneni jednoho klienta

X=(7500/2*99+(7500+17500)/2*42+29*(17500+32500)/2+28*(32500+67500)/2+17*(67500+125000)/2+9*(125000+300000)/2+3*400000)/227
X    #prumerna hodnota z dat (bez uvazovani jakehokoliv modelu)

##############################################
# Pridame odhadnutou hustotu do "histogramu"
###############################################
br=c(0,7500,17500,32500,67500,125000,300000,400001)  #hrany histogramu 
data=c(rep(7000,99),rep(15000,42),rep(30000,29),rep(50000,28),rep(100000,17),rep(200000,9),rep(400000,3)) #pomocna data
hist(data,breaks=br) #"histogram" nasich dat
curve(dexp(x,1/opt),add=T,col=3) #pridame hustotu odhadnuteho exp. rozdeleni



###############################################
# Metoda minimalniho chi^2
###############################################
Y=c(99,42,29,28,17,9,3)
c=c(0,7500,17500,32500,67500,125000,300000,Inf) #opet uvazujeme data vyse pojistnych plneni
n=sum(Y)

chi=function(x){ #vypocet Pearsonovy chi^2 statistiky
p=rep(0,7) #teoreticke cetnosti
for (i in 1:7){ 
p[i]=pexp(c[i+1],1/x)-pexp(c[i],1/x)}

chi=sum((Y-n*p)^2/(n*p))
return(chi) 
}

opt2=optimize(chi,c(20000,100000))$min
opt2

##############################################
# Pridame odhadnutou hustotu do "histogramu"
###############################################
br=c(0,7500,17500,32500,67500,125000,300000,400001)
data=c(rep(7000,99),rep(15000,42),rep(30000,29),rep(50000,28),rep(100000,17),rep(200000,9),rep(400000,3))
hist(data,breaks=br)
curve(dexp(x,1/opt2),add=T,col=3)


#############################################
# Modifikace predchoziho prikladu
#############################################
c=c(0,7500,17500,32500,67500,125000,300000,400001) #temer puvodni data
Y=c(99,42,29,28,17,9,3)	#puvodni data

#vygenerujeme si nova data, ktera nebudou intervalova  
X=runif(Y[1],c[1],c[2])
for (i in 2:7){
x=runif(Y[i],c[i],c[i+1])
X=c(X,x)
}
n=length(X)

X  # nova data
hist(X,breaks=br)  #jejich "histogram" odpovida puvodnim datum

k=round(15*(n/100)^(2/5)) #doporucena volba k
k 
XX=mean(X)  # maximalne verohodny odhad lambda, vezmeme ho jako startovaci odhad

# najdeme intervaly B_1,...B_k
cc=0
for (i in 1:k){
cc[i]=qexp(i/k,1/XX)
}
cc=c(0,cc)

#spocitame cetnosti pozorovani v danych intervalech Y_1,...Y_k
YY=0
for (i in 1:k){
YY[i]=sum(X<cc[i+1])-sum(X<cc[i])
}

# nyni opet provedeme odhad parametru lambda metodou minimalniho chi^2
chi=function(x){
p=0
for (i in 1:k){ 
p[i]=pexp(cc[i+1],1/x)-pexp(cc[i],1/x)}
chi=sum((YY-n*p)^2/n/p)
return(chi)
}

opt3=optimize(chi,c(20000,100000))$min
opt3

##############################################
# Pridame odhadnutou hustotu do "histogramu"
###############################################
br=c(0,7500,17500,32500,67500,125000,300000,400001)
data=c(rep(7000,99),rep(15000,42),rep(30000,29),rep(50000,28),rep(100000,17),rep(200000,9),rep(400000,3))
hist(data,breaks=br)
curve(dexp(x,1/opt3),add=T,col=3)
























