library(Cairo)
library(VGAM)
library(rgl)

setwd("/home/martin/Lectures/JourFixe")

## Die Binomialverteilung
plot(	c(0:10),
	dbinom(c(0:10),10,0.3),
	t="h",
	xlab="k",
	ylab="Prob(k|n=10,p=0.3)")

## Welche Urne war es?
dbinom(	4,
	10,
	c(0.3,0.5,0.7),
	)

## Die Binomiale Likelihoodfunktion
curve(	dbinom(4,10,x),
	from=0,
	to=1,
	xlab="p",
	ylab="Lik(p|k=4,n=10)")
abline(v=0.4,col="red")

## Schätzung bei K=c(2,2,3,4)
K<-c(2,2,3,4)
lik.binom<-function(p,K,n) prod(dbinom(K,n,p))
p<-seq(0,1,0.01)
lik<-sapply(p,lik.binom,K,10)
lik.max<-p[lik==max(lik)]

plot(	p,
	lik,
	t="l",
	xlab="p",
	ylab="Lik(p|K=(2,2,3,4),n=10)")
abline(v=0.28,col="red")
text(x=0.28,y=c(0,max(lik)),labels=c(paste("p=",lik.max),paste("max L(p|K,10)=",round(max(lik),6))),pos=4)

## Nicht-konstanter Fehler bei der Binomialverteilung
freq<-function(issues,range) sapply(lapply(range,"==",issues),sum)
plot(c(0:10),dbinom(c(0:10),10,0.3),t="h",ylim=c(0,0.5))
result<-c()
for (s in c(1:100)){
	d<-rbinom(100,10,0.3)
	result<-rbind(result,freq(d,c(0:10)))
}
for (f in c(0:10)){
	points(x=rep(f,100),y=result[,f+1]/100,col="grey")
}

## Likelihoodfläche, wenn das Binomiale p und n geschätzt werden soll
##---

## Die Betafunktion
curve(dbeta(x,0.1,2),from=0,to=1,col="red")
curve(dbeta(x,2,12),from=0,to=1,add=TRUE,col="green")
curve(dbeta(x,1,2),from=0,to=1,add=TRUE,col="blue")
curve(dbeta(x,2,4),from=0,to=1,add=TRUE,col="orange")
curve(dbeta(x,2,6),from=0,to=1,add=TRUE,col="violet")

## Die Zeckendaten
zecken1<-c(rep(0,5),rep(1,4),rep(2,6),rep(3,3),4,5)
zecken2<-rbetabin.ab(20,10,2,4)
zecken<-zecken2
## Die Likelihoodfunktion der betabin Verteilung

lik.betabin<-function(n,K,a,b) prod(dbetabin.ab(K,n,a,b))

## Die Likelihoodoberfläche der Zeckendaten
lik<-c()
## brute force
for (a in seq(2,6,0.1)){
	for (b in seq(5,25,0.2)){
		lik<-rbind(lik,c(a,b,lik.betabin(10,zecken,a,b)))
	}
}

max.lik<-lik[lik[,3]==max(lik[,3]),]
print(lik[lik[,3]==max(lik[,3]),])

plot3d(lik,t="l")
points3d(x=max.lik[1],y=max.lik[2],z=max.lik[3],col="red",size=5)
text3d(x=max.lik[1],y=max.lik[2],z=max.lik[3],col="red",pos=4,texts="MAX")

## Intelligent
## Die negative Loglikelihood
nloglik.betabin<-function(p,n,K) -sum(log(dbetabin.ab(K,n,p[1],p[2])))
nloglik.binom<-function(p,n,K) -sum(log(dbinom(K,n,p)))
fitbb.zecken<-optim(
	fn=nloglik.betabin,
	par=c(0,0),
	K=zecken,
	n=10,
	lower=c(0.1),
	method="L-BFGS-B"
	)

fitb.zecken<-optim(
	fn=nloglik.binom,
	par=0.01,
	K=zecken,
	n=10,
	lower=0.1,
	upper=0.9,
	method="L-BFGS-B"
	)

## Vergleich des AIC

AIC<-function(L,k,n,corr=TRUE){
	AIC<-2*k-2*L
	cAIC=AIC+((2*k)*(k+1))/(n-k-1)
	return<-ifelse(corr==TRUE,cAIC,AIC)
	return
}

AIC(-fitb.zecken$value,1,20)
AIC(-fitbb.zecken$value,2,20)