Michi
new fiels
Michi commited 6b93c51 at 2012-06-03 16:52:14
# 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')
 
 
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX