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

x0=1
#fwhm = 5 = 2.35*sigma
w=2.13
n0=33.333
A=1
B=2

#First generate random data 
#Datum is Gaussian, and measurement Poisson

xmin=-7
xmax=7
x=seq(xmin,xmax,by=1)

d=n0*(A*exp(-0.5*(x-x0)^2/w^2)+B)
N=d
imax=length(d)
i=1
while (i <= imax) {N[i]=rpois(1,d[i])
	i=i+1}
	
# The minimum chi2 for later	
logepmin=sum(N*log(d)-d)	
	
quartz()

# Plot the data	
par(font.lab=2)
par(font.axis=2)	
plot(x,N,type='s',xlim=c(xmin,xmax),lwd=2)


# posterior function for arbitrary A and B
posterior<- function(A,B) {
	d=n0*(A*exp(-0.5*(x-x0)^2/w^2)+B)
	dummy=exp(sum(N*log(d)-d)-logepmin)
	return(dummy)}

# scan the parameter space
Atest=seq(0,3,by = 0.02)	
Btest=seq(0,3,by = 0.02)

jmax=length(Atest)
kmax=length(Btest)
j=1
k=1
# generate a posterior matrix for the contour plot
postplot = matrix(nrow=jmax,ncol=kmax)
while (j<=jmax){ k=1
	while (k<=kmax) {
	postplot[j,k]=posterior(Atest[j],Btest[k])
	k=k+1}
	j=j+1}
quartz()
contour(Atest,Btest,postplot,levels=c(0.9,0.7,0.5,0.3,0.1),drawlabels=FALSE,xlab='A',ylab='B',xlim=c(0,3),ylim=c(0,3))