# 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"))