##############################################
# metoda momentu - denni srazky
##############################################
set.seed(15)
X=rgamma(1000,0.35,scale=10)
plot(X,type="l")
hist(X,prob=T,n=50)
curve(dgamma(x,0.35,scale=10),add=T,col=2)

alpha=mean(X)^2/var(X)
beta=var(X)/mean(X)
alpha
beta

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



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

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(20000,100000),maximum=T)$max
opt

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*300000)/227


##############################################
# histogram
###############################################
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/opt),add=T,col=3)

###############################################
# metoda minimalniho chi^2
###############################################
Y=c(99,42,29,28,17,9,3)
n=sum(Y)
c=c(0,7500,17500,32500,67500,125000,300000,Inf)


chi=function(x){
p=0
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

##############################################
# histogram
###############################################
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
#############################################
c=c(0,7500,17500,32500,67500,125000,300000,400000)
Y=c(99,42,29,28,17,9,3)
X=0
for (i in 1:7){
x=runif(Y[i],c[i],c[i+1])
X=c(X,x)
}
X=X[-1]
n=length(X)

X
hist(X,breaks=br)

k=round(15*(n/100)^(2/5))
XX=mean(X)

cc=0
for (i in 1:k){
cc[i]=qexp(i/k,1/XX)
}
cc=c(0,cc)

YY=0
for (i in 1:k){
YY[i]=sum(X<cc[i+1])-sum(X<cc[i])
}

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)
}

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


























