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


# The Lighthouse Problem
alpha0 = 1 #km true position at coastline
beta = 1 #km position out at see

# Sequence of positions (use already prior range from -10 to 10 km)
alpha=seq(-10,10,by = 0.05)

# number of trials
N=512

# generate random sample distributed according to a Cauchy distribution (see lecture)
# with N trial
pos=rcauchy(N,alpha0,beta)

# prior distribution 
# uniform between -10 and +10 km
prior = 1.0/20.0

sum = 0.0
i = 1

logpost <- function(alpha) { while (i<=N){sum=sum+log(beta^2+(pos[i]-alpha)^2)
		i=i+1 }; sum}

#find minimum of log likelihood
lmin=min(logpost(alpha))

# normalized posterior
post <- function(alpha) exp(-logpost(alpha)+lmin)

# plot posterior 

par(font.lab=2)
par(font.axis=2)
plot(alpha,post(alpha),type='l', xlab=expression(paste(alpha,"(km)")),ylab=expression(prob(~alpha~"|{"~x[k]~"}",~beta,I)), lwd=2, main=paste("N=",N),xlim=c(-12,12),ylim=c(0,1.2))

# plot points of samples at vertical position 1.1

normpoint = (pos-pos)+1.1
points(pos,normpoint)

# vertical line at mean of sample
abline(v=mean(pos),lwd=2,lty=2)