# EXAMPLE 3: Introduction to Statistics for Astrophysicists
# SS 2012
# JOCHEN WELLER




# Sequence of probabilities (Hypothesis)
H=seq(0,1,by = 0.01)

# number of trials
n=32

# bias 
bias=0.25

# generate random sample with possible outcome 0,1 and bias: toss the (un-)fair coin
coin=sample(c(0,1),n,replace=TRUE,c(1.0-bias,bias))

# 

# count heads of sample
heads=sum(coin)


# prior distribution 
# uniform
prior = 1.0


#Gaussian approximation
H0 = heads/n
sigma=sqrt(H0*(1-H0)/n)
likeapprox <- function(H) 1.0/sqrt(2.0*pi)/sigma*exp(-0.5*(H-H0)**2/sigma**2)

#calculate normalization of binomial distribution
sum = 0.0
i = 0
while (i < n-heads) {sum=sum+choose(n-heads,i)*((-1)^i)/(heads+i+1) 
	i=i+1 }
norm = abs(1.0/sum)

# plot posterior (likelihood times prior)

par(font.lab=2)
par(font.axis=2)
plot(H,norm*dbinom(heads,n,H)*prior,type='l', xlab='H',ylab='prob(H|{data},I)', lwd=2, main=paste("N=",n))
lines(H,likeapprox(H)*prior,lty=2,lwd=2)
abline(v=H0,lty=3,lwd=2,col='blue')
abline(v=H0-sigma,lty=3,lwd=2,col='blue')
abline(v=H0+sigma,lty=3,lwd=2,col='blue')
legend("topright",legend=c("Binomial","Gaussian approximation","maximum, sigma"),lty=c(1,2,3),lwd=c(2,2,2),col=c("black","black","blue"))