examples/Lecture2a.R (6b93c51) - SS2012-StatisticalMethods.git

Lecture2a.R

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

SS2012-StatisticalMethods.git

new fiels